# include "gastro.h"

# define NDECBANDS 12
double DecBand[] = {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,
		    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};

int NDecLines[] =  {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, 0, 
		    597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48, 0};

int NRaBlocks [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0, 
		    48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0};

char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", 
			"n5230", "n6000", "n6730", "n7500", "n8230", "weirdness", 
			"s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", 
			"s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"};

/* given RA[0] - RA[1] & DEC[0] - DEC[1], return all that fall in range */
char **gregions (CatStats *catstats, int *nnames) {
  
  CatalogRegion area;
  double dtmp;
  int i, Nnames;
  char **names, file[256];

  Nnames = 0;
  ALLOCATE (names, char *, 1);

  /* force RA[i] to be in range 0 - 360 */
  while (catstats[0].RA[0] > 360.0) { catstats[0].RA[0] -= 360.0; }
  while (catstats[0].RA[0] < 0.0)   { catstats[0].RA[0] += 360.0; }
     
  while (catstats[0].RA[1] > 360.0) { catstats[0].RA[1] -= 360.0; }
  while (catstats[0].RA[1] < 0.0)   { catstats[0].RA[1] += 360.0; }
     
  /* force DEC[0] < DEC[1] */
  if (catstats[0].DEC[0] > catstats[0].DEC[1]) {
    dtmp = catstats[0].DEC[0];
    catstats[0].DEC[0] = catstats[0].DEC[1];
    catstats[0].DEC[1] = dtmp;
  }

  /* include the pole if dec > 86.25 is included (ADD SOUTH POLE!)  */
  if (catstats[0].DEC[1] > 86.25) {
    sprintf (file, "n8230/pole");
    names[Nnames] = strcreate (file);
    Nnames ++;
  }
  
  /* if RA[0] > RA[1], split in two passes */
  if (catstats[0].RA[0] > catstats[0].RA[1]) {
    area.RA[0]  = 0.0;
    area.RA[1]  = catstats[0].RA[1];
    area.DEC[0] = catstats[0].DEC[0]; 
    area.DEC[1] = catstats[0].DEC[1];
    names = find_dec_bands (area, names, &Nnames);

    area.RA[0]  = catstats[0].RA[0];
    area.RA[1]  = 360.0;
    area.DEC[0] = catstats[0].DEC[0]; 
    area.DEC[1] = catstats[0].DEC[1];
    names = find_dec_bands (area, names, &Nnames);

  } else {

    area.RA[0]  = catstats[0].RA[0];
    area.RA[1]  = catstats[0].RA[1];
    area.DEC[0] = catstats[0].DEC[0]; 
    area.DEC[1] = catstats[0].DEC[1];
    names = find_dec_bands (area, names, &Nnames);
  }

  /* add the full path to the short names */
  for (i = 0; i < Nnames; i++) {
    sprintf (file, "%s/%s.cpt", GSC_DIR, names[i]);
    free (names[i]);
    names[i] = strcreate (file);
  }

  *nnames = Nnames;
  return (names);

}

/* given an area [RA[0] - RA[1], DEC[0] - DEC[1],
   where 0.0 < RA[0] < RA[1] < 360.0 
   and -90.0 < (DEC[0], DEC[1]) < 90.0 
   and DEC[0] < DEC[1]
   search through the dec bands to identify ones which overlap the area */

/* must pass in an allocted pointer names, and Nnames contains the number so far */
char **find_dec_bands (CatalogRegion area, char **names, int *Nnames) {

  int i;

  /* for NORTH region (DecBand[i] < DecBand[i+1]) */
  for (i = 0; i < NDECBANDS; i++) {
    if ((area.DEC[1] > DecBand[i]) || (area.DEC[0] < DecBand[i+1])) {
      names = load_ra_blocks (i, area, names, Nnames);
    }
  }

  /* for SOUTH region (DecBand[i] > DecBand[i+1]) */
  for (i = NDECBANDS+1; i < 2*NDECBANDS; i++) {
    if ((area.DEC[0] < DecBand[i]) || (area.DEC[1] > DecBand[i+1])) {
      names = load_ra_blocks (i, area, names, Nnames);
    }
  }

  return (names);
}


/* given the number of a DecBlock and a defined area, return all area names in that block */
/* area is defined so that RA[0] < RA[1], DEC[0] < DEC[1] */
char **load_ra_blocks (int Ndec, CatalogRegion area, char **names, int *Ntnames) {

  CatalogRegion region;
  int i, Nlines, Nnames;
  char *buffer, file[256];
  FILE *f;

  /* open file */ 
  f = fopen (GSCFILE, "r");
  if (f == NULL) {
    fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE);
    exit (0);
  }

  /* count lines before Dec Band */
  Nlines = 0;
  for (i = 0; i < Ndec; i++) Nlines += NDecLines[i];

  /* seek to start of Dec Band */
  fseek (f, 5*2880 + 48*Nlines, SEEK_SET);
  
  /* load data into buffer */
  Nlines = NDecLines[Ndec];
  ALLOCATE (buffer, char, 48*Nlines);
  fread (buffer, 48*Nlines, 1, f);
  fclose (f);

  /* must pass in a pointer which is already allocated */
  Nnames = *Ntnames;
  REALLOCATE (names, char *, MAX (Nnames + Nlines, 1));

  for (i = 0; i < Nlines; i++) {
    parse_GSC_line (&region, &buffer[i*48]);
    if (area.RA[0] > region.RA[1]) continue;
    if (area.RA[1] < region.RA[0]) continue;
    if (area.DEC[0] > region.DEC[1]) continue;
    if (area.DEC[1] < region.DEC[0]) continue;
    sprintf (file, "%s/%s", Dec2Sections[Ndec], region.name);
    names[Nnames] = strcreate (file);
    Nnames ++;
  }

  free (buffer);
  REALLOCATE (names, char *, MAX (Nnames, 1));
  *Ntnames = Nnames;
  return (names);

}

/* parse the GSC line, placing values in 'region'.
   force RA[0] < RA[1], DEC[0] < DEC[1] */
int parse_GSC_line (CatalogRegion *region, char *line) {
  
  double v0, s, dtmp;
  char tmp[10];

  /* copy name (4 chars) */
  bzero (region[0].name, 10);
  strncpy (region[0].name, &line[1], 4);

  /* get RA[0] */
  bzero (tmp, 10);
  strncpy (tmp, &line[7], 2);
  v0 = atof (tmp);
  
  bzero (tmp, 10);
  strncpy (tmp, &line[10], 2);
  v0 += atof (tmp) / 60.0;
  
  bzero (tmp, 10);
  strncpy (tmp, &line[13], 5);
  v0 += atof (tmp) / 3600.0;
  region[0].RA[0] = v0 * 15.0;

  /* get RA[1] */
  bzero (tmp, 10);
  strncpy (tmp, &line[19], 2);
  v0 = atof (tmp);
  
  bzero (tmp, 10);
  strncpy (tmp, &line[22], 2);
  v0 += atof (tmp) / 60.0;
  
  bzero (tmp, 10);
  strncpy (tmp, &line[25], 5);
  v0 += atof (tmp) / 3600.0;
  region[0].RA[1] = v0 * 15.0;

  /* this only happens at the 0,360 boundary */
  if (region[0].RA[1] < region[0].RA[0]) region[0].RA[1] += 360.0;

  /* get DEC[0] */
  if (line[31] == '-') s = -1;
  else s = +1;

  bzero (tmp, 10);
  strncpy (tmp, &line[32], 2);
  v0 = atof (tmp);
  
  bzero (tmp, 10);
  strncpy (tmp, &line[35], 2);
  v0 += atof (tmp) / 60.0;
  region[0].DEC[0] = s * v0;

  /* get DEC[1] */
  if (line[40] == '-') s = -1;
  else s = +1;

  bzero (tmp, 10);
  strncpy (tmp, &line[41], 2);
  v0 = atof (tmp);
  
  bzero (tmp, 10);
  strncpy (tmp, &line[44], 2);
  v0 += atof (tmp) / 60.0;
  region[0].DEC[1] = s * v0;

  if (region[0].DEC[0] > region[0].DEC[1]) {
    dtmp = region[0].DEC[1];
    region[0].DEC[1] = region[0].DEC[0];
    region[0].DEC[0] = dtmp;
  }

  return (TRUE);

}

