(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialization Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)









































































BeginPackage["jleGroup`options`"];



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



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.";



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.";



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.";



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.";



EndPackage[];





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



SetOptions[lsIterateToConvergence,opCheckImaginary\[Rule]True]





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[spAnalyzeStrip]=
    {opPrintDiagnostics\[Rule]False,
      opPlanetName\[Rule]"Planet",
      opPlanetMagnitude\[Rule]0.0,
      opPlanetNear\[Rule]250,
      opDoAstrometryCentroiding\[Rule]True};

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

Options[spExtractStar]=
    {opOutputFile\[Rule]Null, 
      opPrintDiagnostics\[Rule]False};

Options[spFlatten]=
    {opFlatId\[Rule]Null,
      opDarkId\[Rule]Null,
      opOverscanPixels\[Rule]{},
      opPrintDiagnostics\[Rule]False
      };

Options[spCentroid]=
    {opCentroidFunction\[Rule]"margAnal",
      opPrintDiagnostics\[Rule]False,
      opReturn\[Rule]"CP"};



lsModelContext[spGaussian] = "`spGaussian`";



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





Begin["`Private`"];

















spMakeMeanRow[infile_, frow_, opts___]:=
  Module[
    {outfile,saveToFile,diag,nrowstep,nrows,
      hid,na,na1,na2,na3,
      meanRow,
      i,fList,rList,did,nread,lastrow},
    
    (*Get options*)
    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];
    
    (*Verify that infile frame exists*)
    (*Note:  currently in fitsAide, 
      there is no way to verify that an existing file is a SIMPLE=
        TRUE FITS file, 
      so if the file exists we will assume it is of the correct type.  
          If it's not, this will take a very long time.*)
    
    If[FileType[infile]=!=File,
      Print["Input frame ",infile," does not exist.  Exiting."];
      Return[-1]];
    
    (*Get frame size parameters*)
    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]];
    (*Note:  we will assume 2D data.  However, 
      since the test data are 3D with NAXIS3=1, we will allow that as well.*)

        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]]];
    
    (*Check validity of firstRow variable*)
    If[frow>na2,
      Print["First row (",frow,") is > number of rows in image (",na2,
        ").  Exiting."];
      Return[-1]];
    
    (*Initialize mean row*)
    meanRow=Table[0,{na1}];
    
    (*Loop over rows*)
    
    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}];
    
    (*Get mean row*)
    meanRow/=(na2-frow+1);
    
    (*Store to file if requested*)
    If[saveToFile,
      newHeader=spMakeMeanHeader[infile,hid,na1];
      newId=faNewDataset[newHeader,meanRow];
      If[diag,
        Print["Writing to FITS file ",outfile]];
      faWriteFITS[outfile,newId]];
    
    (*Return mean row*)
    meanRow
    ]



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","",""}}
    ]





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











spExtractStar[infile_, frow_, fcol_,nrow_,ncol_,opts___]:=
  Module[
    {outfile,saveToFile,diag,hid,na,na1,na2,na3,sfid,newdata,newHeader,
      newId},
    
    (*Get options*)
    outfile=opOutputFile/.{opts}/.Options[spExtractStar];
    saveToFile=If[outfile===Null,False,True];
    diag=opPrintDiagnostics/.{opts}/.Options[spExtractStar];
    
    (*Verify that infile frame exists*)
    (*Note:  currently in fitsAide, 
      there is no way to verify that an existing file is a SIMPLE=
        TRUE FITS file, 
      so if the file exists we will assume it is of the correct type.  
          If it's not, this will take a very long time.*)
    
    If[FileType[infile]=!=File,
      Print["Input frame ",infile," does not exist.  Exiting."];
      Return[-1]];
    
    (*Get frame size parameters*)
    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]];
    (*Note:  we will assume 2D data.  However, 
      since the test data are 3D with NAXIS3=1, we will allow that as well.*)

        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]]];
    
    (*Check validity of requested extraction*)
    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]];
    
    (*Extract subimage*)
    If[na===3,
      sfid=faPartialReadFITS[infile,1,{fcol,frow,1},{ncol,nrow,1},{1}];
      newdata=faData[sfid][[1]],
      (*else*)
      
      sfid=faPartialReadFITS[infile,1,{fcol,frow},{ncol,nrow},{1}];
      newdata=faData[sfid]];
    
    (*Store to file if requested*)
    
    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]];
    
    (*Return extracted subimage header and data id*)
    newId
    ]



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","",""}}
    ]





















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[row_,overscanPixels_,opts___]:=
  Module[
    {bias},
    
    bias=Mean[row[[overscanPixels]]];
    
    row-bias
    ]











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[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*)
        }]
    ]





spaussian[p_,c_]:=
  Module[
    {x0,sigma,area},
    
    {x0,sigma,area,bkg}=p;
    
    N[(area/(sigma Sqrt[2Pi]))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\[Rule]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\[Rule]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[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\[Rule]"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[2Pi]);
    
    (*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\[Rule]disp];*)
 
       iteration=1;
    maxFracParamError=10^-3;
    maxIterations=10;
    For[
      Null,
      iteration==1||
        maxFracError>maxFracParamError&&iteration\[LessEqual]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\[LeftDoubleBracket]lsFitList\[RightDoubleBracket]]
        ]
      ];
    
    (*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 := \[IndentingNewLine]Module[\[IndentingNewLine]{index, 
        fracErrorList = {}}, \[IndentingNewLine]\[IndentingNewLine]Do[
        index = lsFitList\[LeftDoubleBracket]
            j\[RightDoubleBracket]; \(AppendTo[fracErrorList, 
            Abs[\(lsParamHist\[LeftDoubleBracket]\(-1\), index\
\[RightDoubleBracket] - lsParamHist\[LeftDoubleBracket]\(-2\), index\
\[RightDoubleBracket]\)\/jleGroup`leastSquares`lsError\[LeftDoubleBracket]j\
\[RightDoubleBracket]]];\), {j, Length[lsFitList]}]; 
      Max[fracErrorList]]\)\(\[IndentingNewLine]\)\)\)



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\[Rule]{Round[ctr]-1,{2,2}}*)]];
    
    delta=astromTable[[refStarNumber,{2,4}]]-ctr[[{1,2}]];
    If[diag,
      Print["\tOffset (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["\tCentroiding #",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\[Rule]"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<>"\"","-"\[Rule]"/"];
      timeline="DATATIME= "<>ut;
      astromStrings=spOutputStrings[astromTable,astromResults];
      outtable=Flatten[{dateline<>"\n",timeline<>"\n",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<>"\n",timeline<>"\n",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<>"\"","-"\[Rule]"/"];
      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["\tCentroiding #",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\[Rule]"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["\tplanetRowColumn=",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["\tplanetRowColumn=",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\[Rule]"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["\tFitted 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<>"\n",timeline<>"\n",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<>"\n",timeline<>"\n",targetStrings2}];
    spWriteFititFile[outFileBase<>".stars",outtable];
    
    (*Return*)
    {astromResults,targetResults,planetResults}
    ]



spGetPlanetPosition[jplEphFile_,ymd_,ut_]:=
  Module[
    {obsmjd,pid,pdata,mjd,radecrange,ra,dec,rai,deci},
    
    (*Start time of strip*)
    
    obsmjd=tcDateTimeStringToMJD[
        StringReplace[ymd<>ut,{"-"\[Rule]" ",":"\[Rule]" "}]];
    
    (*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[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]]<>"\n")&/@stable
    ]

toString[expr_,{n_,f_}]:=
  Module[
    {s,len},
    s=ToString[
        NumberForm[N[expr],{n,f},NumberPadding\[Rule]{"","0"},
          ExponentFunction\[Rule](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\[Rule]{"","0"},
          ExponentFunction\[Rule](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[fileName_]:=
  Module[
    {fp,table},
    
    fp=OpenRead[fileName];
    Skip[fp,String];
    Skip[fp,String];
    table=Import[fp,"Table"];
    Close[fp];
    
    table
    ]



spWriteFititFile[outputFileName_,outputTable_,opts___]:=
  Module[
    {fp},
    
    fp=OpenWrite[outputFileName];
    WriteString[fp,#]&/@outputTable;
    Close[fp];
    ]





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[(* "`Private`" *)];

Protect["asbGroup`stripPipeline`*"];

EndPackage[];





























































































































































































































































































































































































































