# 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;
  float *RA, *DEC;
  int i, j, k, m, N, *N1, Nsecfilt, NPTS, QUIET;
  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);
  }

  GetMeasures = FALSE;
  if ((N = get_argument (argc, argv, "-m"))) {
    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, float, Nstars);
  ALLOCATE (DEC, float, 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_lists (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 %11.7f %3d of %3d  %4.1f %5d\n", catalog.average[k].R, catalog.average[k].D, 
		 catalog.average[k].Nm, catalog.average[k].Nn + catalog.average[k].Nm,
		 0.01*catalog.average[k].Xp, catalog.average[k].code);
      
	/* 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].offset;
	for (j = 0; j < catalog.average[k].Nm; j++, m++) {

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

	  if (GetMeasures && !QUIET) {
	    date = sec_to_date (catalog.measure[m].t);
	    gprint (GP_LOG, "%6.3f %6.3f %5.3f  %20s  %5.2f %5.2f %2d %3x %3d %-20s\n", 
		     Mcat, Mrel, catalog.measure[m].dM,
		     date, catalog.measure[m].dR, catalog.measure[m].dD,
		     catalog.measure[m].dophot, catalog.measure[m].flags,
		     catalog.measure[m].photcode, GetPhotcodeNamebyCode (catalog.measure[m].photcode));
	    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_value (double value, short int ival) {
  if (ival == NO_MAG) 
    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;
  }

  isort_pair (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_value (NO_MAG, NO_MAG);
      } else {
	print_value (secfilt[seq].M, secfilt[seq].M);
      }
      break;

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

    case 2: /* average mag chisq */
      if (seq == -1) {
	print_value (NO_MAG, NO_MAG);
      } else {
	print_value (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;
  }
}
