# include "gastro2.h"
# define MAX_NAME 256

static int Nregion, NREGION;
static CatStats *region;
static char *CatDir;

/* data which defines the layout of the regions */
# 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 patch on sky with RA[0] - RA[1] & DEC[0] - DEC[1], return all that fall in range */
CatStats *gregions (CatStats *patch, int *nregion) {
  
  CatStats area;
  double dtmp;
  char file[MAX_NAME];

  init_regions ();

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

  /* include the pole if dec > 86.25 is included (ADD SOUTH POLE!)  */
  if (patch[0].DEC[1] > 86.25) {
    area.RA[0]  =   0.00;
    area.RA[1]  = 360.00;
    area.DEC[0] =  86.25;
    area.DEC[1] =  90.00;
    sprintf (file, "n8230/pole");
    area.name = strcreate (file);
    add_to_regions (&area);
    free (area.name);
  }
  
  /* if RA[0] > RA[1], split in two passes */
  if (patch[0].RA[0] > patch[0].RA[1]) {
    area = patch[0];
    area.RA[0]  = 0.0;
    find_dec_bands (&area);

    area = patch[0];
    area.RA[1]  = 360.0;
    find_dec_bands (&area);
  } else {
    area = patch[0];
    find_dec_bands (&area);
  }

  *nregion = Nregion;
  return (region);
}

/* 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 */

int find_dec_bands (CatStats *area) {

  int i;

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

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

/* 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] */
int load_ra_blocks (int Ndec, CatStats *area) {

  CatStats tregion;
  int i, Nlines, Nskip;
  char *buffer, file[MAX_NAME];
  FILE *f;

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

  ALLOCATE (buffer, char, 48*Nlines);
  ALLOCATE (tregion.name, char, MAX_NAME);

  /* get data from GSC file */ 
  f = fopen (GSCFILE, "r");
  if (f == NULL) {
    fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE);
    exit (1);
  }
  fseek (f, 5*2880 + 48*Nskip, SEEK_SET);
  fread (buffer, 48*Nlines, 1, f);
  fclose (f);

  /* find overlapping regions */
  for (i = 0; i < Nlines; i++) {
    parse_GSC_line (&tregion, &buffer[i*48]);
    if (area[0].RA[0] > tregion.RA[1]) continue;
    if (area[0].RA[1] < tregion.RA[0]) continue;
    if (area[0].DEC[0] > tregion.DEC[1]) continue;
    if (area[0].DEC[1] < tregion.DEC[0]) continue;
    sprintf (file, "%s/%s", Dec2Sections[Ndec], tregion.name);
    strcpy (tregion.name, file);
    add_to_regions (&tregion);
  }

  free (buffer);
  free (tregion.name);
  return (TRUE);
}

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

  /* copy name (4 chars) */
  bzero (tregion[0].name, 10);
  strncpy (tregion[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;
  tregion[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;
  tregion[0].RA[1] = v0 * 15.0;

  /* this only happens at the 0,360 boundary */
  if (tregion[0].RA[1] < tregion[0].RA[0]) tregion[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;
  tregion[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;
  tregion[0].DEC[1] = s * v0;

  if (tregion[0].DEC[0] > tregion[0].DEC[1]) {
    dtmp = tregion[0].DEC[1];
    tregion[0].DEC[1] = tregion[0].DEC[0];
    tregion[0].DEC[0] = dtmp;
  }
  return (TRUE);
}

void init_regions () {

  Nregion = 0;
  NREGION = 100;

  ALLOCATE (region, CatStats, NREGION);

}

void add_to_regions (CatStats *area) {

  region[Nregion] = area[0];
  ALLOCATE (region[Nregion].name, char, MAX_NAME);
  sprintf (region[Nregion].name, "%s/%s.%s", CatDir, area[0].name, "cpt");

  area_of_region (&region[Nregion]);

  Nregion ++;

  if (Nregion == NREGION) {
    NREGION += 100;
    REALLOCATE (region, CatStats, NREGION);
  }

}

void set_catalog (char *catdir) {
  CatDir = catdir;
}
