# include "dvo.h"

void print_value (FILE *f, double value, short int ival);

int gstar (int argc, char **argv) {
  
  char filename[128], catdir[256], *date;
  double Ra, Dec, Radius, Radius2, r, dec0, dec1;
  double t, Mcat, Mrel;
  float *RA, *DEC;
  int i, j, k, m, N, *N1, Nregions, Nsec, NPTS;
  int Nstars, found, GetMeasures, Nlo, Nhi;
  int SaveVectors;
  Vector *vec1, *vec2, *vec3, *vec4;
  RegionFile *regions;
  Catalog catalog;
  PhotCode *code;

  VarConfig ("CATDIR", "%s", catdir);
  if (!InitPhotcodes ()) return (FALSE);
  Nsec = GetPhotcodeNsecfilt ();

  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) {
    fprintf (stderr, "USAGE: gstar RA DEC Radius [-m]\n");
    return (FALSE);
  }
  
  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;
  
  regions = find_regions (Ra, Dec, Radius, &Nregions);

  if (Nregions > 1) {
    fprintf (stderr, "warning, radius overlaps region boundary, not yet implemented\n");
  }

  /* set filename, read in header */
  sprintf (filename, "%s/%s", catdir, regions[0].name);
  /* lock, load, unlock catalog */
  catalog.filename = filename;
  switch (lock_catalog (&catalog, LCK_SOFT)) {
  case 2:
    unlock_catalog (&catalog);
  case 0:
    return (FALSE);
  }
  if (!load_catalog (&catalog, LOAD_AVES | LOAD_MEAS | LOAD_SECF, TRUE)) {
    unlock_catalog (&catalog);
    return (FALSE);
  }
  unlock_catalog (&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);
  }

  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];
      fprintf (GetOutfile(), "star: %d\n", k);
      fprintf (GetOutfile(), "%9.5f %9.5f %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 -- primary code is 0 in this function */
      for (j = 0; j < Nsec + 1; j++) {
	code = GetPhotcodebyNsec (j);
	fprintf (GetOutfile (), "%s ", code[0].name);
      }
      fprintf (GetOutfile (), "\n");

      /* average mags */
      print_value (GetOutfile(), 0.001*catalog.average[k].M, catalog.average[k].M);
      for (j = 0; j < Nsec; j++) print_value (GetOutfile(), 0.001*catalog.secfilt[j + Nsec*k].M, catalog.secfilt[j + Nsec*k].M);
      fprintf (GetOutfile(), "\n");

      /* average mag errors */
      print_value (GetOutfile(), 0.001*catalog.average[k].dM, catalog.average[k].dM);
      for (j = 0; j < Nsec; j++) print_value (GetOutfile(), 0.001*catalog.secfilt[j + Nsec*k].dM, catalog.secfilt[j + Nsec*k].dM);
      fprintf (GetOutfile(), "\n");

      /* average mag chisq */
      print_value (GetOutfile(), pow (10.0, 0.01*catalog.average[k].Xm), catalog.average[k].Xm);
      for (j = 0; j < Nsec; j++) print_value (GetOutfile(), pow (10.0, 0.01*catalog.secfilt[j + Nsec*k].Xm), catalog.secfilt[j + Nsec*k].Xm);
      fprintf (GetOutfile(), "\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*Nsec]);

	  if (GetMeasures) {
	    date = sec_to_date (catalog.measure[m].t);
	    fprintf (GetOutfile(), "%6.3f %6.3f %5.3f  %20s  %5.2f %5.2f %2d %3x %3d %-20s\n", 
		     Mcat, Mrel, 0.001*catalog.measure[m].dM,
		     date, 0.01*catalog.measure[m].dR, 0.01*catalog.measure[m].dD,
		     catalog.measure[m].dophot, catalog.measure[m].flags,
		     catalog.measure[m].source, GetPhotcodeNamebyCode (catalog.measure[m].source));
	    free (date);
	  }
	  
	  if (SaveVectors) {
	    vec1[0].elements[N] = Mcat;
	    vec2[0].elements[N] = t;
	    vec3[0].elements[N] = 0.001*catalog.measure[m].airmass;
	    vec4[0].elements[N] = catalog.measure[m].source;
	    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);
  if (catalog.average != 0) free (catalog.average);
  if (catalog.measure != 0) free (catalog.measure);
  if (catalog.secfilt != 0) free (catalog.secfilt);

  return (TRUE);

}


void print_value (FILE *f, double value, short int ival) {
  if (ival == NO_MAG) 
    fprintf (f, "NaN    ");
  else 
    fprintf (f, "%6.3f ", value);
}
