# include "dvoshell.h"

void initPhotcodeSequence (int Nsecfilt);
void freePhotcodeSequence ();
void printPhotcodeSequence (Average *average, SecFilt *secfilt, int entry, int type);

int gstar (int argc, char **argv) {
  
  char *date;
  double Ra, Dec, Radius, Radius2, r, dec0, dec1;
  double Mcat, Mrel;
  double *RA, *DEC;
  int i, j, k, m, N, *N1, Nsecfilt, NPTS, QUIET, FULL_OUTPUT;
  int Nstars, found, GetMeasures, Nlo, Nhi;
  int SaveVectors;
  Vector *vec1, *vec2, *vec3, *vec4;
  SkyTable *sky;
  SkyList *skylist;
  Catalog catalog;
  int TimeFormat;
  time_t TimeReference;

  if (!InitPhotcodes ()) return (FALSE);
  Nsecfilt = GetPhotcodeNsecfilt ();

  QUIET = FALSE;
  if ((N = get_argument (argc, argv, "-q"))) {
    QUIET = TRUE;
    remove_argument (N, &argc, argv);
  }

  NPTS = 0;
  vec1 = vec2 = vec3 = vec4 = NULL;
  SaveVectors = FALSE;
  if ((N = get_argument (argc, argv, "-save"))) {
    remove_argument (N, &argc, argv);
    SaveVectors = TRUE;
    if ((vec1 = SelectVector ("gs:m", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    if ((vec2 = SelectVector ("gs:t", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    if ((vec3 = SelectVector ("gs:z", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    if ((vec4 = SelectVector ("gs:f", ANYVECTOR, TRUE)) == NULL) return (FALSE);
  }

  FULL_OUTPUT = FALSE;
  if ((N = get_argument (argc, argv, "-full"))) {
    FULL_OUTPUT = TRUE;
    remove_argument (N, &argc, argv);
  }

  GetMeasures = FALSE;
  if ((N = get_argument (argc, argv, "-m"))) {
    GetMeasures = TRUE;
    remove_argument (N, &argc, argv);
  }
  if ((N = get_argument (argc, argv, "-meas"))) {
    GetMeasures = TRUE;
    remove_argument (N, &argc, argv);
  }

  if (argc != 4) {
    gprint (GP_ERR, "USAGE: gstar RA DEC Radius [-m]\n");
    return (FALSE);
  }
  
  GetTimeFormat (&TimeReference, &TimeFormat);

  Ra = atof (argv[1]);
  Dec = atof (argv[2]);
  Radius = atof (argv[3]);

  while (Ra < 0.0) Ra += 360.0;
  while (Ra > 360.0) Ra -= 360.0;
  
  /* load sky from correct table */
  sky = GetSkyTable ();
  skylist = SkyListByRadius (sky, -1, Ra, Dec, Radius);

  if (skylist[0].Nregions > 1) {
    gprint (GP_ERR, "warning, radius overlaps region boundary, not yet implemented\n");
  }

  /* lock, load, unlock catalog */
  catalog.filename = skylist[0].filename[0];
  catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
  catalog.Nsecfilt = 0;

  // an error exit status here is a significant error
  if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
      fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename);
      exit (2);
  }
  dvo_catalog_unlock (&catalog);

  Nstars = catalog.Naverage;
  ALLOCATE (RA, double, Nstars);
  ALLOCATE (DEC, double, Nstars);
  ALLOCATE (N1, int, Nstars);

  /* find star(s) in RA, DEC list -- use a dumb algorithm for now, improve later */
  /* stars are not guaranteed to be sorted in RA or in DEC, so first sort the list */
  for (i = 0; i < Nstars; i++) {
    RA[i] = catalog.average[i].R;
    DEC[i] = catalog.average[i].D;
    N1[i] = i;
  }
  /* sort list by DEC */
  if (Nstars > 1) sort_coords_index (DEC, RA, N1, Nstars);
  /* at this point, RA, DEC, and N1 are sorted by DEC.  
     catalog.average[N1[i]].R = RA[i] */

  /* bracket the RA range of interest */
  dec0 = Dec - Radius;
  dec1 = Dec + Radius;

  Nlo = 0; Nhi = catalog.Naverage;
  while (Nhi - Nlo > 10) {
    N = 0.5*(Nlo + Nhi);
    if (DEC[N] < dec0) {
      Nlo = N;
    } else {
      Nhi = N + 1;
    }
  }
  /* DEC[Nlo] is guaranteed to be just lower than dec0 */

  Radius2 = Radius*Radius;
  found = FALSE;

  /* data has been loaded, get ready to plot it */
  if (SaveVectors) {
    N = 0;
    NPTS = 1000;
    ALLOCATE (vec1[0].elements, float, NPTS);
    ALLOCATE (vec2[0].elements, float, NPTS);
    ALLOCATE (vec3[0].elements, float, NPTS);
    ALLOCATE (vec4[0].elements, float, NPTS);
  }

  initPhotcodeSequence (Nsecfilt);

  for (i = Nlo; (i < catalog.Naverage) && !found; i++) {

    if (dec0 > DEC[i]) continue;
    if (dec1 < DEC[i]) found = TRUE;
    
    r = SQ(Dec - DEC[i]) + SQ(Ra - RA[i]);
    if (r < Radius2) {
      k = N1[i];
      if (!QUIET) {
	gprint (GP_LOG, "star: %d\n", k);
	gprint (GP_LOG, "%11.7f ", catalog.average[k].R);
	gprint (GP_LOG, "%11.7f ", catalog.average[k].D);
	gprint (GP_LOG, "%3d   ",  catalog.average[k].Nmeasure);
	gprint (GP_LOG, "%4.1f ",  0.01*catalog.average[k].Xp);
	gprint (GP_LOG, "%5d",     catalog.average[k].code);

	if (FULL_OUTPUT) {
	    gprint (GP_LOG, "%f",     catalog.average[k].dR);
	    gprint (GP_LOG, "%f",     catalog.average[k].dD);
	    gprint (GP_LOG, "%f",     catalog.average[k].uR);
	    gprint (GP_LOG, "%f",     catalog.average[k].uD);
	    gprint (GP_LOG, "%f",     catalog.average[k].duR);
	    gprint (GP_LOG, "%f",     catalog.average[k].duD);
	    gprint (GP_LOG, "%f",     catalog.average[k].P);
	    gprint (GP_LOG, "%f",     catalog.average[k].dP);
	    gprint (GP_LOG, "%x",     catalog.average[k].objID);
	    gprint (GP_LOG, "%x",     catalog.average[k].catID);
	}

	gprint (GP_LOG, "\n");
      
	/* filter names */
	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 3);
	gprint (GP_LOG, "\n");

	/* average mags */
	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 0);
	gprint (GP_LOG, "\n");

	/* average mag errors */
	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 1);
	gprint (GP_LOG, "\n");

	/* average mag chisq */
	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 2);
	gprint (GP_LOG, "\n");
      }

      if (GetMeasures || SaveVectors) {
	m = catalog.average[k].measureOffset;
	for (j = 0; j < catalog.average[k].Nmeasure; j++, m++) {

	  Mcat = PhotCat (&catalog.measure[m]);
	  Mrel = PhotRel (&catalog.measure[m], &catalog.average[k], &catalog.secfilt[k*Nsecfilt]);

	  if (GetMeasures && !QUIET) {
	    date = ohana_sec_to_date (catalog.measure[m].t);
	    gprint (GP_LOG, "%6.3f ",  Mcat);
	    gprint (GP_LOG, "%6.3f ",  Mrel);
	    gprint (GP_LOG, "%5.3f  ", catalog.measure[m].dM);
	    gprint (GP_LOG, "%20s  ",  date);
	    gprint (GP_LOG, "%7.4f ",  catalog.measure[m].dR);
	    gprint (GP_LOG, "%7.4f",   catalog.measure[m].dD);
	    gprint (GP_LOG, "%2d ",    catalog.measure[m].dophot);
	    gprint (GP_LOG, "%3x ",    catalog.measure[m].dbFlags);
	    gprint (GP_LOG, "%3d ",    catalog.measure[m].photcode);
	    gprint (GP_LOG, "%-20s  ", GetPhotcodeNamebyCode (catalog.measure[m].photcode));

	    if (FULL_OUTPUT) {
		gprint (GP_LOG, "%f", catalog.measure[m].Mcal);
		gprint (GP_LOG, "%f", catalog.measure[m].Map);
		gprint (GP_LOG, "%f", pow(10.0, 0.4*catalog.measure[m].dt));
		gprint (GP_LOG, "%f", 1.0 + catalog.measure[m].airmass);
		gprint (GP_LOG, "%f", catalog.measure[m].az);
		gprint (GP_LOG, "%f", catalog.measure[m].Xccd);
		gprint (GP_LOG, "%f", catalog.measure[m].Yccd);
		gprint (GP_LOG, "%d", catalog.measure[m].dXccd);
		gprint (GP_LOG, "%d", catalog.measure[m].dYccd);
		gprint (GP_LOG, "%f", catalog.measure[m].Sky);
		gprint (GP_LOG, "%f", catalog.measure[m].dSky);
		gprint (GP_LOG, "%d", catalog.measure[m].averef);
		gprint (GP_LOG, "%d", catalog.measure[m].detID);
		gprint (GP_LOG, "%d", catalog.measure[m].imageID);
		gprint (GP_LOG, "%f", catalog.measure[m].qPSF);
		gprint (GP_LOG, "%f", catalog.measure[m].psfChisq);
		gprint (GP_LOG, "%f", catalog.measure[m].crNsigma);
		gprint (GP_LOG, "%f", catalog.measure[m].extNsigma);
		gprint (GP_LOG, "%f", 0.01*catalog.measure[m].FWx);
		gprint (GP_LOG, "%f", 0.01*catalog.measure[m].FWy);
		gprint (GP_LOG, "%f", (360.0/(float)0xffff)*catalog.measure[m].theta);

		gprint (GP_LOG, "%x", catalog.measure[m].photFlags);
		gprint (GP_LOG, "%d", catalog.measure[m].stargal);
	    }
	    gprint (GP_LOG, "\n");

	    free (date);
	  }
	  
	  if (SaveVectors) {
	    vec1[0].elements[N] = Mcat;
	    vec2[0].elements[N] = TimeValue (catalog.measure[m].t, TimeReference, TimeFormat);
	    vec3[0].elements[N] = catalog.measure[m].airmass;
	    vec4[0].elements[N] = catalog.measure[m].photcode;
	    N ++;
	    if (N == NPTS - 1) {
	      NPTS += 2000;
	      REALLOCATE (vec1[0].elements, float, NPTS);
	      REALLOCATE (vec2[0].elements, float, NPTS);
	      REALLOCATE (vec3[0].elements, float, NPTS);
	      REALLOCATE (vec4[0].elements, float, NPTS);
	    }
	  }
	}
      }
    }
  }

  if (SaveVectors) {
    vec1[0].Nelements = N;
    vec2[0].Nelements = N;
    vec3[0].Nelements = N;
    vec4[0].Nelements = N;
  }

  free (RA);
  free (DEC);
  free (N1);
  dvo_catalog_free (&catalog);

  freePhotcodeSequence ();
  return (TRUE);

}

void print_double (double value) {
  if (isnan(value)) 
    gprint (GP_LOG, "NaN    ");
  else 
    gprint (GP_LOG, "%6.3f ", value);
}

void print_short (double value, short int ival) {
  if (ival == NAN_S_SHORT) 
    gprint (GP_LOG, "NaN    ");
  else 
    gprint (GP_LOG, "%6.3f ", value);
}

// XXX fix printing to be in photcode numerical order for PRI/SEC data
static int *sequence = NULL;

void initPhotcodeSequence (int Nsecfilt) {

  int j;
  int *codeNumber;
  PhotCode *code;

  // sequence contains, in desired order, secfilt number (0 == pri)
  ALLOCATE (sequence, int, Nsecfilt);
  ALLOCATE (codeNumber, int, Nsecfilt);
  
  /* filter names -- primary code is 0 in this function */
  for (j = 0; j < Nsecfilt; j++) {
    code = GetPhotcodebyNsec (j);
    codeNumber[j] = code[0].code;
    sequence[j] = j;
  }

  isortpair (codeNumber, sequence, Nsecfilt);
  free (codeNumber);
}

void freePhotcodeSequence () {
  free (sequence);
}


void printPhotcodeSequence (Average *average, SecFilt *secfilt, int entry, int type) {

  int seq;
  PhotCode *code;

  seq = sequence[entry];

  switch (type) {
    case 0: /* average mags */
      if (seq == -1) {
	print_double (NAN);
      } else {
	print_double (secfilt[seq].M);
      }
      break;

    case 1: /* average mags errors */
      if (seq == -1) {
	print_double (NAN);
      } else {
	print_double (secfilt[seq].dM);
      }
      break;

    case 2: /* average mag chisq */
      if (seq == -1) {
	print_short (NAN_S_SHORT, NAN_S_SHORT);
      } else {
	print_short (pow (10.0, 0.01*secfilt[seq].Xm), secfilt[seq].Xm);
      }
      break;

    case 3: /* filter names */
      code = GetPhotcodebyNsec (seq);
      gprint (GP_LOG, "%-6s ", code[0].name);
      break;
  }
}
