Strip Astrometry Pipeline Reduction

7 March 2001ff
A. S. Bosh

Introduction

Package History

Revision Log

v2.4.1 (30 May 2002, asb).  Incorporate spRecycleFITSid to allow fitting more strips at once.

v2.4 (29 May 2002, asb).  Add option to input centroiding results from the astrometric star list, rather than centroiding them.  This allows one to attempt different methods of getting a good Pluto fit, without having to fit all the stars all over again.  Add removal of duds (stars that didn't fit) and output to a different file.

v2.3 (29 May 2002, asb).  Add option to check which operating system you're on, and add "../.." or not as appropriate.

v2.2 (25 May 2002, asb).  Add file output, extracting everything from a strip.

v2.1 (25 May 2002, asb).  Basic functionality present (except output of list), first conversion to a package.

v1.0 (7 March 2001, asb).  Begin this package, start it as a notebook and add package functionality later.

Overview

Inputs

Data file in SIMPLE=TRUE FITS format.  Image data.  However, for testing purposes, begin with SNAPSHOT SIMPLE=FALSE *.md file, converted to SIMPLE=TRUE FITS format.

Flat file in SIMPLE=TRUE FITS format.  Image data.

Dark file in SIMPLE=TRUE FITS format.  Image data.

Astrometric standard stars file in SIMPLE=TRUE FITS format, with TABLE extension.  Table data, giving star name/number, RA, Dec, magnitude, other quantities?

Occulting body ephemeris file in SIMPLE=TRUE FITS format, with TABLE extension.  Table data, giving date (MJD), RA, Dec (both geocentric of date), magnitude, other quantities?  Note:  the eventual plan is to replace this with a direct call to JPL SPICE routines to get the required position.

Outputs

Positions file in SIMPLE=TRUE FITS format, with TABLE extension.  Table data, giving star name/number, measured x and y coordinates and errors, fitted RA and Dec and errors, magnitude and error (other quantities)?  In addition, header keywords will include the astrometric solution parameters (RA offset, Dec offset, Dec plate scale, angle between RA and Dec axes).

Calculation Plan

Modification to analyze flanking diffraction images (when using bear-cage grating) to determine position of bright stars.

Questions

What machine will this be primarily run on, and how much memory does it have?  [If it has a lot of memory, I could read in the entire flat frame rather than doing partial reads.]

Is there supposed to be a low-light flat included in the procedure?  [There isn't one in flat1.1.c.]

Initialization (for testing/development only)

<< jleGroup`

This loads the jleGroupFastBinaryFilesLINK program.  Confirm that we have the correct version of ReadListBinary.  Its context should be FastBinaryFiles`.

Context[ReadListBinary]

Initial Package Stuff

Add to Options Package

BeginPackage["jleGroup`options`"];

Options for all functions

opPrintDiagnostics::usage = "Option for all sp functions, default \
is False.";

Options for spCentroid

opReturn::usage = "Option for spCentroid, default is \"CP\" \
(centers and photometry).  The other option is \"All\" to return all the \
results from centroiding (the format of this depends on the type of \
centroiding selected)";

opCentroidFunction::usage = "Option for spCentroid, default is \
\"margAnal\".  Other options are \"margGaussian\" and \"ellipticalLorentz\".  \
Quotes are required in all cases.";

Options for spMakeMeanRow

opOutputFile::usage = "Option for spMakeMeanRow and spExtractStar, \
default is Null.  To save the mean row or extracted star to a file, provide a \
filename with full pathname here.";

opRowStep::usage = "Option for spMakeMeanRow, default is 100.  \
Number of rows to read in and average at a time.";

opRowsRead::usage = "Option for spMakeMeanRow, default is All (no \
quotes).  Total number of rows to read from the image.";

Options for spFlatten

opFlatId::usage = "Option for spFlatten, default is Null.  To apply \
a flat, supply the FITS id of the file.";

opDarkId::usage = "Option for spFlatten, default is Null.  To apply \
a dark, supply the FITS id of the file.";

opOverscanPixels::usage = "Option for spFlatten, default is {}.  \
Not implemented.";

Options for spAnalyzeStrip

opPlanetName::usage = 
"Option for \
spAnalyzeStrip, default is \"Planet\".";

opPlanetName::usage = 
"Option for \
spAnalyzeStrip, default is \"Planet\".";

opPlanetNear::usage = 
"Option for \
spAnalyzeStrip, default is 250 pixels.  When finding the planet on the strip, \
the location predicted via combining the ra/dec from the ephemeris and the \
ra/dec plate scales is improved by finding a nearby star and offsetting from \
it.  This parameter sets how close that nearby star must be.  The function \
does not find the NEAREST star within this range, but instead just selects \
the FIRST star that meets this criterion.";

opDoAstrometryCentroiding::usage = 
"Option for \
spAnalyzeStrip, default is True.  When True, the spAnalyzeStrip function will \
read in the astrometry reference file and perform centroiding on all stars \
within it.  When False, the spAnalyzeStrip function will assume that the \
astrometry reference file is the result of prior centroiding on the data file \
that is passed.  When this option is passed as False, this allows testing of \
different ways of fitting the target stars and planet without needing to \
re-centroid all the stars in the astrometric star list.";

End

EndPackage[];

Begin the stripPipeline Package

Begin package and declare all other packages needed

BeginPackage[ "asbGroup`stripPipeline`",{
 
		"jleGroup`leastSquares`", 
		"jleGroup`options`",
		"jleGroup`fitsAide`",
		"jleGroup`imageProcessing`",
		"jleGroup`pointSources`",
		"jleGroup`timeAndCoords`"
		
		}];

Set lsIterateToConvergence options

SetOptions[lsIterateToConvergence, 
      opCheckImaginary -> True]

Public declarations

Usage statements

\(\(spMakeMeanRow::usage = 
"spMakeMeanRow[inputFilename, firstRow, options] reads a file and makes a mean row.  The input filename should contain the full pathname.  The firstRow parameter allows the user to skip over the ramp-up region of strips.  The function returns the 1-D mean row.  Options include opPrintDiagnostics, opOutputFile, opRowStep, and opRowsRead.  See usage statements for current defaults.";\)\)

\(\(spExtractStar::usage = 
"spExtractStar[inputFilename, firstRow, firstColumn, numberOfRows, numberOfColumns, options] extracts a segment around (presumably) a star, but there is no checking to ensure that a star is in the region.  The input filename should contain the full pathname.  Provide the first row, first column, number of rows, and number of columns for the extraction.  The function returns the FITS id for the extracted region.  Options include opPrintDiagnostics and opOutputFile.  See usage statements for current defaults.";\)\)

spFlatten::usage = 
"spFlatten[imageFITSid, \
options] applies dark and flat to an image or section of an image.  Biases \
are not currently implemented.  Provide the FITS id of the data section.  The \
function returns the 2-D flattened data array.  Options include \
opPrintDiagnostics, opFlatId, opDarkId, opOverscanPixels.  See usage \
statements for current defaults.  Note that opFlatId and/or opDarkId must be \
modified from their defaults in order for anything to happen.  Note also that \
the flattening is intended to be done on the extracted segments, not the \
whole frame at once, although it should be possible to do the whole \
frame.";

spCentroid::usage = 
"spCentroid[headerFITSid, \
dataArray, options] finds the center of the supplied image.  Provide the FITS \
id of the header and the 2-D data array.  The function returns {row center, \
column center, image total, row sigma, column sigma, image sigma}; depending \
on the model chosen the sigmas may be set to Null.  Options include \
opPrintDiagnostics, opReturn (this modifies the set of information that is \
returned by the function), opCentroidFunction.  See usage statements for \
current defaults.";

\(\(spAnalyzeStrip::usage = 
"spAnalyzeStrip[dataFileName, planetEphFileName, astromFile, targetFile, outFileBase, refStarNumber, plateScaleRowCol, options] determines the centers of all astrometry and target stars, and occulting planet on the strip, and stores these results to two files, outFileBase.astrom.dat (astrometric reference stars) and outFileBase.targ.dat (target stars and planet).  All file names must include the full pathname.  The reference star number is the index number of the star from astromFile to be used for initial centering of the strip.  The plateScaleRowCol must be in arcseconds per pixel.  Options include opPrintDiagnostics, opPlanetName, and opPlanetMagnitude.  See usage statements for current defaults.";\)\)

Options

\(\(Options[
          spAnalyzeStrip] = 
{opPrintDiagnostics -> 
            False, 
opPlanetName -> "Planet", 
opPlanetMagnitude -> 
            0.0, 
opPlanetNear -> 
            250, 
opDoAstrometryCentroiding -> 
            True};\)\)

Options[
          spMakeMeanRow] = 
{opOutputFile -> 
            Null,  
opPrintDiagnostics -> 
            False, 
opRowStep -> 
            100, 
opRowsRead -> All};

Options[
          spExtractStar] = 
{opOutputFile -> 
            Null,  
opPrintDiagnostics -> 
            False};

Options[
          spFlatten] = 
{opFlatId -> 
            Null, 
opDarkId -> 
            Null, 
opOverscanPixels -> {}, 
opPrintDiagnostics -> 
            False
};

\(\(Options[
          spCentroid] = 
{opCentroidFunction -> "margAnal", 
opPrintDiagnostics -> 
            False, 
opReturn -> "CP"};\)\)

lsModelContext

lsModelContext[spGaussian] = "`spGaussian`";

lsInitialNames

lsInitialNames[spGaussian] :={"x0","sigma","area","bkg"};

lsintialStep

Functions

Begin["`Private`"];

Note about naxis* symbols

It seems that jleGroup`fitsAide` makes symbols called naxis, naxis1, etc.  I was previously using these (defined in the Module) in spMakeMeanRow, but would always get blue messages about Context conflicts.  The symbols are created via faKeywordUse, and seem to be well embedded in the package so I didn't try to change anything.  However, in future major revisions to this package, we might consider revising this so that these "extra" symbols are not created on the sly.

spMakeMeanRow

Inputs:
  file (flat or dark)
  first row

Output:
  mean row

Options
  save mean row to file
  print diagnostics for debugging

Old Algorithm:
  Skip over the ramp-up rows until you get to flat first row (frow)
  Initialize sum row (mean[i]) to all 0s
  Loop over # of rows (nrow), adding each row to the sum row (mean[i])
  Divide the sum row (mean[i]) by the # of rows (flatnrow) to get the mean row (mean[i])
  Calculate mean signal level (mean_bar) of the mean row

New Algorithm: (*NO, this is too inefficient at disk access, go with old algorithm, except that we will get 1000 rows at a time*)
  Initialize mean row (mean[[i]]) to all 0s
  Loop over # of columns (NAXIS1):
    faPartialReadFITS, starting at row=first row, going for # of rows, getting one column
    take mean of that column, store in mean row (mean[[i]])
  Calculate mean signal level of the mean row

spMakeMeanRow[infile_,  frow_,  opts___]:=
Module[
{outfile, saveToFile, diag, nrowstep, 
            nrows, 
hid, na, na1, na2, 
            na3, 
meanRow, 
i, fList, 
            rList, did, nread, lastrow},

  
outfile = opOutputFile /. {opts} /. 
                Options[spMakeMeanRow];
saveToFile = 
              If[outfile === Null, False, True];
diag = opPrintDiagnostics /. {opts} /. 
                Options[spMakeMeanRow];
nrowstep = opRowStep /. {opts} /. 
                Options[spMakeMeanRow];
nrows = opRowsRead /. {opts} /. 
                Options[spMakeMeanRow];

  
  
\(If[
              FileType[infile] =!= 
                File, 
Print["Input frame ", 
                infile, " does not exist.  Exiting."]; 
Return[\(-1\)]]\);

  
hid = First[faReadHeaders[infile, 1]];
na = 
              faNumberValue[hid, "NAXIS"];
If[
              na =!= 2  && 
                na =!= 3, 
Print["Input frame ", 
                infile, " is not 2D or 3D.  (NAXIS=", 
                na, ").  Exiting."]; Return[-1]];
  
na1 = faNumberValue[hid, "NAXIS1"];
na2 = 
              faNumberValue[hid, "NAXIS2"];
If[
              na === 3, 
na3 = 
                faNumberValue[hid, "NAXIS3"]; 
If[
                na3 =!= 1, 
Print["NAXIS3=", 
                  na3, ", must be equal to 1.  Cannot do true 3D data.  \
Exiting."]; 
Return[-1]]];
\(If[
              diag, 
Print["Frame header loaded."]; 
Print["  NAXIS=", 
                na]; 
Print["  NAXIS1=", 
                na1]; 
Print["  NAXIS2=", 
                na2]; 
If[na === 3, 
                Print["  NAXIS3=", na3]]]\);

  
\(If[
              frow > na2, 
Print["First row (", 
                frow, ") is > number of rows in image (", 
                na2, ").  Exiting."]; 
Return[\(-1\)]]\);

  
meanRow = Table[0, {na1}];

  
lastrow = 
              If[nrows === All, na2, 
                If[frow + nrows - 1 > na2, na2, frow + nrows - 1]];
If[
              diag, 
Print["Entering loop over rows.  \
Will skip ", frow - 1, " rows to begin.  Will read through row ", 
                lastrow, "."]];
Table[
nread = 
                If[i + nrowstep - 1 > na2, na2 - i + 1, 
                  nrowstep]; 
If[
                na === 3, 
fList = {1, i, 1}; 
                rList = {na1, nread, 1}, 
fList = {1, i}; 
                rList = {na1, nread}]; 
did = 
                faPartialReadFITS[infile, 1, fList, 
                  rList, {}]; 
If[diag, 
                Print["fList=", fList, "\trList=", 
                  rList, "\tData ID=", did, "\tDimensions=", 
                  Dimensions[faData[did]]]]; 
If[
                na === 3, 
(meanRow = 
                        meanRow + #) &  /@  
                  N[faData[
                        did][[1]]], 
(meanRow = 
                        meanRow + #) &  /@  
                  N[faData[did]]]; 
spRecycleFITSid[
                did], 
{i, frow, lastrow, nrowstep}];

  
meanRow /= (na2 - frow + 1);

  
If[
              saveToFile, 
newHeader = 
                spMakeMeanHeader[infile, hid, 
                  na1]; 
newId = 
                faNewDataset[newHeader, meanRow]; 
If[
                diag, 
Print["Writing to FITS file ", 
                  outfile]]; 
faWriteFITS[outfile, 
                newId]];

  
meanRow
]

spMakeMeanHeader (*private*)

\(spMakeMeanHeader[infile_, inId_, 
        na1_] := 
Module[
{}, 

{{"SIMPLE", "T", ""}, 
{"BITPIX", "-64", ""}, 
{"NAXIS", "1", ""}, 
{"NAXIS1", 
            ToString[
              na1], ""}, 
{"BSCALE", "1", ""}, 
{"BZERO", "0", ""}, 
{"BUNIT", 
            faStringValue[
              inId, "BUNIT"], ""}, 
{"CRPIX1", 
            faStringValue[
              inId, "CRPIX1"], ""}, 
{"CDELT1", 
            faStringValue[
              inId, "CDELT1"], ""}, 
{"CRVAL1", 
            faStringValue[
              inId, "CRVAL1"], ""}, 
{"CFINT1", 
            faStringValue[
              inId, "CFINT1"], ""}, 
{"CLABL1", 
            faStringValue[
              inId, "CLABL1"], ""}, 
{"CUNIT1", 
            faStringValue[
              inId, "CUNIT1"], ""}, 
{"CDSCR1", 
            faStringValue[
              inId, "CDSCR1"], ""}, 
 (*CHANGE\ \(\(THIS!!\)!\)*) 
{"ORIGFILE", 
            If[StringLength[infile] > 60, StringTake[infile, \(-60\)], 
              infile], ""}, 
{"COMP1", $MachineName, ""}, 
{"PRGRM1", "Mathematica", ""}, 
{"DATE1", 
            ToString[
              Date[]], ""}, 
{"END", "", ""}}
]\)

(*Can later add package name and version by using built-ins, once they're set up*)

spRecycleFITSid (*private*)

spRecycleFITSid[
        id_] := 
jleGroup`fitsAide`Private`faId = 
        id

spExtractStar

Inputs:
  file (data)
  first row
  first column
  number of rows
  number of columns

Output:
  subframe

Options
  save subimage to file
  print diagnostics for debugging

Algorithm:
  Test that data file (input) exists
  Check for options (saving extracted image to a file, diagnostics)
  Perform faPartialReadFITS to get subimage
  Write subimage, if desired
  Return subimage
  (*See notes under spMakeSubHeader for info about CRVAL etc.*)

spExtractStar[infile_,  frow_,  fcol_, nrow_, ncol_, 
        opts___]:=
Module[
{outfile, saveToFile, diag, hid, na, na1, na2, na3, sfid, 
            newdata, newHeader, newId},

  
outfile = opOutputFile /. {opts} /. 
                Options[spExtractStar];
saveToFile = 
              If[outfile === Null, False, True];
diag = opPrintDiagnostics /. {opts} /. 
                Options[spExtractStar];

  
  
\(If[
              FileType[infile] =!= 
                File, 
Print["Input frame ", 
                infile, " does not exist.  Exiting."]; 
Return[\(-1\)]]\);

  
hid = First[faReadHeaders[infile, 1]];
na = 
              faNumberValue[hid, "NAXIS"];
If[
              na =!= 2  && 
                na =!= 3, 
Print["Input frame ", 
                infile, " is not 2D or 3D.  (NAXIS=", 
                na, ").  Exiting."]; Return[-1]];
  
na1 = faNumberValue[hid, "NAXIS1"];
na2 = 
              faNumberValue[hid, "NAXIS2"];
If[
              na === 3, 
na3 = 
                faNumberValue[hid, "NAXIS3"]; 
If[
                na3 =!= 1, 
Print["NAXIS3=", 
                  na3, ", must be equal to 1.  Cannot do true 3D data.  \
Exiting."]; 
Return[-1]]];
\(If[
              diag, 
Print["Frame header loaded."]; 
Print["  NAXIS=", 
                na]; 
Print["  NAXIS1=", 
                na1]; 
Print["  NAXIS2=", 
                na2]; 
If[na === 3, 
                Print["  NAXIS3=", na3]]]\);

  
\(If[
              frow > na2, 
Print["First row (", 
                frow, ") is > number of rows in image (", 
                na2, ").  Exiting."]; 
Return[\(-1\)]]\);
\(If[
              frow + nrow - 1 > 
                na2, 
Print["Last row (", 
                frow + nrow - 1, ") is > number of rows in image (", 
                na2, ").  Exiting."]; 
Return[\(-1\)]]\);
\(If[
              fcol > na1, 
Print["First column (", 
                fcol, ") is > number of columns in image (", 
                na1, ").  Exiting."]; 
Return[\(-1\)]]\);
\(If[
              fcol + ncol - 1 > 
                na1, 
Print["Last column (", 
                fcol + ncol - 1, ") is > number of columns in image (", 
                na1, ").  Exiting."]; 
Return[\(-1\)]]\);

  
If[
              na === 3, 
sfid = 
                faPartialReadFITS[infile, 
                  1, {fcol, frow, 1}, {ncol, nrow, 
                    1}, {1}]; 
newdata = faData[
                    sfid][[1]], 
  
sfid = 
                faPartialReadFITS[infile, 
                  1, {fcol, frow}, {ncol, 
                    nrow}, {1}]; 
newdata = faData[sfid]];

  
newHeader = 
              spMakeSubHeader[infile, hid, ncol, nrow, fcol, frow];
newId = 
              faNewDataset[newHeader, newdata];
If[
              saveToFile, 
If[
                diag, 
Print["Writing to FITS file ", 
                  outfile]]; 
faWriteFITS[outfile, 
                newId]];

  
newId
]

spMakeSubHeader (*private*)

\(spMakeSubHeader[infile_, inId_, ncol_, nrow_, fcol_, 
        frow_] := 
Module[
{}, 

{{"SIMPLE", "T", ""}, 
{"BITPIX", 
            faStringValue[
              inId, "BITPIX"], ""}, 
{"NAXIS", "2", ""}, 
{"NAXIS1", 
            ToString[ncol], ""}, 
{"NAXIS2", 
            ToString[nrow], ""}, 
{"BSCALE", 
            faStringValue[
              inId, "BSCALE"], ""}, 
{"BZERO", 
            faStringValue[
              inId, "BZERO"], ""}, 
{"BUNIT", 
            faStringValue[
              inId, "BUNIT"], ""}, 
{"CRPIX1", "1", ""}, 
{"CDELT1", "1.0", ""}, 
{"CRVAL1", 
            ToString[
              fcol], ""}, 
{"CFINT1", "1.0", ""}, 
{"CLABL1", "Column Number", ""}, 
{"CUNIT1", "pixel", ""}, 
{"CDSCR1", "on the CCD chip", ""}, 
{"CRPIX2", "1", ""}, 
{"CDELT2", "1.0", ""}, 
{"CRVAL2", 
            ToString[
              frow], ""}, 
{"CFINT2", "1.0", ""}, 
{"CLABL2", "Row Number", ""}, 
{"CUNIT2", "pixel", ""}, 
{"CDSCR2", "on the CCD chip", ""}, 
 (*CHANGE\ \(\(THIS!!\)!\)*) 
{"ORIGFILE", 
            If[StringLength[infile] > 60, StringTake[infile, \(-60\)], 
              infile], ""}, 
{"COMP1", $MachineName, ""}, 
{"PRGRM1", "Mathematica", ""}, 
{"DATE1", 
            ToString[
              Date[]], ""}, 
{"END", "", ""}}
]\)

(*Can later add package name and version by using built-ins, once they're set up*)

Header keyword information

CR keyword usage is described in the file gbMacOrUnixLocalFileName["astron.mit.edu","/pscf/jle-grp/software/lib/fits_templates/md_4.1.ft"].  Here is an excerpt:

spFlatten

Inputs:
  data FITS id

Output:
  flattened data block

Options
  FITS id for flat (default Null)
  FITS id for dark (default Null)
  print diagnostics for debugging

Algorithm:
  Hello

Note:  can't subtract overscans because we're reading in a subframe for the data, which isn't going to have the overscan pixels attached.  Difficult to do at the extraction level as well, because there we're doing a partial load at the position, so again the overscan pixels don't get read in.  Skip for now.

\(spFlatten[imageId_, 
        opts___] := 
Module[
{flatId, 
          darkId, overscan, diag, data, crval1, cdelt1, crpix1, ncol, fcol, 
          lcol, dark, flat, applyFlat, applyDark, 
          flbar}, 

 (*Get\ options*) 
applyFlat = False; 
applyDark = 
          False; 
flatId = \(opFlatId /. {opts}\) /. 
            Options[spFlatten]; 
darkId = \(opDarkId /. {opts}\) /. Options[
              spFlatten]; 
overscan = \(opOverscanPixels /. {opts}\) /. Options[
              spFlatten]; 
diag = \(opPrintDiagnostics /. {opts}\) /. Options[
              spFlatten]; 

 (*Get\ data*) 
data = 
          faData[imageId]; 

 (*Find\ first\ column*) 
crval1 = 
          faNumberValue[imageId, "CRVAL1"]; 
cdelt1 = 
          faNumberValue[imageId, "CDELT1"]; 
crpix1 = 
          faNumberValue[imageId, "CRPIX1"]; 
ncol = 
          faNumberValue[imageId, "NAXIS1"]; 
fcol = 
          crval1 + cdelt1*\((1 - crpix1)\); 
lcol = 
          fcol + ncol - 1; 
If[diag, 
          Print["First column in image is ", 
            fcol]; 
Print["Last column in image is ", 
            lcol]]; 

 (*Flattening\ subimage*) 
dark = 
          If[darkId =!= Null, applyDark = True; faData[darkId], 
            0]; 
flat = 
          If[flatId =!= Null, applyFlat = True; faData[flatId], 
            1]; 

 (*Subtract\ bias*) 
 (*\(If[
              overscan =!= {}, 
data = \(spSubtractOverscan[#, overscan] &\) /@ data; 
If[
                applyDark, 
                dark = spSubtractOverscan[dark, 
                    overscan]]; 
If[applyFlat, 
                flat = spSubtractOverscan[flat, 
                    
                    overscan]]
];\)*) 

 (*Subtract\ dark*) 
If[
          applyDark, 
 (*select\ appropriate\ columns*) 
dark = Take[dark, {fcol, lcol}]; 
If[
            diag, Print["Length of dark subline is ", 
              Length[dark]]; 
Print["Length of data line is ", Length[data[\([1]\)]]]]; 
data = \(\((# - 
                    dark)\) &\) /@ 
              data
]; 

 (*Divide\ by\ flat*) 
If[
          applyFlat, 
 (*select\ appropriate\ columns*) 
flat = 
            Take[flat, {fcol, lcol}] - dark; 
flbar = 
            Mean[flat]; 
If[diag, 
            Print["Length of flat subline is ", 
              Length[flat]]; 
Print["Length of data line is ", Length[data[\([1]\)]]]]; 
data = \(\((#/\((flat/
                        flbar)\))\) &\) /@ 
              data
]; 

 (*Return\ flattened\ data*) 
data
]\)

spSubtractOverscan (*private*)

spSubtractOverscan[row_, overscanPixels_, 
        opts___] := 
Module[
{bias}, 

bias = 
          Mean[row[[overscanPixels]]]; 

row - bias
]

spCentroid

Inputs:
  header FITS id
  data block (*separate from header id in case data have been flattened*)

Output:
  star row center
  star column center
  star total counts (background-subtracted)
  error in row center
  error in column center
  error in star total counts

Options
  print diagnostics for debugging

Algorithm:
  Hello

\(spCentroid[headerId_, data_, 
        opts___] := 
Module[
{centr, 
          diag}, 

 (*Get\ options*) 
centr = \(opCentroidFunction /. {opts}\) /. 
            Options[spCentroid]; 
diag = \(opPrintDiagnostics /. {opts}\) /. 
            Options[spCentroid]; 

Switch[
centr, 
"margAnal", 
          spMargAnal[headerId, 
            data, {opts}], 
"margGaussian", 
          spMargGaussian[headerId, 
            data, {opts}], 
"ellipticalLorentz", 
          spEllipticalLorentz[headerId, data, {opts}], 
_, 
          Print["Error, centroid function ", 
            centr, " not recognized.  Exiting."]; 
          Return[\(-1\)]]
]\)

spMargAnal (*private*)

If opReturn->"CP" (default), then returns {row center, column center, image total, Null, Null, Null}.  If opReturn->All, then returns:  { { {rowctr,   colctr}, geometric mean diameter,  image total, bkg per pixel},   {rowctr, row diameter, row image total, row  bkg per pixel},   {colctr, col diameter, col image total, col bkg per pixel}  }

\(spMargAnal[headerId_, data_, 
        opts___] := 
Module[
{diag, ret, 
          crval1, crval2, cdelt1, cdelt2, crpix1, crpix2, ncol, nrow, fcol, 
          frow, idHistory, box, frame, 
          results}, 

 (*Get\ options*) 
diag = \(opPrintDiagnostics /. Flatten[{opts}]\) /. 
            Options[spCentroid]; 
ret = \(opReturn /. 
              Flatten[{opts}]\) /. 
            Options[spCentroid]; 

 (*Find\ first\ column*) 
crval1 = 
          faNumberValue[headerId, "CRVAL1"]; 
cdelt1 = 
          faNumberValue[headerId, "CDELT1"]; 
crpix1 = 
          faNumberValue[headerId, "CRPIX1"]; 
ncol = 
          faNumberValue[headerId, "NAXIS1"]; 
fcol = 
          crval1 + cdelt1*\((1 - crpix1)\); 
If[diag, 
          Print["First column is ", fcol]; 
          Print["Number of columns is ", 
            ncol]]; 

 (*Find\ first\ row*) 
crval2 = 
          faNumberValue[headerId, "CRVAL2"]; 
cdelt2 = 
          faNumberValue[headerId, "CDELT2"]; 
crpix2 = 
          faNumberValue[headerId, "CRPIX2"]; 
nrow = 
          faNumberValue[headerId, "NAXIS2"]; 
frow = 
          crval2 + cdelt2*\((1 - crpix2)\); 
If[diag, 
          Print["First row is ", frow]; 
          Print["Number of rows is ", 
            nrow]]; 

 (*Create\ "frame"\ object\ for\ ipMargAnal*) 
idHistory = {}; 
box = {{frow, fcol}, {nrow, 
              ncol}}; 
frame = {idHistory, box, 
            data}; 

 (*Perform\ ipMargAnal*) 
results = 
          N[ipMargAnal[
              frame]]; 

 (*Return\ center, \ image\ total, \ and\ errors*) 
If[
          
          ret === "All", 
results, 
{results[\([1, 1, 
                1]\)] (*row\ center*) , 
results[\([1, 1, 
                2]\)] (*column\ center*) , 
results[\([1, 
                3]\)] (*image\ total*) , 
Null (*error\ in\ row\ center*) , 
Null (*error\ in\ column\ center*) , 
Null (*error\ in\ image\ total*) 
}]
]\)

spMargGaussian (*private*)

If opReturn->"CP" (default), then returns {row center, column center, image total, row sigma, column sigma, total sigma}.  If opReturn->All, then returns {{row center, row diameter, row total, row bkg}, {column center, column diameter, column total, column bkg}}

spaussian[p_, 
        c_] := 
Module[
{x0, sigma, 
          area}, 

{x0, sigma, area, 
            bkg} = p; 

N[(area/(sigma Sqrt[2  Pi])) 
                  Exp[-0.5 ((# - x0)/sigma)^2] + bkg] & /@ 
          c
]

\(spMargGaussian[headerId_, data_, 
        opts___] := 
Module[
{diag, ret, 
          disp, crval1, crval2, cdelt1, cdelt2, crpix1, crpix2, ncol, nrow, 
          fcol, frow, mcol, mrow, res, ctr, rowres, 
          colres}, 

 (*Get\ options*) 
diag = \(opPrintDiagnostics /. Flatten[{opts}]\) /. 
            Options[spCentroid]; 
ret = \(opReturn /. 
              Flatten[{opts}]\) /. 
            Options[spCentroid]; 

 (*Find\ first\ column*) 
crval1 = 
          faNumberValue[headerId, "CRVAL1"]; 
cdelt1 = 
          faNumberValue[headerId, "CDELT1"]; 
crpix1 = 
          faNumberValue[headerId, "CRPIX1"]; 
ncol = 
          faNumberValue[headerId, "NAXIS1"]; 
fcol = 
          crval1 + cdelt1*\((1 - crpix1)\); 
If[diag, 
          Print["First column is ", fcol]; 
          Print["Number of columns is ", 
            ncol]]; 

 (*Find\ first\ row*) 
crval2 = 
          faNumberValue[headerId, "CRVAL2"]; 
cdelt2 = 
          faNumberValue[headerId, "CDELT2"]; 
crpix2 = 
          faNumberValue[headerId, "CRPIX2"]; 
nrow = 
          faNumberValue[headerId, "NAXIS2"]; 
frow = 
          crval2 + cdelt2*\((1 - crpix2)\); 
If[diag, 
          Print["First row is ", frow]; 
          Print["Number of rows is ", 
            nrow]]; 

 (*Create\ marginal\ column\ and\ row*) 
mcol = \(Apply[Plus, #] &\) /@ 
            data; 
mrow = \(Apply[Plus, #] &\) /@ 
            Transpose[
              data]; 

 (*Get\ approximate\ center*) 
res = 
          spMargAnal[headerId, data]; 
ctr = 
          res[\([{1, 2}]\)]; \  (*{row, \ 
            column}*) 

 (*Fit\ marginal\ distributions\ for\ Gaussian\ parameters*) 
lsInitialNames[
            spGaussian] = {"x0", "sigma", "area", "bkg"}; 
lsUseModel[spGaussian]; 
disp = 
          If[diag, True, 
            False]; 

 (*row\ center*) 
lsInitialValues = {ctr[\([1]\)], 2. , res[\([3]\)], 
            mcol[\([1]\)]}; 
lsParamHist = {lsParams = 
              lsInitialValues}; 
{lsCoor, 
            lsData} = {Range[frow, frow + nrow - 1], 
            mcol}; 
lsFitList = {1, 2, 3, 
            4}; 
If[diag, 
          Print["Fitting for row center, Gaussian parameters:  ", 
            lsNames[\([lsFitList]\)]]]; 
lsIterateToConvergence[20, 10^\(-6\), lsFitList, 
          opDisplay -> 
            disp]; 
rowres = {lsParamHist[\([\(-1\)]\)], 
            lsError}; 

 (*column\ center*) 
lsInitialValues = {ctr[\([2]\)], 2. , res[\([3]\)], 
            mrow[\([1]\)]}; 
lsParamHist = {lsParams = 
              lsInitialValues}; 
{lsCoor, 
            lsData} = {Range[fcol, fcol + ncol - 1], 
            mrow}; 
lsFitList = {1, 2, 3, 
            4}; 
If[diag, 
          Print["Fitting for column center, Gaussian parameters:  ", 
            lsNames[\([lsFitList]\)]]]; 
lsIterateToConvergence[20, 10^\(-6\), lsFitList, 
          opDisplay -> 
            disp]; 
colres = {lsParamHist[\([\(-1\)]\)], 
            lsError}; 

If[
          ret === "All", 
{rowres, 
            colres}, 
{rowres[\([1, 1]\)], 
            colres[\([1, 1]\)], \((rowres[\([1, 3]\)] + colres[\([1, 3]\)])\)/
              2, 
rowres[\([2, 1]\)], colres[\([2, 1]\)], 
            Sqrt[rowres[\([2, 3]\)]^2 + colres[\([2, 3]\)]^2]/
              2}
]
]\)

spEllipticalLorentz (*private*)

If opReturn->"CP" (default), then returns {row center, column center, image total, row sigma, column sigma, Null}.  If opReturn->All, then returns {{bkg, peak, row center, column center, diameter, shapeIndex, axisAngle, axisRatio}

\(spEllipticalLorentz[headerId_, data_, 
        opts___] := 
Module[
{diag, ret, 
          disp, crval1, crval2, cdelt1, cdelt2, crpix1, crpix2, ncol, nrow, 
          fcol, frow, res, peakr, peakc, peak, bkgr, bkgc, bkg, diamr, diamc, 
          diam, axisRatio, row, col, results, rtnCode, iteration, 
          maxIterations, 
          maxFracParamError}, 

 (*Get\ options*) 
diag = \(opPrintDiagnostics /. 
              Flatten[{opts}]\) /. 
            Options[spCentroid]; 
ret = \(opReturn /. 
              Flatten[{opts}]\) /. 
            Options[spCentroid]; 

 (*Find\ first\ column*) 
crval1 = 
          faNumberValue[headerId, "CRVAL1"]; 
cdelt1 = 
          faNumberValue[headerId, "CDELT1"]; 
crpix1 = 
          faNumberValue[headerId, "CRPIX1"]; 
ncol = 
          faNumberValue[headerId, "NAXIS1"]; 
fcol = 
          crval1 + cdelt1*\((1 - crpix1)\); 
If[diag, 
          Print["First column is ", fcol]; 
          Print["Number of columns is ", 
            ncol]]; 

 (*Find\ first\ row*) 
crval2 = 
          faNumberValue[headerId, "CRVAL2"]; 
cdelt2 = 
          faNumberValue[headerId, "CDELT2"]; 
crpix2 = 
          faNumberValue[headerId, "CRPIX2"]; 
nrow = 
          faNumberValue[headerId, "NAXIS2"]; 
frow = 
          crval2 + cdelt2*\((1 - crpix2)\); 
If[diag, 
          Print["First row is ", frow]; 
          Print["Number of rows is ", 
            nrow]]; 

 (*Get\ approximate\ parameters*) 
If[diag, 
          Print["Getting marginal analysis results for initial parameter values..."]]; 
res = 
          spMargAnal[headerId, data, 
            opReturn -> "All"]; 
bkg = 
          res[\([1, 4]\)]; 
diam = 
          res[\([1, 2]\)]; 
{diamr, diamc} = 
          res[\([{2, 3}, 2]\)]; 
axisRatio = 
          If[diamr > diamc, diamr/diamc, 
            diamc/diamr]; 
{row, col} = 
          res[\([1, 1]\)]; 
total = 
          res[\([1, 3]\)]; 
peak = 
          total/\((diam*
                Sqrt[2  
                    Pi])\); 

 (*Fit\ for\ elliptical\ Lorentz\ psf*) 
lsNames = 
          lsInitialNames[
            psSingleSource[psConstantBackground, 
              psEllipticalLorentzStar]]; 
lsUseModel[
          psSingleSource[psConstantBackground, 
            psEllipticalLorentzStar]]; 
disp = 
          If[diag, True, 
            False]; 

 (*Set\ up\ for\ fitting*) 
lsInitialValues = 
          N[{bkg, peak, row, col, diam, 2, 0, 
              axisRatio}]; 
If[diag, 
          Print["Initial values = ", 
            lsInitialValues]]; 
lsParamHist = {lsParams = 
              lsInitialValues}; 
lsChiSquaredHist = {}; 
lsChangeFrac = 1.0; 
lsWeight = 
          1; 
lsUsePoints = 1; 
lsData = 
          Flatten[data]; 
        lsCoor = 
          Flatten[Table[{i, j}, {i, frow, frow + nrow - 1}, {j, fcol, 
                fcol + ncol - 1}], 1]; 
lsFitList = {1, 2, 
            3, 4, 5, 6, 7, 
            8}; 

 (*Fitting*) 
If[diag, 
          Print["Fitting for all elliptical Lorentz parameters:  ", 
            lsNames[\([lsFitList]\)]]]; 
 (*\(rtnCode = 
              lsIterateToConvergence[10, 10^\(-3\), lsFitList, 
                opDisplay -> disp];\)*) 
iteration = 
          1; 
maxFracParamError = 
          10^\(-3\); 
maxIterations = 
          10; 
For[
Null, 
iteration == 1 || 
            maxFracError > maxFracParamError && 
              iteration ≤ 
                maxIterations, 
\(iteration++\), 
lsUseLastParams; 

 (*Checking\ that\ parameters\ are\ reasonable*) 
If[
            MemberQ[\(\((# < 0)\) &\) /@ 
                lsParams[\([{1, 2, 3, 4, 5, 6, 8}]\)], True], 
            Print["Error in lsIterateToConvergence, parameters are non-physical.  Returning."]; 
rtnCode = \(-1\); 
            iteration = 
              maxIterations;
]; 

lsIterate; 
lsUpdateStep[
            0.005]; 

 (*Checking\ that\ fit\ is\ working\ ok*) 
If[Im[lsChiSquared] =!= 0, 
            Print["Error in lsIterateToConvergence, function derivatives are imaginary.  Returning."]; 
rtnCode = \(-1\); 
            iteration = 
              maxIterations;
]; 

If[
            disp, 
Print[
              CForm[lsChiSquared], " ", \((N[#1, 
                      Max[4, Ceiling[
                          Log[10, 
                            If[Abs[#1] == 0, $MachineEpsilon, 
                              Abs[#1]]]]]] &)\) /@ 
                lsParams[[lsFitList]]]
]
]; 

 (*Return\ result*) 
results = 
          If[rtnCode === \(-1\), 
 (*then*) {{0, 0, 0, 0, 
                0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
                0}}, 
 (*else*) {lsParamHist[\([\(-1\)]\)], 
              lsError}]; 
phot = 
          If[rtnCode === \(-1\), 
 (*then*) 
              0, 
 (*else*) Apply[Plus, 
              lsModel - results[\([1, 1]\)]]]; 
If[
          ret === "All", 
results, 
Flatten[{results[\([1, {3, 4}]\)], phot, 
              results[\([2, {3, 4}]\)], 
              Null}]
]
]\)

maxFracError := 
Module[
{index, 
          fracErrorList = {}}, 

Do[
          index = lsFitList[[
              j]]; AppendTo[fracErrorList, 
              Abs[((lsParamHist[[-1, index]] - lsParamHist[[-2, index]])/jleGroup)`leastSquares`lsError[[j]]]];, {j, Length[lsFitList]}]; 
        Max[fracErrorList]]

spAnalyzeStrip

\(spAnalyzeStrip[dataFileName_, planetEphFileName_, astromFile_, 
        targetFile_, outFileBase_, refStarNumber_, plateScaleRowCol_, 
        opts___] := 
Module[
{diag, 
          planetName, planetMag, rc, id, ctr, delta, i, prefix, planetNear, 
          doAstromCentr, 
hid, ymd, ut, dateline, timeline, 
          outtable, refNear, res, 
astromTable, 
          astromResults, astromStrings, 
astromTable2, 
          astromResults2, astromStrings2, 
targetTable, 
          targetResults, targetStrings, 
targetTable2, 
          targetResults2, targetStrings2, 
planetTable, 
          planetResults, planetStrings, 
planetTable2, 
          planetResults2, planetStrings2, 
planetRADec, 
          refRADec, refRowColumn, 
          planetRowColumn
}, 

 (*Get\ options*) 
diag = \(opPrintDiagnostics /. {opts}\) /. 
            Options[spAnalyzeStrip]; 
planetName = \(opPlanetName /. {opts}\) /. 
            Options[spAnalyzeStrip]; 
planetMag = \(opPlanetMagnitude /. {opts}\) /. 
            Options[spAnalyzeStrip]; 
planetNear = \(opPlanetNear /. {opts}\) /. 
            Options[spAnalyzeStrip]; 
doAstromCentr = \(opDoAstrometryCentroiding /. {opts}\) /. 
            Options[spAnalyzeStrip]; 

 (*Add\ annoying\ pathname\ prefix\ if\ running\ Mac\ OSX . \ \ For\ some\ reason\ this\ is\ necessary\ for\ anything\ read\ in\ via\ FastBinaryFiles*) 
prefix = 
          If[$MachineType === "PowerMac" && $OperatingSystem === "MacOSX", 
"../..", 
""]; 

 (*Read\ in\ reference\ file\ of\ astrometry\ star\ positions*) 
If[diag, 
          Print["Reading in reference file of astrometry star positions"]]; 
astromTable = 
          spReadFititFile[
            astromFile]; 

 (*Determine\ offset\ for\ new\ strip*) 
If[diag, 
          Print["Determining row, column offsets for new strip"]]; 
rc = 
          Round[astromTable[\([refStarNumber, {2, 4}]\)] - 
              25]; 
id = 
          spExtractStar[prefix <> dataFileName, rc[\([1]\)], rc[\([2]\)], 50, 
            50]; 
ctr = 
          spCentroid[id, faData[id]]; 
If[diag, 
          ipShowFrame[{1, {rc, {50, 50}}, 
              faData[id]} (*\(,\)\(opBox -> {Round[ctr] - 1, {2, 
                    2}}\)*) ]]; 

delta = 
          astromTable[\([refStarNumber, {2, 4}]\)] - 
            ctr[\([{1, 2}]\)]; 
If[diag, 
          Print["	Offset (reference-current) = ", 
            delta, "{rows, columns}"]]; 

 (*Get\ new\ astrometry\ star\ centers*) 
If[
          doAstromCentr, 
If[diag, 
            Print["Getting new centers for astrometry stars"]]; 
len = 
            Length[astromTable]; 
astromResults = 
Table[
rc = 
                Round[astromTable[\([i, {2, 4}]\)] - delta - 
                    20]; 
If[diag, 
                Print["	Centroiding #", i, " of ", 
                  len, ".  Extract at rc = ", 
                  rc, " for {40,40}"]]; 
id = 
                spExtractStar[prefix <> dataFileName, rc[\([1]\)], 
                  rc[\([2]\)], 40, 40]; 
res = 
                spCentroid[id, faData[id], 
                  opCentroidFunction -> "ellipticalLorentz"]; 
If[
                MemberQ[NumberQ /@ Drop[res, \(-1\)], False] || 
                  MemberQ[\(\((Head[#] === Complex)\) &\) /@ res, True], 
                res = {0, 0, 0, 0, 0, 0}]; 
spRecycleFITSid[
                id]; 
res, 
{i, 
                len}
]; 

 (*Write\ out\ astrometry\ star\ results*) 
If[diag, 
            Print["Writing out astrometry star results"]]; 
hid = 
            faReadHeaders[prefix <> dataFileName, 
              1]; 
ymd = 
            faStringValue[
              hid[\([1]\)], "DATE-OBS"]; 
ut = 
            faStringValue[
              hid[\([1]\)], "UT"]; 
dateline = 
            StringReplace["DATADATE= "" <> 
                ymd <> """, "-" -> "/"]; 
timeline = "DATATIME= " <> 
              ut; 
astromStrings = 
            spOutputStrings[astromTable, 
              astromResults]; 
outtable = 
            Flatten[{dateline <> "
", timeline <> "
", 
                astromStrings}]; 
spWriteFititFile[
            outFileBase <> ".astrom.dat", 
            outtable]; 

 (*Write\ out\ dud - free\ astrometry\ star\ results*) 
{astromTable2, 
              astromResults2} = 
            spRemoveDuds[astromTable, 
              astromResults]; 
astromStrings2 = 
            spOutputStrings[astromTable2, 
              astromResults2]; 
outtable = 
            Flatten[{dateline <> "
", timeline <> "
", 
                astromStrings2}]; 
spWriteFititFile[
            outFileBase <> ".use", 
            outtable]
]; 

If[\(! doAstromCentr\), 
If[diag, 
            Print["Skipping centroiding of astrometric reference stars, going to target stars."]]; 
hid = 
            faReadHeaders[prefix <> dataFileName, 
              1]; 
ymd = 
            faStringValue[
              hid[\([1]\)], "DATE-OBS"]; 
ut = 
            faStringValue[
              hid[\([1]\)], "UT"]; 
dateline = 
            StringReplace["DATADATE= "" <> 
                ymd <> """, "-" -> "/"]; 
timeline = "DATATIME= " <> 
              ut;
]; 

 (*Read\ in\ reference\ file\ of\ target\ star\ positions*) 
If[diag, 
          
          Print["Reading in reference file of target star positions"]]; 
targetTable = 
          spReadFititFile[
            targetFile]; 
 (*\(If[diag, 
              Print[targetTable]];\)*) 

 (*Get\ new\ target\ star\ centers*) 
If[diag, 
          Print["Getting new centers for target stars"]]; 
len = 
          Length[targetTable]; 
targetResults = 
Table[
rc = 
              Round[targetTable[\([i, {2, 4}]\)] - delta - 
                  20]; 
If[diag, 
              Print["	Centroiding #", i, " of ", 
                len, ".  Extract at rc = ", 
                rc, " for {40,40}"]]; 
id = 
              spExtractStar[prefix <> dataFileName, rc[\([1]\)], rc[\([2]\)], 
                40, 40]; 
res = 
              spCentroid[id, faData[id], 
                opCentroidFunction -> "ellipticalLorentz"]; 
If[
              MemberQ[NumberQ /@ Drop[res, \(-1\)], False] || 
                MemberQ[\(\((Head[#] === Complex)\) &\) /@ res, True], 
              res = {0, 0, 0, 0, 0, 0}]; 
spRecycleFITSid[
              id]; 
res, 
{i, 
              len}
]; 
 (*\(If[diag, 
              Print[targetResults]];\)*) 

 (*Get\ planet\ position*) 
If[diag, 
          Print["Getting planet position information"]]; 
planetRADec = 
          N[spGetPlanetPosition[planetEphFileName, ymd, ut]/
              Degree]; 
refRADec = 
          astromTable[\([refStarNumber, {9, 
                10}]\)]; 
refRowColumn = 
          ctr[\([{1, 2}]\)]; 
planetRowColumn = 
          N[\((planetRADec - refRADec)\)*3600/plateScaleRowCol] + 
            refRowColumn; 
If[diag, 
          Print["	planetRowColumn=", 
            planetRowColumn]]; 
 (*Find\ nearer\ reference\ star, \ to\ get\ more\ accurate\ center\ for\ planet*) 
refNear = 
          Null; \(If[
              Abs[planetRowColumn[\([1]\)] - #[\([2]\)]] < planetNear && 
                Abs[planetRowColumn[\([2]\)] - #[\([4]\)]] < 
                  planetNear, 
refNear = #] &\) /@ 
          astromTable; 
 (*If\ it\ did\ not\ find\ a\ nearby\ star ... *) 
If[
          refNear === Null, 
If[diag, 
            Print["Did not find a nearby star, using default extraction center"]]; 
          refNear = 
            astromTable[\([refStarNumber]\)]]; 
refRADec = 
          refNear[\([{9, 10}]\)]; 
refRowColumn = 
          refNear[\([{2, 4}]\)] - delta; 
planetRowColumn = 
          N[\((planetRADec - refRADec)\)*3600/plateScaleRowCol] + 
            refRowColumn; 
If[diag, 
          Print["	planetRowColumn=", 
            planetRowColumn]]; 
planetTable = {{9000, 
              planetRowColumn[\([1]\)], 0, planetRowColumn[\([2]\)], 0, 0, 0, 
              planetName, planetRADec[\([1]\)], planetRADec[\([2]\)], 
              planetMag, 0}}; 
 (*\(If[diag, 
              Print[planetTable]];\)*) 

 (*Get\ planet\ center*) 
rc = 
          Round[planetRowColumn - 10]; 
If[diag, 
          Print["Getting planet center from strip.  Extract at rc = ", 
            rc, " for {20,20}"]]; 
id = 
          spExtractStar[prefix <> dataFileName, rc[\([1]\)], rc[\([2]\)], 20, 
            20]; 
res = 
          spCentroid[id, faData[id], 
            opCentroidFunction -> "ellipticalLorentz"]; 
If[
          MemberQ[NumberQ /@ Drop[res, \(-1\)], False] || 
            MemberQ[\(\((Head[#] === Complex)\) &\) /@ res, True], 
          res = {0, 0, 0, 0, 0, 
              0}]; 
planetResults = {res}; 
If[diag, 
          Print["	Fitted planet center is:  ", 
            planetResults]]; 

 (*Write\ out\ target\ star\ and\ planet\ results*) 
If[diag, 
          Print["Writing out results from target star(s) and planet"]]; 
targetStrings = 
          spOutputStrings[targetTable, 
            targetResults]; 
planetStrings = 
          spOutputStrings[planetTable, 
            planetResults]; 
outtable = 
          Flatten[{dateline <> "
", timeline <> "
", targetStrings, 
              planetStrings}]; 
spWriteFititFile[
          outFileBase <> ".targ.dat", 
          outtable]; 

 (*Write\ out\ dud - free\ planet\ and\ target\ star\ results*) 
{targetTable2, targetResults2} = 
          spRemoveDuds[Join[targetTable, planetTable], 
            Join[targetResults, 
              planetResults]]; 
targetStrings2 = 
          spOutputStrings[targetTable2, 
            targetResults2]; 
outtable = 
          Flatten[{dateline <> "
", timeline <> "
", 
              targetStrings2}]; 
spWriteFititFile[
          outFileBase <> ".stars", 
          outtable]; 

 (*Return*) 
{astromResults, targetResults, 
          planetResults}
]\)

spGetPlanetPosition (*private*)

\(spGetPlanetPosition[jplEphFile_, ymd_, 
        ut_] := 
Module[
{obsmjd, pid, 
          pdata, mjd, radecrange, ra, dec, rai, 
          deci}, 

 (*Start\ time\ of\ strip*) 
obsmjd = 
          tcDateTimeStringToMJD[
            StringReplace[
              ymd <> ut, {"-" -> " ", ":" -> " "}]]; 

 (*Read\ in\ ephemeris*) 
pid = 
          faColumnsInput[jplEphFile]; 
pdata = 
          faData[pid[\([\)\(1\)\(]\)]]; 

 (*Make\ ra\ and\ dec\ interpolation\ functions*) 
mjd = \(Transpose[
              pdata]\)[\([\)\(1\)\(]\)]; 
radecrange = \(tcRectangularToSpherical[#[\([\)\({2, 3, 4}\)\(]\)]] &\) /@ 
            pdata; 
{ra, \ 
            dec} = \(Transpose[radecrange]\)[\([\)\({1, 
              2}\)\(]\)]; 
rai = 
          Interpolation[Transpose[{mjd, ra}]]; 
deci = 
          Interpolation[
            Transpose[{mjd, 
                dec}]]; 

 (*Return\ ra\ and\ dec\ at\ observed\ time*) 
N[{rai[obsmjd], 
            deci[obsmjd]}]
]\)

spOutputStrings (*private*)

The "full" list is read in from the file of standard stars, starting centers.  The "new" list is the fitted centers on the current frame.  The format of "full" is {star number, row center, row sigma, column center, column sigma, instrumental magnitude, mag sigma (always 0), star name, ra (deg), dec (deg), catalog magnitude, mean catalog position error}.  The format of "new" is {row center, column center, total signal, row sigma, column sigma, signal sigma (always Null)}.  Both "full" and "new" are 2D lists

\(spOutputStrings[full_, \ 
        new_] := 
Module[
{tfull, tnew, 
          s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, 
          stable}, 

If[
          Length[full] =!= 
            Length[new], 
Print["Error, the two input lists must be the same length.  ", Length[full], " =!= ", 
            Length[new], ".  Exiting."]; 
Return[\(-1\)]]; 

tfull = 
          Transpose[full]; 
tnew = 
          Transpose[
            new]; 

 (*star\ number*) 
s1 = \(toStringInt[#, 5] &\) /@ 
            tfull[\([1]\)]; 

 (*row\ center*) 
s2 = \(toString[#, {8, 3}] &\) /@ 
            tnew[\([1]\)]; 

 (*row\ error*) 
s3 = \(toString[#, {4, 3}] &\) /@ 
            tnew[\([4]\)]; 

 (*column\ center*) 
s4 = \(toString[#, {8, 3}] &\) /@ 
            tnew[\([2]\)]; 

 (*column\ error*) 
s5 = \(toString[#, {4, 3}] &\) /@ 
            tnew[\([5]\)]; 

 (*instrumental\ magnitude; \ 
          0\ total\ signal\ is\ an\ error\ flag*) 
s6 = 
\(If[# === 0, 
                toString[
                  0, {6, 3}], 
toString[\(-2.5\)*
                    Log[10, #], {6, 3}]] &\) /@ 
            tnew[\([3]\)]; 

 (*magnitude\ error*) 
s7 = \(toString[#, {5, 3}] &\) /@ 
            tfull[\([7]\)]; 

 (*star\ name*) 
s8 = \(toStringString[#, 21] &\) /@ 
            tfull[\([8]\)]; 

 (*catalog\ RA, \ what\ epoch\ should\ it\ \(\(be\)\(?\)\)*) 
s9 = \(toString[#, {12, 9}] &\) /@ 
            tfull[\([9]\)]; 

 (*catalog\ Dec, \ what\ epoch\ should\ it\ \(\(be\)\(?\)\)*) 
s10 = \(toString[#, {13, 9}] &\) /@ 
            tfull[\([10]\)]; 

 (*catalog\ magnitude*) 
s11 = \(toString[#, {3, 1}] &\) /@ 
            tfull[\([11]\)]; 

 (*average\ position\ error*) 
s12 = \(toString[#, {4, 3}] &\) /@ 
            tfull[\([12]\)]; 

 (*put\ them\ all\ together*) 
stable = 
          Transpose[{s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, 
              s12}]; 

\(\((#[\([1]\)] <> " " <> #[\([2]\)] <> " " <> #[\([3]\)] <> " " <> #[\([4]\)] <> " " <> #[\([5]\)] <> "  " <> #[\([6]\)] <> "  " <> #[\([7]\)] <> "    " <> #[\([8]\)] <> #[\([9]\)] <> " " <> #[\([10]\)] <> "   " <> #[\([11]\)] <> "      " <> #[\([12]\)] <> "
")\) &\) /@ 
          stable
]\)

\(toString[
        expr_, {n_, f_}] := 
Module[
{s, 
          len}, 
s = 
          ToString[
            NumberForm[N[expr], {n, f}, 
              NumberPadding -> {"", "0"}, 
              ExponentFunction -> \((Null &)\)]]; 
len = StringLength[s]; 
If[\((len - 1)\) < n, 
          StringJoin[Table[" ", {n - \((len - 1)\)}], s], 
          s]
]\)

\(toStringInt[expr_, 
        n_] := 
Module[
{s, 
          len}, 
s = 
          ToString[
            NumberForm[expr, n, NumberPadding -> {"", "0"}, 
              ExponentFunction -> \((Null &)\)]]; 
len = StringLength[s]; 
If[len < n, 
          StringJoin[Table[" ", {n - len}], s], 
          s]
]\)

\(toStringString[expr_, 
        n_] := 
Module[
{len}, 
len = StringLength[expr]; 
If[len < n, 
          StringJoin[expr, Table[" ", {n - len}]], 
          expr]
]\)

spReadFititFile (*private*)

\(spReadFititFile[
        fileName_] := 
Module[
{fp, 
          table}, 

fp = 
          OpenRead[fileName]; 
Skip[fp, 
          String]; 
Skip[fp, 
          String]; 
table = 
          Import[fp, "Table"]; 
Close[
          fp]; 

table
]\)

spWriteFititFile (*private*)

spWriteFititFile[outputFileName_, outputTable_, 
        opts___] := 
Module[
{fp}, 

fp = 
          OpenWrite[
            outputFileName]; 
WriteString[fp, #] & /@ 
          outputTable; 
Close[
          fp];
]

spRemoveDuds (*private*)

If the star didn't fit correctly, its results entry will be {0,0,0,0,0,0}

spRemoveDuds[table_, results_, 
        opts___] := 
Module[
{newtable, 
          newresults}, 

 (*Initialize new,  dud - 
            free table and results*) 
newtable = {}; 
newresults = {}; 

 (*Check each entry*) 
If[
              results[[#, {1, 2}]] =!= {0, 
                  0}, 
AppendTo[newtable, 
                table[[#]]]; 
AppendTo[newresults, 
                results[[#]]];
] & /@ 
          Range[Length[
              results]]; 

 (*Return dud - free table*) 
{newtable, 
          newresults}
]

End the Package

End[(* "`Private`" *)];

Protect["asbGroup`stripPipeline`*"];

EndPackage[];

Testing procedure (*not current*)

Make mean flat row

flatfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu", "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:Strip_astrometry:Converted:093053.a1.fits"];

flatrow = 
          spMakeMeanRow[flatfile,  1000, opPrintDiagnostics -> True, 
            opRowStep -> 500, 
            opOutputFile -> "ursula:Local Files:Projects:Strip astrometry \
pipeline (EWD):Test data:meanFlatRow20010621.fits"]; // Timing

ListPlot[flatrow, PlotJoined -> True, 
        PlotRange -> All];

Make mean dark row

darkfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu", "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:Strip_astrometry:Converted:091319.a1.fits"];

darkrow = 
          spMakeMeanRow[darkfile,  1000, opPrintDiagnostics -> True, 
            opRowStep -> 500, 
            opOutputFile -> "ursula:Local Files:Projects:Strip astrometry \
pipeline (EWD):Test data:meanDarkRow20010621.fits"]; // Timing

ListPlot[darkrow, PlotJoined -> True, 
        PlotRange -> All];

Extract a star subimage

datafile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu", "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:Strip_astrometry:Converted:083418.a1.fits"];

starfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu",  "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:star20010621.fits"];

id = 
      spExtractStar[datafile,  1050,  1165, 50, 50, 
        opPrintDiagnostics -> True]

id = 
      spExtractStar[datafile,  1050,  1165, 50, 50, 
        opOutputFile -> starfile, 
        opPrintDiagnostics -> True]

faHeader[id] // ColumnForm

ListDensityPlot[faData[id], Mesh -> False];

Test reading back in, make sure everything is the same...

id2 = faReadFITS[starfile]

faHeader[id2[[1]]] // ColumnForm

faHeader[id2[[1]]] == faHeader[id]

ListDensityPlot[faData[id2[[1]]], 
        Mesh -> False];

faData[id2[[1]]] == faData[id]

OK, everything is the same!

Notes about flattening

Attempts to flatten data:  the raw data and darks have a diagonal moire pattern that is not visible on the flat image (although it may be there).  In addition, the flat has bad columns (signal dropouts at columns 109, 484, 987, 1396) that are not present in the data; perhaps this is due to some effect that only occurs over a certain pedestal value.  I attempted to make low-light flats by taking the median of a subimage of 200 rows (all 2060 columns), but this has faintly dark "stripes" that correlate with stars in the image.  For now, work into the function the ability to not flatten or to flatten.

The anomalous up-spike in the mean dark row is due to a bright star (?) at row ~2885, column~1550.  So these aren't really darks?

In version 1.1 of this notebook, there are various plots/images of the flat and dark and their effects.

Flatten

flatfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu", "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:meanFlatRow20010621.fits"];

darkfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu", "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:meanDarkRow20010621.fits"];

starfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu",  "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:star20010621.fits"];

flatId = faReadFITS[flatfile]

darkId = faReadFITS[darkfile]

starId = faReadFITS[starfile]

Median[faData[flatId[[1]]]]

Mean[faData[flatId[[1]]]]

The original data.

ListDensityPlot[faData[starId[[1]]], 
        Mesh -> False];

faData[starId[[1]]][[1]]

Dark subtraction and flattening.  Although this should be the best, it is not.  This is mainly because the flat does not accurately represent the data.

tmp = 
        spFlatten[starId[[1]], opFlatId -> flatId[[1]], 
          opDarkId -> darkId[[1]], 
          opPrintDiagnostics -> True];

Dimensions[tmp]

ListDensityPlot[tmp, Mesh -> False];

tmp[[1]]

Dark subtraction but no flattening.

tmp2 = 
        spFlatten[
          starId[[1]],   
            opDarkId -> darkId[[1]], 
          opPrintDiagnostics -> True];

Dimensions[tmp2]

tmp2[[1]]

ListDensityPlot[tmp2, Mesh -> False];

Flattened, but no dark subtraction.  While this appears to work well on the small area seen here, on a global scale it actually makes the data worse.  This is because the flat has features the data don't.  Problems with creating a low-light flat as well.

tmp3 = 
        spFlatten[starId[[1]], 
          opFlatId -> 
            flatId[[1]],   
            opPrintDiagnostics -> True];

Dimensions[tmp3]

tmp3[[1]]

ListDensityPlot[tmp3, Mesh -> False];

No dark subtraction, no flattening.  Should be the same as input data, and it is.

tmp4 = 
        spFlatten[
          starId[[1]],   opPrintDiagnostics -> 
            True];

Dimensions[tmp4]

tmp4[[1]]

ListDensityPlot[tmp4, Mesh -> False];

tmp4[[1]] == faData[starId[[1]]][[1]]

Centroid

spMargAnal

starfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu",  "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:star20010621.fits"];

starId = faReadFITS[starfile]

res = 
      spMargAnal[starId[[1]], faData[starId[[1]]], 
        opPrintDiagnostics -> True]

{26.9591, 25.635}

1165 + 26.9591 - 1

1050 + 25.635 - 1

res[[1, 1, 1]]

res[[1, 1, 2]]

res[[1, 3]]

res = spMargAnal[starId[[1]], faData[starId[[1]]]]

res = 
      spMargAnal[starId[[1]], faData[starId[[1]]], 
        opPrintDiagnostics -> False]

spMargGaussian

Make marginal distributions.
Fit gaussian.

starfile = 
        gbMacOrUnixLocalFileName["ursula.lowell.edu",  "ursula:Local \
Files:Projects:Strip astrometry pipeline (EWD):Test \
data:star20010621.fits"];

id = faReadFITS[starfile];

data = faData[id[[1]]];

tmp = 
      spMargGaussian[id[[1]], data, 
        opPrintDiagnostics -> True]

Read in file giving reference star names, RA, Dec, approximate strip position

Somehow get planet position into here

Write out results to file that can be read by Larry's program

This will be a *.use-type file, with an extra 2 lines at the top to give date and time.