# include <pslib.h>

// ignoring things like allocation, error tests, etc

main (int argc, char **argv) {

  psArray headers;
  psMetadata *base, *site, *camera, *recipe;

  // the set of code below (up to END CONFIG) is generic boilerplate which can be used for most analysis stages
  // the specific issues are defined by the recipefile loaded for this analysis
  // I've written this keeping the external metadata separated (base, site, camera, recipe),
  // but we could put all of this in a single metdata collection, called eg config

  // USAGE: phase2 (input) (outroot) [-camera (camera)] [-recipe (recipe)]
  input = argv[1];
  outroot = argv[2];

  fits = psFitsAlloc (input);

  // load only the headers from the given file
  headers = psFitsReadHeaderSet (fits);
 
  // identify camera (may be specified on command line)
  if (cameraname == NULL) {
    cameraname = pmCameraFromHeader (psListGet (headers->list, 0));
  }
  camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname);

  // load the camera details for the specific camera
  camera = psMetadataParseConfig (NULL, camerafile, FALSE);

  // recipefile is defined in the camera metadata, if not specified on command line
  if (recipefile == NULL) {
    recipefile = psMetadataLookup (camera, "RECIPE.PHASE2");
    recipefile = DetrendDatabaseLookup ("recipe", camera, header);   // an alternative
    // make this a database lookup & a function of date
  }
  recipe = psMetadataParseConfig (NULL, recipefile, FALSE);  

  // END CONFIG loading code

  // check the validity of the given set of headers based on expectactions for this camera
  status = pmCameraValidateHeaders (headers, camera);

  // construct a complete, empty (ie, no pixel data) FPA structure for this camera / header set
  // we need to construct the output container.  is it big enough to fit in memory?  assume so?
  FPA = pmFPAfromHeader (headers, camera);
  
  // I'm pretending the extname value for each cell in part of psCell
  // Do we need to construct a list of the extname values associated with each readout?
  // Is EXTNAME invarient?  (probably not; does not always exist)

  /* This will work, but it is not efficient: 
     we should work with the extensions in sequence */
  for (i = 0; i < FPA[0].chips.n; i++) {
    chip = FPA[0].chips.data[i];

    if (Base == CELL) {
      f = psFitsMoveExtName(f, extname);
    }
    
    for (j = 0; j < chip[0].cells.n; j++) {
      cell = chip[0].cells.data[j];

      if (Base == CELL) {
	f = psFitsMoveExtName(f, extname);
      }

      /* these need to work with the psFits pointer correctly 
      bias = load_bias_image (recipe, cell[0].metadata, cell[0].extname);
      dark = load_dark_image (recipe, cell[0].metadata, cell[0].extname);
      flat = load_flat_image (recipe, cell[0].metadata, cell[0].extname);
      badpix = load_badpix_image (recipe, cell[0].metadata, cell[0].extname);
      fringe = load_fringe_image (recipe, cell[0].metadata, cell[0].extname);
      */

      for (k = 0; k < cell[0].readouts.n; k++) {
	readout = cell[0].readouts[k];

	// load the specific readout pixels.  put this in a function?  
	status = pmReadoutLoad (fits, cell, k);

	kernel = make_kernel ();
	mask = make_mask ();

	bias_subtract_readout (readout, bias, recipe);
	dark_subtract_readout (readout, dark, recipe);
	nonlinear_correct_readout (readout, mask, recipe);

	trim_readout (readout, mask, recipe);
	mask_readout (readout, mask, badpix, kernel, recipe);

	// where do we convolve the flat with the kernel?
	pmFlatField (readout, mask, flat);

	defringe_readout (readout, mask, fringe, recipe);
      }
    }
  }	
}

// bias subtract the given readout according to the recipe (bias image may be NULL)
// no mask needed yet
bias_subtract_readout (psReadout *readout, psImage *bias, psMetadata *recipe) {
  
  char *biassec;
  psImage *overscan;
  psList *overscanlist;
  pmOverscanAxis direction;
  psRectangle biassec;

  // determine overscan option
  biassec = pmRecipeDefineSection (recipe, readout[0].metadata, "BIAS.OVERSCAN"); 
  overscan = psImageSubsection (readout[0].image, biassec);

  psListAdd (overscanlist, overscan, PS_LIST_HEAD);             // list of overscans (allocate list first)

  // the direction is implied by noting the long dimension of the biassec
  direction = DirectionFromOverscan (biassec);  // weak, but do we want a function to deal with this or add it inline?

  // the type of statistic is defined by the recipe: BIAS.OVERSCAN.STATS
  statname = psMetadataLookup (recipe, "BIAS.OVERSCAN.STATS");  // string defining overscan stats
  stats = DefineStats (statname);

  // identify the type of fit requested and relevant parameters
  // the type of fit is defined by the recipe: BIAS.OVERSCAN.STATS
  fitname = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT");  // string defining overscan stats

  if (!strcmp (fitname, "NONE")) {
    fitType = PM_FIT_NONE;
    fitSpec = NULL;
    nBin = 0;
  }

  if (!strcmp (fitname, "SPLINE")) {
    fitType = PM_FIT_SPLINE;
    nBin = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NBIN");
    nPts = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NPTS");
    order = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.ORDER");  // only 1 or 3 is allowed for 'order'
    fitSpec = psSpline1DAlloc (nPts, 3, 0, Npix);                  // Npix is either NAXIS1 or NAXIS2 (Nx or Ny);
 }

  if (!strcmp (fitname, "POLYNOMIAL")) {
    fitType = PM_FIT_POLYNOMIAL;
    nBin = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NBIN");
    nPts = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.NPTS");
    order = psMetadataLookup (recipe, "BIAS.OVERSCAN.FIT.ORDER");
    fitSpec = psPolynomial1DAlloc (nPts, PS_POLYNOMIAL_CHEB);
  }

  readout[0].image = pmSubtractBias (readout[0].image, fitSpec, overscanlist, direction, stats, nBin, PM_FIT_NONE, bias);

  // add these to the readout[0].metadata? 
  // or to the chip / cell metadata?

}

// subtract dark image using the appropriate scaling factor for this exposure time
// no mask needed yet
dark_subtract_readout (psReadout *readout, psImage *dark, psMetadata *recipe) {
  
  char *biassec;
  psImage *overscan;
  psList *overscanlist;
  pmOverscanAxis direction;

  if (dark == NULL) {
    return;
  }
    
  // need to get a lookup table from the database, apply exptime for readout and dark
  factor_source = psMetadataLookup (recipe, "DARK.FACTOR"); // use the recipe to determine the dark scaling
  factor = getScalingFactor (readout, dark, factor_source); 

  // out = in - dark * factor
  tmpdark = psBinaryOp (NULL, dark, "*", factor);
  readout[0].image = psBinaryOp (readout[0].image, readout[0].image, "-", tmpdark);

  psFree (tmpdark);
}

// this corrects the non-linear response, but probably should use data determined at the cell level
// we need to set a mask value for pixels which have were out of range (TBD)
nonlinearity_correct_readout (psReadout *readout, psReadout *mask, psMetadata *recipe) {
  
  char *method;
  psPolynomial1D correction;
  psVector inFlux, outFlux;

  method = psMetadataLookup (recipe, "NONLINEAR.METHOD");

  // use a polynomial; get coeffs from recipe
  if (!strcmp (method, "POLYNOMIAL")) {
    correction.n = psMetadataLookup (recipe, "NONLINEAR.ORDER"); 
    correction.coeffs = psMetadataLookup (recipe, "NONLINEAR.COEFFS"); 
    pmNonLinearityPolynomial (readout, correction);
  }

  if (!strcmp (method, "LOOKUP")) {
    lookup_table = DetrendDatabaseLookup (time, cell, camera, etc?);
    inFlux = psFITSTableRead (lookup_table, "IN_FLUX"); 
    outFlux = psFITSTableRead (lookup_table, "OUT_FLUX"); 
    pmNonLinearityLookup (readout, inFlux, outFlux);
  }
}

// trim the image and mask based on the requested region
trim_readout (psReadout *readout, psReadout *mask, psMetadata *recipe) {
  
  datasec = pmRecipeDefineSection (recipe, readout[0].metadata, "TRIM.DATASEC"); 

  psImageTrim (readout[0].image, datasec);
  psImageTrim (mask[0].image, datasec);

}

// apply the bad pixel mask to the image mask
mask_readout (psReadout *readout, psImage *mask, psImage *badpix, psKernel *kernel, psMetadata *recipe) {
  
}

psImage *load_bias_image (psMetadata *recipe, psMetadata *header, char *extname) {

  // load in the complete corresponding bias image (this should be done once per cell)
  // we are implying that the bias image is guaranteed to have the corresponding EXTNAME
  bias_source = psMetadataLookup (recipe, "BIAS.IMAGE");
  biasname = NULL;
  if (!strncmp (bias_source, "FILE:", 5)) {
    biasname = &bias_source[5];
  }
  if (!strcmp (bias_source, "DB:BEST")) {
    biasname = DetrendDatabaseLookup ("best", header);   // these APIs need to be seriously fleshed out!!
  }
  if (!strcmp (bias_source, "DB:CLOSE")) {
    biasname = DetrendDatabaseLookup ("close", header);   // these APIs need to be seriously fleshed out!!
  }
  if (biasname != NULL) {
    bias     = psImageReadSection (NULL, 0, 0, 0, 0, 0, extname, 0, biasname);
  }

  return (bias);
}

psImage *load_dark_image (psMetadata *recipe, psMetadata *header, char *extname) {

  // load in the complete corresponding bias image (this should be done once per cell)
  // we are implying that the bias image is guaranteed to have the corresponding EXTNAME
  dark_source = psMetadataLookup (recipe, "DARK.IMAGE");
  darkname = NULL;
  if (!strncmp (dark_source, "FILE:", 5)) {
    darkname = &dark_source[5];
  }
  if (!strcmp (dark_source, "DB:BEST")) {
    darkname = DetrendDatabaseLookup ("best", header);   // these APIs need to be seriously fleshed out!!
  }
  if (!strcmp (dark_source, "DB:CLOSE")) {
    darkname = DetrendDatabaseLookup ("close", header);   // these APIs need to be seriously fleshed out!!
  }
  if (darkname != NULL) {
    dark     = psImageReadSection (NULL, 0, 0, 0, 0, 0, extname, 0, darkname);
  }

  return (dark);
}

psRectangle *pmRecipeDefineSection (psMetadata *recipe, psMetadata *header, char *name) {

  char *source, *section;
  psRectangle *Section;

  source = psMetadataLookup (recipe, name);

  if (!strcasecmp (source, "NONE")) {
    Section = psRectangleAlloc (0, 0, 0, 0); // the NULL section
    return (Section);
  }

  if (!strncasecmp (source, "HEADER:", 7)) {
    char *keyword;
    
    keyword = &source[7]; 
    section = psMetadataLookup (readout[0].metadata, keyword);
    Section = pmSectiontoRectangle (section);
    psFree (section);
    psFree (source);
    return (Section);
  }

  if (!strncasecmp (source, "RECIPE:", 7)) {
    section = &source[7]; 
    Section = pmSectiontoRectangle (section);
    psFree (source);
    return (Section);
  }

  psError ("invalid section");
  return (NULL);

}


bool pmConfigLoadBasic (psMetadata **site, psMetadata **base) {

  // base file needs to be hardwired, from environment, or based on program name
  base   = psMetadataParseConfig (NULL, basefile, FALSE);  

  // load site-specific information (such as database servers, etc)
  sitefile = psMetadataLookup (base, "SITE");
  site   = psMetadataParseConfig (NULL, sitefile, FALSE);  

}

bool pmConfigLoadCamera (psMetadata **camera, psMetadata *site, psMetadata *header) {
  
  if (site == NULL) return (FALSE);

  // identify camera (may be specified on command line)
  cameraname = getargs ("-camera");  // pseudo-code, not a valid API

  // if there is no source for camera info, fail
  if ((cameraname == NULL) && (header == NULL)) {
    return (FALSE);
  }

  // try to determine camera from header 
  if (cameraname == NULL) {
    cameraname = pmCameraFromHeader (header);
  }

  // find appropriate camera definition file
  camerafile = psMetadataLookup (site, "CAMERA.%s", cameraname);

  // load the camera details for the specific camera
  *camera = psMetadataParseConfig (NULL, camerafile, FALSE);

  return (TRUE);
}

bool pmConfigLoadRecipe (psMetadata **recipe, psMetadata *camera, char *script) {

  char recipename[64];

  if (site   == NULL) return (FALSE);
  if (camera == NULL) return (FALSE);

  // define command-line options
  recipefile = getargs ("-recipe");  // pseudo-code, not a valid API

  // recipefile is defined in the camera metadata, if not specified on command line
  if (recipefile == NULL) {
    sprintf (recipename, "RECIPE.%s", script);

    // option 1: recipe is in the camera definition 
    recipefile = psMetadataLookup (camera, recipename);

    // option 2: recipe is in the metadata database (equiv to a detrend file)
    // this option requires a different API (header data, etc)
    recipefile = DetrendDatabaseLookup (recipename, site, camera, header);
  }

  *recipe = psMetadataParseConfig (NULL, recipefile, FALSE);  
  return (TRUE);
}
