# include "relphot.h"

// 25.0 is the nominal zero point for measurements in the database
// for all measurements except PHOT_REF:
// measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0 
// for measurements using PHOT_REF:
// measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known)

int relphot_images (SkyList *skylist) {

  int i, status, Ncatalog;
  Catalog *catalog = NULL;
  FITS_DB db;

  INITTIME;

  /* register database handle with shutdown procedure */
  set_db (&db);
  gfits_db_init (&db);
  db.mode   = dvo_catalog_catmode (CATMODE);
  db.format = dvo_catalog_catformat (CATFORMAT);

  /* lock and load the image db table */
  status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT));
  if (!status && UPDATE) {
    fprintf (stderr, "error\n");
    Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
  }

  // if the file is missing, db.dbstate will have a value of either:
  // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE)
  if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) {
    Shutdown ("ERROR: database %s contains no image data", CATDIR);
  }

  // read data from Image.dat file
  if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
  MARKTIME("-- load image data: %f sec\n", dtime);

  /* load regions and images based on specified sky patch and/or catalog */
  load_images (&db, skylist, &UserPatch, FALSE, USE_ALL_IMAGES);
  MARKTIME("-- load images: %f sec\n", dtime);

  /* unlock, if we can (else, unlocked below) */
  if (!UPDATE) dvo_image_unlock (&db); 

  // load the flat correction table (if defined)
# if (0)
  char flatcorrfile[256];
  sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
# endif
  FlatCorrectionTable *flatcorr = NULL;

  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
    /* load catalog data from region files (hostID is 0 since we are not a client */
    catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL);
    MARKTIME("-- load catalog data: %f sec\n", dtime);
  
    /* add in a loop over the catalogs calling dvo_catalog_chipcoords */

    /* match measurements with images, mosaics */
    initImageBins  (catalog, Ncatalog, TRUE);
    MARKTIME("-- make image bins: %f sec\n", dtime);

    initMosaicBins (catalog, Ncatalog, TRUE);
    initGridBins   (catalog, Ncatalog);
    initMrel (catalog, Ncatalog);

    findImages (catalog, Ncatalog, TRUE);
    MARKTIME("-- set up image indexes: %f sec\n", dtime);

    findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);

    SAVEPLOT = FALSE;

    setExclusions (catalog, Ncatalog, TRUE);

    global_stats (catalog, Ncatalog, flatcorr, 0);

    if (PLOTSTUFF) {
      plot_star_coords (catalog, Ncatalog);
      // plot_mosaic_fields (catalog);
    }

    // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
    if (USE_GRID) {
      int star_toofew;

      // until we finish the grid analysis, do not reject stars out-of-hand based on ID_OBJ_FEW
      // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
      star_toofew = STAR_TOOFEW;
      STAR_TOOFEW = 0;
      for (i = 0; i < NGRID; i++) {
	STAR_BAD = ID_OBJ_POOR;
	setMrel  (catalog, Ncatalog, flatcorr);
	STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
	setMgrid (catalog, flatcorr);
      }
      STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
      STAR_TOOFEW = star_toofew;
    }

    /* determine fit values */
    for (i = 0; i < NLOOP; i++) {
      if (PLOTSTUFF) {
	plot_scatter (catalog, Ncatalog, flatcorr); 
      }
      setMrel  (catalog, Ncatalog, flatcorr); // threaded
      if (PLOTSTUFF) {
	plot_scatter (catalog, Ncatalog, flatcorr); 
      }
      setMcal  (catalog, FALSE, flatcorr);
      setMmos  (catalog, FALSE, flatcorr);
      setMgrid (catalog, flatcorr);
      MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    
      if (PLOTSTUFF) {
	plot_scatter (catalog, Ncatalog, flatcorr); 
	plot_grid (catalog, flatcorr); 
	plot_mosaics ();
	plot_images ();
	plot_stars (catalog, Ncatalog);
	plot_chisq (catalog, Ncatalog);
      }
      // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
      // if (i % 6 == 1) rationalize_images ();

      // NOTE : in the past, I was not iterating enough before cleaning.  make sure we do
      // at least 8 loops first -- that should get the systematic errors down to the ~1%
      // level, even in cases where we have an even split between photometric data and
      // data with 1 mag of extinction.

      if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); 
      if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
      if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
      if ((i > 8) && (i % 8 == 5)) clean_images ();

      global_stats (catalog, Ncatalog, flatcorr, i);
      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    }

    if (PLOTSTUFF) {
      plot_scatter (catalog, Ncatalog, flatcorr); 
      plot_grid (catalog, flatcorr); 
      plot_mosaics ();
      plot_images ();
      plot_stars (catalog, Ncatalog);
      plot_chisq (catalog, Ncatalog);
    }
  
    if (USE_GRID) dump_grid ();

    /* set Mcal & Mmos for bad images */
    int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
    setMcal  (catalog, onlyPoorImages, flatcorr);
    setMmos  (catalog, onlyPoorImages, flatcorr);
    MARKTIME("-- finalize Mcal values: %f sec\n", dtime);

    setMcalFinal (); // copy per-mosaic calibrations to the images

    if (SAVE_IMAGE_UPDATES) {

      FITS_DB dbX;
      off_t *LineNumber;
      Image *image, *subset;
      off_t Nimage, Nsubset, i, Nx;
      char filename[1024];

      gfits_db_init (&dbX);
      dbX.lockstate = LCK_XCLD;
      dbX.timeout   = 60.0;
      dbX.mode      = db.mode;
      dbX.format    = db.format;

      snprintf (filename, 1024, "%s.bck", ImageCat);
      if (!gfits_db_lock (&dbX, filename)) {
	fprintf (stderr, "can't lock backup image image catalog\n");
	return (FALSE);
      }
    
      // apply the changes from the image subset to the full table:

      // convert database table to internal structure (binary to Image)
      // 'image' points to the same memory as db->ftable->buffer
      image = gfits_table_get_Image (&db.ftable, &Nimage, &db.scaledValue, &db.nativeOrder);
      if (!image) {
	fprintf (stderr, "ERROR: failed to read images\n");
	exit (2);
      }
      gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);

      subset = getimages (&Nsubset, &LineNumber);
      for (i = 0; i < Nsubset; i++) {
	if (LineNumber[i] == -1) continue;
	memcpy (&image[LineNumber[i]], &subset[i], Nx);
      }

      // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
      // the entire Image table, but with the updates in place I think this says: re-work the
      // ftable/vtable usage in this program to be more sensible...
      gfits_copy_header (&db.header,  &dbX.header);
      gfits_copy_matrix (&db.matrix,  &dbX.matrix);
      gfits_copy_header (&db.theader, &dbX.theader);
      gfits_copy_ftable (&db.ftable,  &dbX.ftable);

      dbX.ftable.header = &dbX.theader;
      dbX.virtual = FALSE;

      // save Images.dat using the copied structure
      if (UPDATE_CATFORMAT) {
	// ensure the db format is updated
	dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
	char photcodeFile[1024];
	sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
	SavePhotcodesFITS (photcodeFile);
      }
      dvo_image_save (&dbX, VERBOSE);
      dvo_image_unlock (&dbX); 
      gfits_db_free (&dbX);
      MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
    }

    /* at this point, we have correct cal coeffs in the image/mosaic structures */
    for (i = 0; i < Ncatalog; i++) {
      // these tiny values are set by BrightCatalogSplit from load_catalogs
      free_tiny_values (&catalog[i]);
      dvo_catalog_free (&catalog[i]);
    }
    free (catalog);
    freeImageBins (Ncatalog, TRUE);
    freeMosaicBins (Ncatalog, TRUE);
    freeGridBins (Ncatalog);

    // end of if (NLOOP > 0) block : this loop determines the offsets per chip 
  } else { 
    // If nloop == 0, the above pass is not performed, in which case
    // the Mcal values passed to the mosaics are not returned to the images...
    setMcalFinal (); // copy per-mosaic calibrations to the images
  }

  reload_images (&db);

  // only change the real database files if -update is requested
  if (!UPDATE) exit (0);
  
  /* Load catalog data from region files, update Mrel include all data.  In a parallel
     context, this function writes the image parameters as a subset table for the remote
     clients */
  if (!CALIBRATE_STACKS_AND_WARPS) {
    reload_catalogs (skylist, flatcorr, 0, NULL);
    MARKTIME("-- updated all catalogs: %f sec\n", dtime);
  }

  if (UPDATE_CATFORMAT) {
    // ensure the db format is updated
    db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
  }
  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
    // do not save changes if we did not make changes. 
    dvo_image_update (&db, VERBOSE);
  }
  dvo_image_unlock (&db); 

  freeImages (db.ftable.buffer);
  gfits_db_free (&db);

  return TRUE;
}

