# include "dvoshell.h"

int lightcurve (int argc, char **argv) {
  
  double Ra, Dec, Radius, Radius2, r;
  float *RA, *DEC;
  int Nstars, found, PhotCodeSelect;
  int i, j, k, m, N, NPTS, Nsecfilt, RELPHOT, *N1, TimeFormat;
  time_t TimeReference;

  PhotCode *code;
  Catalog catalog;
  SkyTable *sky;
  SkyList *skylist;
  Vector *tvec, *mvec, *dmvec;

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

  if ((tvec = SelectVector ("tc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((mvec = SelectVector ("mc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((dmvec = SelectVector ("dmc", ANYVECTOR, TRUE)) == NULL) return (FALSE);

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

  code = NULL;
  PhotCodeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-photcode"))) {
    remove_argument (N, &argc, argv);
    if ((code = GetPhotcodebyName (argv[N])) == NULL) {
      gprint (GP_ERR, "ERROR: photcode not found in photcode table\n");
      return (FALSE);
    }
    PhotCodeSelect = TRUE;
    remove_argument (N, &argc, argv);
  }

  if (argc < 4) {
    gprint (GP_ERR, "USAGE: lightcurve RA DEC Radius\n");
    return (FALSE);
  }
  
  Ra = atof (argv[1]);
  Dec = atof (argv[2]);
  Radius = atof (argv[3]);

  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");
  }

  /* set filename, read in header */
  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] */

  N = 0;
  NPTS = 100;
  REALLOCATE (tvec[0].elements, float, NPTS);
  REALLOCATE (mvec[0].elements, float, NPTS);
  REALLOCATE (dmvec[0].elements, float, NPTS);
  
  GetTimeFormat (&TimeReference, &TimeFormat);

  Radius2 = Radius*Radius;
  found = FALSE;
  for (i = 0; (i < catalog.Naverage) && !found; i++) {

    /* this can be improved by using a couple of jumps to get within range */
    if (Dec > DEC[i] + Radius)
      continue;
    
    r = SQ(Dec - DEC[i]) + SQ(Ra - RA[i]);
    if (r < Radius2) {
      k = N1[i];
      /* found star, extract measurements */
      m = catalog.average[k].offset;
      for (j = 0; j < catalog.average[k].Nm; j++, m++) {

	if (PhotCodeSelect) {
	  if ((code[0].type == PHOT_REF) || (code[0].type == PHOT_DEP)) {
	    if (code[0].code != catalog.measure[m].photcode) continue;
	  } 
	  if (code[0].type == PHOT_SEC) {
	    if (code[0].code != GetPhotcodeEquivCodebyCode (catalog.measure[m].photcode)) continue;
	  } 
	}      

	tvec[0].elements[N] = TimeValue (catalog.measure[m].t, TimeReference, TimeFormat);
	dmvec[0].elements[N] = catalog.measure[m].dM;
	if (RELPHOT) {
	  mvec[0].elements[N] = PhotCat (&catalog.measure[m]);
	} else {
	  mvec[0].elements[N] = PhotRel (&catalog.measure[m], &catalog.average[k], &catalog.secfilt[k*Nsecfilt]);
	}
	N++; 
	if (N == NPTS) {
	  NPTS += 100;
	  REALLOCATE (tvec[0].elements, float, NPTS);
	  REALLOCATE (mvec[0].elements, float, NPTS);
	  REALLOCATE (dmvec[0].elements, float, NPTS);
	}
      }      
    }
  }
  sortthree (tvec[0].elements, mvec[0].elements, dmvec[0].elements, N);
  tvec[0].Nelements = mvec[0].Nelements = dmvec[0].Nelements = N;

  free (RA);
  free (DEC);
  free (N1);
  dvo_catalog_free (&catalog);
  SkyListFree (skylist, FALSE);
  return (TRUE);
}
