# include "gastro2.h"

# define NDECLIST 24
double DecList[] = {-90.0, -82.5, -75.0, -67.5, -60.0, -52.5, -45.0, -37.5, -30.0, -22.5, -15.0, -7.5, 0.0, 
		    +7.5, +15.0, +22.5, +30.0, +37.5, +45.0, +52.5, +60.0, +67.5, +75.0, +82.5, +90.0};

typedef struct {
  char name[256];
  float rmin, rmax;
  int Nstar, Nbyte, Ncbyte;
} TwoMassAcc;

int main (int argc, char **argv) {

  int i, j, k, Nregions;
  TwoMassAcc acc[49];
  CatStats patch, *regions;
  char line[1024], filename[256], saveline[1024], *path, *Fullpath, *Basename, *Pathname;
  FILE *f, *fout;
  int status, done, Nacc, Nstar, NSTAR, Nout, Nscan;
  double *ra, *dec, *mag, DecMin, DecMax;

  ConfigInit (&argc, argv);

  if (argc != 5) {
    fprintf (stderr, "USAGE: extr2mass (input) (acc) (outdir) (gsctbl)\n");
    exit (2);
  }
  
  /* load accelerator file */
  f = fopen (argv[2], "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "can't load accelerator file %s\n", argv[2]);
    return (FALSE);
  }
  for (i = 0; i < 49; i++) {
    status = fscanf (f, "%s %f %f %d %d %d", acc[i].name, &acc[i].rmin, &acc[i].rmax,
	    &acc[i].Nstar, &acc[i].Nbyte, &acc[i].Ncbyte);
    if (status != 6) { 
      fprintf (stderr, "ERROR: can't read 2Mass accelerator file\n");
      return (FALSE);
    }
    acc[i].rmin *= 15.0;
    acc[i].rmax *= 15.0;
  }
  fclose (f);

  /* find file in acc file */
  Nacc = -1;
  for (i = 0; (Nacc == -1) && (i < 49); i++) {
    if (!strcmp (argv[1], acc[i].name)) Nacc = i; 
  }
  if (Nacc == -1) {
    fprintf (stderr, "can't find file in acc file\n");
    exit (1);
  }

  fprintf (stderr, "reading from %s, %d stars\n", acc[Nacc].name, acc[Nacc].Nstar);

  /* load 2mass data file */
  /*
  sprintf (line, "gzip -cd %s > %s", acc[Nacc].name, acc[Nacc].name);
  system (line);
  */
  
  sprintf (line, "gzip -cd %s", acc[Nacc].name);
  f = popen (line, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "ERROR: can't read 2Mass data file: %s\n", line);
    return (FALSE);
  }

  /* all data in 2mass file are sequential in DEC.  read stars while in dec range */

  /* load data in DEC blocks */
  saveline[0] = 0;
  for (i = 0; i < NDECLIST; i++) {
    
    Nstar = 0;
    NSTAR = 1000;
    ALLOCATE (ra, double, NSTAR);
    ALLOCATE (dec, double, NSTAR);
    ALLOCATE (mag, double, NSTAR);
    
    DecMin = DecList[i];
    DecMax = DecList[i+1];

    done = FALSE;
    while (!done) {
      if (saveline[0]) {
	strcpy (line, saveline);
	saveline[0] = 0;
      } else {
	if (scan_line (f, line) == EOF) { goto escape; }
      }
      Nscan = sscanf (line, "%lf %lf %*s %*s %*s %*s %lf", &ra[Nstar], &dec[Nstar], &mag[Nstar]);
      if (Nscan != 3) { 
	fprintf (stderr, "read wrong number of entries: %d\n", Nscan);
	fprintf (stderr, "%s\n", line);
	exit (1);
      }
      if (dec[Nstar] > DecMax) {
	strcpy (saveline, line);
	goto escape; 
      } 
      Nstar ++;
      if (Nstar == NSTAR) {
	NSTAR += 1000;
	REALLOCATE (ra, double, NSTAR);
	REALLOCATE (dec, double, NSTAR);
	REALLOCATE (mag, double, NSTAR);
      }
      if (Nstar > acc[Nacc].Nstar) {
	fprintf (stderr, "problem!\n %s\n", line);
	exit (1);
      }
    }
  escape:
      
    fprintf (stderr, "read %d stars for DEC block %f %f\n", Nstar, DecMin, DecMax);

    /* now we write the Nstar objects to the appropriate files */
    patch.RA[0] = acc[Nacc].rmin;
    patch.RA[1] = acc[Nacc].rmax;
    patch.DEC[0] = DecMin;
    patch.DEC[1] = DecMax;
    
    regions = gregions (&patch, &Nregions);
    
    for (j = 0; j < Nregions; j++) {
      /* name has complete path to gsc, strip off path - 1 */
      Fullpath = pathname (regions[j].name);
      Basename = filebasename (regions[j].name);
      Pathname = filebasename (Fullpath);
      sprintf (filename, "%s/%s/%s", argv[3], Pathname, Basename);

      path = pathname (filename);
      mkdirhier (path);
      
      fout = fopen (filename, "a+");
      if (fout == (FILE *) NULL) {
	fprintf (stderr, "can't open output file %s\n", filename);
	exit (1);
      }

      Nout = 0;
      for (k = 0; k < Nstar; k++) {
	if (ra[k] < regions[j].RA[0]) continue;
	if (ra[k] >= regions[j].RA[1]) continue;
	if (dec[k] < regions[j].DEC[0]) continue;
	if (dec[k] >= regions[j].DEC[1]) continue;
	fprintf (fout, "%8.4f %8.4f %4.1f\n", ra[k], dec[k], mag[k]);
	Nout ++;
      }
      if (Nout > 0) {
	fprintf (stderr, "wrote %d stars to %s\n", Nout, filename);
      }
      
      free (path);
      free (Fullpath);
      free (Basename);
      free (Pathname);
      fclose (fout);
    }

    free (ra);
    free (dec);
    free (mag);
  }

  fclose (f);
  exit (0);
}


void area_of_region (CatStats *region) {
  
  double area;

  area = DEG_RAD*(region[0].RA[1] - region[0].RA[0])*(sin(region[0].DEC[1]*RAD_DEG) - sin(region[0].DEC[0]*RAD_DEG));
  region[0].Area = area;
}

