/***************************************************************************
 *
 * findref - Search for an archive for reference star candidates.
 *
 * Usage:
 *    findref Archive OutFile [-r radius] [-s sigma] [-v]
 *
 * where:
 *    Archive = name of the photometry archive to search
 *    OutFile = name of candidate file to create (ascii text)
 *    -r radius = set search radius in arcsec [default: 30 arcsec]
 *    -s sigma = set rejection threshold [default 0.015 mag]
 *    -v = verbose output
 *
 * Searches an archive for all Type 11 stars within a given radius
 * of the microlensing event that satisfies a maximum rms variation
 * criterion.  The default search radius is 30", and the default
 * maximum rms variation is 0.015 magnitudes.  The search radius and
 * rejection threshold may be changed by using the -r and -s flags
 * respectively.  
 *
 * The pixel scale is derived from pixscale=timeinfo.Seeing/timeinfo.FWHM,
 * we store every other fool thing in the archives but...  Sigh.  If
 * pixscale turns up bogus, it defaults to 1.0.
 *
 * Produces a candidate list in an ascii format data file.
 *
 * R. Pogge, OSU Astronomy
 * pogge@astronomy.ohio-state.edu
 * 2000 April 30
 *
 * Modification History:
 *    2000 May 1: changed logic to search in arcsec, using the pixel
 *                scale teased out of timeinfo.  Also adds the mean
 *                magnitude of the lens to the report. [rwp/osu]
 *    2000 Aug 2: Fixed the TypMean problem, whereby JP thought he
 *                would save space by coding the observation number and
 *                the mean DoPhot type as one 2-byte integer.  [rwp/osu]
 *
 ***************************************************************************/

#include "fsui.h"

#undef __DEBUG

void Usage();

main(int narg, char *argv[]) 
{ 
  SUIVIFIP *sfp;   /* archive file pointer      */
  char ArchRoot[80];
  char ArchFile[256];

  FILE *fpRef;     /* output ascii file pointer */
  char RefRoot[80];     
  char RefFile[256]; 

  /* arguments from the command line */

  int a;
  float SearchRad = 30.0;
  float SearchPix = 30.0;
  float MaxRMS = 0.015;
  int Verbose = 0;
  int iarg = 1;

  /* archive data structures */

  STARINFO *star_i;
  TIMEINFO timeinfo;
  MESURE mesure;

  int j;
  
  int NStars = 0;
  int NObs = 0;
  int NFound = 0;

  int LensID;
  int NObsLens;
  int StarID; 
  int NObsStar;
  float Mag;
  float MagErr;
  float MagRMS;
  int MeanTyp;

  float X, Y;
  float Xc, Yc;
  float R, dx, dy;
  float pixscale = 1.0;

  /* Some initializations */

  star_i = &StarInfo; 

  /* the command line must contain at least 3 arguments (cmd+2) */

  if (narg < 3) {
    Usage();
    exit(1);
  }

  /* default search directories; cwd if undefined */

  if ((InDir=getenv("IN_FILE")) == NULL) InDir  = "";
  if ((OutDir=getenv("OUT_FILE")) == NULL) OutDir = "";

  /* input archive and output candidate list names */
  
  strcpy(ArchRoot,argv[iarg++]);          
  sprintf(ArchFile ,"%s%s",InDir,ArchRoot);

  strcpy(RefRoot,argv[iarg++]);          
  sprintf(RefFile ,"%s%s",OutDir,RefRoot);

#ifdef __DEBUG
  fprintf(stdout,"Archive File: %s\n",ArchFile);
  fprintf(stdout,"Candidate List File: %s\n",RefFile);
#endif

  /* open the archive file readonly */

  if ((sfp = SuiviOpen(ArchFile,SUOF_RO)) == NULL) {
    fprintf(stderr,"Error opening photometry archive %s\n",ArchFile); 
    exit(1); 
  }

  /* test to make sure that the data is not byte swapped */

  NStars = SuiviGetNbStars(sfp); 
  if ((NStars < 0) || (NStars > 100000)) { 
    fprintf(stdout,"Archive data are byte swapped.\n");
    fprintf(stdout,"please convert %s to unswapped 98 format.\n",ArchRoot);
    exit(1);
  }

  /* Parse the rest of the command-line arguments.  These
   * are options to change the search criteria, etc. */

  for (a=iarg; a<narg; a++) {
    switch(tolower(argv[a][1])) {
    case 'v': 
      Verbose = 1;
      break;  
    case 'r':
      SearchRad = atof(argv[++a]);
      break;
    case 's':
      MaxRMS = atof(argv[++a]);
      break;
    default: 
      fprintf(stderr,"Syntax error: unknown option %s\n", argv[a]); 
      Usage();
      exit(1);
      break;
    }
  }

  /* Open the output ascii photometry file */

  if ((fpRef =  fopen(RefFile,"a")) == NULL) {
    fprintf(stderr,"Error opening candidate list file %s\n",RefFile);
    exit(1);   
  }  

  /* 
   * Read in the Archive comment header fields.  A good acid test of the
   * validity of the archive
   */

  if(SuiviReadComment(sfp,CommStrg) != 0) { 
    fprintf(stderr,"Error: cannot read comment fields from archive\n");
    exit(1); 
  }

  /* Get the number of observations and stars/observation in the archive */

  NObs = SuiviGetNbMesures(sfp);      
  NStars = SuiviGetNbStars(sfp);

#ifdef __DEBUG
  fprintf(stdout,"Archive %s contains %d observations of %d stars\n",
	  ArchRoot,NObs,NStars);
#endif

  /* Get the ID of the microlensing star */

  GetMicrostarNumPL(sfp,&LensID);

  /* Extract the data for the microlensing event */
  
  if (SuiviReadStarInfo(sfp,LensID,(char *)(star_i)) != 0) {
    fprintf(stderr,"Error reading microlensing info from %s\n",ArchRoot); 
    exit(1);
  }

  Mag      = (float) star_i->FlxMean;
  MagErr   = (float) star_i->SigRef/1000.0;
  MagRMS   = (float) star_i->SigMean/1000.;
  NObsLens = (int) (star_i->TypMean/100.0) ;
  MeanTyp  = star_i->TypMean-(NObsLens*100.0);

  /* Sigh.  JP tried to "saved space" by coding the observation number
   * and mean DoPhot type together in a single 2-byte integer.  This
   * breaks if the number of observations exceeds 327, resulting in
   * negative values for both (the 2-byte signed integer coding offset
   * is -65536, the first 3 for the Obs# and the last 2 for the DoPhot
   * type).  This fixes it. */

  if (NObsLens < 0)
    NObsLens = NObsLens + 655;
  if (MeanTyp < 0) 
    MeanTyp = MeanTyp + 36;

  /* Get the position of the microlensing event */

  Xc = star_i->X;
  Yc = star_i->Y;

  /* extract the timeinfo records for the first record.  Pray it is right */

  if (SuiviReadTimeInfo(sfp,1,1,(char *)&timeinfo) != 0) {
    fprintf(stderr,"Error: cannot read TimeInfo data from %s\n",ArchRoot);
    exit(1);
  }

  if (timeinfo.Seeing == 0.0 || timeinfo.FWHM == 0.00) {
    pixscale = 1.0;
  } else {
    pixscale = timeinfo.Seeing / timeinfo.FWHM;
  }

  SearchPix = SearchRad / pixscale;

  if (Verbose) {
    fprintf(stdout,"Lens: %d %f %f %f %f %f %d %d \n",LensID, Xc,Yc,
	    Mag,MagErr,MagRMS,NObsLens,MeanTyp); 
    PrintStarInfo(star_i);
  }

  /* write the candidate list file header */

  fprintf(fpRef,"#\n");
  fprintf(fpRef,"# Archive: %s\n",ArchRoot);
  fprintf(fpRef,"# Lens: Star %d @ X=%.2f Y=%.2f Mean=%5.2f\n",
	  LensID,Xc,Yc,Mag);
  fprintf(fpRef,"# Search Radius: %.1f arcsec (%d pixels)\n",
	  SearchRad, (int)(SearchPix));
  fprintf(fpRef,"# Rejection Threshold: %5.3f\n",MaxRMS);
  fprintf(fpRef,"# DoPhot Type: 11\n");
  fprintf(fpRef,"#\n");

  /* create the header for the rest of the candidate list */

  fprintf(fpRef,"# Star   X       Y     Mag    Err   RMS  Typ\n");
  
  /* Search for suitable reference stars archive */
  
  for (j=0; j<NStars; j++) {

    NObs = SuiviGetNbMesures(sfp);
    StarID = j+1;

    if (StarID != LensID) {
      if (SuiviReadStarInfo(sfp,StarID,(char *)(star_i)) != 0) {
	fprintf(stderr,"Error: cannot read data for star %d in %s\n",
		StarID,ArchRoot);
	exit(1);
      }

      /* get the X,Y coordinates of the star */

      X = star_i->X;
      Y = star_i->Y;

      /* compute the radial distance from the lensing event */

      dx = X - Xc;
      dy = Y - Yc;
      R = (float)(sqrt((double)(dx*dx + dy*dy)));

      /* extract the mean photometric data */

      Mag      = (float) star_i->FlxMean;
      MagErr   = (float) star_i->SigRef/1000.0;
      MagRMS   = (float) star_i->SigMean/1000.;
      NObsStar = (int) (star_i->TypMean/100.0) ;
      MeanTyp  = star_i->TypMean-(NObsStar*100.0);

      /* once again, fix if NObsStar>327... */

      if (NObsStar < 0)
	NObsStar = NObsStar + 655;
      if (MeanTyp < 0) 
	MeanTyp = MeanTyp + 36;

      /* 
       * A star is acceptable as a candidate reference star if:
       *   1) the star is within the SearchPix pixels of the lens.
       *   2) its rms brightness variations are less than MaxRMS
       *   3) it is a DoPhot type 11 star
       * The star must satisfy all three criteria.
       */ 

      if ((R <= SearchPix) && (MagRMS <= MaxRMS) && (MeanTyp == 11)) {
	NFound++ ;
	fprintf(fpRef,"%5d %7.2f %7.2f %6.3f %5.3f %5.3f %3d\n",
		StarID, X, Y, Mag, MagErr, MagRMS, MeanTyp);
      }

    }

  }

  /* report the results */

  if (NFound == 0) {
    fprintf(stdout,"No candidate reference stars found in the archive!\n");
    fprintf(stdout,"  Check your search criteria & try again\n");
  } else {
    fprintf(stdout,"%d reference stars found in archive %s\n",NFound,ArchRoot);
  }

  /* Close the files */

  SuiviClose(sfp);
  fclose(fpRef);

  /* all done */

}

/***************************************************************************
 *
 * Usage() - print the usage message on errors
 *
 ***************************************************************************/

void Usage()
{
  fprintf(stdout,"\nfindref - search an archive for reference star candidates.\n\n");
  fprintf(stdout,"Usage: findref Archive RefFile [-r rad] [-s sig] [-v]\n\n");
  fprintf(stdout,"   Archive: name of the PLANET archive file to read\n");
  fprintf(stdout,"   RefFile: name of file of reference star candidates to create\n");
  fprintf(stdout,"   -r rad: search radius in arcsec [default: 30 arcsec]\n");
  fprintf(stdout,"   -s sigma: rms rejection threshold [default: 0.015 mag]\n");
  fprintf(stdout,"   -v: verbose output to tty\n\n");
  return;
}