# include "gastro.h"
# define dcos(a) (cos((double)((a)*(RAD_DEG))))
# define dsin(a) (sin((double)((a)*(RAD_DEG))))

/* by necesity hard wired */
# define D_NSTARS 1000
# define BYTES_STAR 31
# define BLOCK 1000
# include <sys/time.h>
# include <time.h>

SStars *gstars (file, NSTARS, coords, NX, NY, dNdM) 
char file[];
int *NSTARS, *NX, *NY;
double *dNdM;
Coords *coords;
{

  char line[64], *p, side[64];
  Header header, theader;
  FILE *f;
  int i, j, sign, Ninstar, nstar, rnumber, N, Nstars, nbytes, Nbytes, status, Nccd;
  SStars *stars;
  char *buffer;
  double X, Y, T1, T2, T3, type, csign;
  double PD, PR, DE, RE;
  double ra, dec;

  /* read in image header, open image data region */
  if (!fits_read_header (file, &header)) {
    fprintf (stderr, "ERROR: can't find image file %s (1)\n", file);
    exit(0);
  }
  /* get complete info from header */
  fits_scan (&header, "NAXIS1", "%d", 1, NX);
  fits_scan (&header, "NAXIS2", "%d", 1, NY);

  if (!strcasecmp (ROUGH_ASTROMETRY, "header")) {
    if (!HEADER[0]) {
      GetCoords (coords, &header);
    } else {
      fits_read_header (HEADER, &theader);
      GetCoords (coords, &theader);
    }
    if (!strcmp (coords[0].ctype, "NONE") || (coords[0].cdelt1 == 0) ||  (coords[0].cdelt2 == 0)) {
      fprintf (stderr, "header coordinates incomplete, trying for rough coordinates\n");
      strcpy (ROUGH_ASTROMETRY, "config");
    } else {
      if (FLIPX) {
	coords[0].pc1_1 *= -1;
	coords[0].crpix1 = *NX - coords[0].crpix1;
      }
      if (FLIPY) {
	coords[0].pc2_2 *= -1;
	coords[0].crpix2 = *NY - coords[0].crpix2;
      }
      ASEC_PIX = fabs (coords[0].cdelt1 * 3600.0);
      csign = coords[0].cdelt1 / fabs (coords[0].cdelt1);
      CCD_PC1_1 = coords[0].pc1_1 * csign;
      CCD_PC2_1 = coords[0].pc2_1 * csign;
      csign = coords[0].cdelt2 / fabs (coords[0].cdelt2);
      CCD_PC1_2 = coords[0].pc1_2 * csign;
      CCD_PC2_2 = coords[0].pc2_2 * csign;
      if (!strcmp (&coords[0].ctype[4], "-PLY")) {
	strcpy (coords[0].ctype, "DEC--TAN");
      }
    }
  }
  
  /* get partial info from header */
  if (!strcasecmp (ROUGH_ASTROMETRY, "config")) {
    /* default values for coords */
    strcpy (coords[0].ctype, "RA---TAN");
    coords[0].pc1_1 = CCD_PC1_1; coords[0].pc1_2 = CCD_PC1_2;
    coords[0].pc2_1 = CCD_PC2_1; coords[0].pc2_2 = CCD_PC2_2;
    coords[0].cdelt1 = coords[0].cdelt2 = ASEC_PIX / 3600.0;
    coords[0].Npolyterms = 0;
    coords[0].crpix1 = 0.5*(*NX);
    coords[0].crpix2 = 0.5*(*NY);
    /* RA */
    fits_scan (&header, "RA", "%s", 1, line);
    coords[0].crval1 = atof(&line[0]);
    p = strchr (line, ':');
    coords[0].crval1 += atof((char *) (p+1)) / 60.0;
    p = strchr (p+1, ':');
    coords[0].crval1 += atof((char *) (p+1)) / 3600.0;
    coords[0].crval1 = coords[0].crval1 * 15.0;
    /* DEC */
    fits_scan (&header, "DEC", "%s", 1, line);
    coords[0].crval2 = atof(&line[0]);
    sign = 1;
    if (line[0] == '-') {
      sign = -1;
      coords[0].crval2 = fabs(coords[0].crval2 );
    }
    p = strchr (line, ':');
    coords[0].crval2 += atof((char *) (p+1)) / 60.0;
    p = strchr (p+1, ':');
    coords[0].crval2 += atof((char *) (p+1)) / 3600.0;
    coords[0].crval2 *= sign;
  }
  if (VERBOSE) fprintf (stderr, "coordinates from header: %9.4f %9.4f\n", coords[0].crval1, coords[0].crval2);

  /* the following two sections are LONEOS derived and may not be needed elsewhere */
  if (LONEOS_COORDS) {
    fits_scan (&header, "COMMENT", "%s", 1, line);
    sscanf (line, "%*s%d%s", &rnumber, side);
    if (get_region_coords (&ra, &dec, rnumber, side)) {
      if (fabs(ra - coords[0].crval1) > 0.1) {
	fprintf (stderr, "large offset from claimed position, using region coords %f %f -> %f %f (%d %s)\n", 
		 coords[0].crval1, coords[0].crval2, ra, dec, rnumber, side);
	coords[0].crval1 = ra;
	coords[0].crval2 = dec;
      }
    }
  }

  /* at this point, we need to correct the crval1, crval2, and ROT_ZERO values
     based on the pole axis angle and the ra, dec offsets */
  if (POLAR_ALIGNMENT) {
    X = coords[0].crval1;
    Y = coords[0].crval2;
    PD = POLE_DEC;  PR = POLE_RA;
    DE = DEC_OFFSET; RE = RA_OFFSET;
    
    T1 = dcos(Y-DE) * dcos(X-RE) * dsin(PD) + dsin(Y-DE) * dcos(PD);
    T2 = dcos(Y-DE) * dsin(X-RE);
    T3 = dsin(Y-DE) * dsin(PD) - dcos(Y-DE) * dcos(X-RE) * dcos(PD);
    
    coords[0].crval1 = (DEG_RAD * atan2 (T2, T1)) + PR;
    coords[0].crval2 = (DEG_RAD * asin (T3));
    while (coords[0].crval1 < 0) coords[0].crval1 += 360.0;
    while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0;
    
    if (VERBOSE) fprintf (stderr, "  after polar alignment: %9.4f %9.4f\n", coords[0].crval1, coords[0].crval2);
  }

  if (FORCE) {
    coords[0].crval1 = F_RA;
    coords[0].crval2 = F_DEC;
    if (VERBOSE) fprintf (stderr, " forcing coordinates to: %9.4f %9.4f\n", coords[0].crval1, coords[0].crval2);
  }    

  Nstars = 0;
  fits_scan (&header, "NSTARS", "%d", 1, &Nstars);
  if (Nstars == -1) {
    fprintf (stderr, "ERROR: failed to find NSTARS\n");
    exit (0);
  }
  ALLOCATE (stars, SStars, Nstars);
  Nbytes = Nstars*BYTES_STAR;

  /* re-open file for stars */
  f = fopen (file, "r");
  if (f == NULL) {
    fprintf (stderr, "ERROR: can't find image file %s (2)\n", file);
    exit(0);
  }
  fseek (f, header.size, SEEK_SET); 

  N = nstar = 0;
  ALLOCATE (buffer, char, (BLOCK*BYTES_STAR));
  
  while ((nbytes = fread (buffer, 1, (BLOCK*BYTES_STAR), f)) != 0) {
    Ninstar = nbytes / BYTES_STAR;
    for (j = 0; j < Ninstar; j++, nstar++) {
      dparse (&stars[N].X,   1, &buffer[j*BYTES_STAR]);
      dparse (&stars[N].Y,   2, &buffer[j*BYTES_STAR]);
      dparse (&stars[N].mag, 3, &buffer[j*BYTES_STAR]);
      dparse (&type,         5, &buffer[j*BYTES_STAR]);
      if ((type == 4) || (type == 6) || (type == 8)) continue;
      N++;
    }
  }
  free (header.buffer);
  free (buffer);
 
  if (nstar != Nstars) {
    fprintf (stderr, "WARNING: only read %d of %d stars\n", nstar, Nstars);
  }
  if (N < 5) { 
    fprintf (stderr, "ERROR: too few stars for reliable solution, only %d\n",
	     N);
    exit (0);
  }

  sort_stars (stars, N);  /* sorting by magnitude */
  Nstars = N;

  if (VERBOSE) fprintf (stderr, "\nread %d stars from data file", Nstars);
  if (*NSTARS < Nstars) {
    REALLOCATE (stars, SStars, *NSTARS);
    if (VERBOSE) fprintf (stderr, ", using %d\n", *NSTARS);
  } else {
    *NSTARS = Nstars;
    if (VERBOSE) fprintf (stderr, "\n");
  }

  *dNdM = *NSTARS / (stars[(*NSTARS-1)].mag - stars[0].mag) ;
  if (VERBOSE) fprintf (stderr, "brightest star in datafile: %f mag\n", stars[0].mag);
  
  return (stars);
}

sort_stars (stars, N) 
SStars *stars;
int N;
{
  int l,j,ir,i;
  int temp;
  SStars tempstar;
  
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      tempstar = stars[--l];
    }
    else {
      tempstar = stars[ir];
      stars[ir] = stars[0];
      if (--ir == 0) {
	stars[0] = tempstar;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && stars[j].mag < stars[j+1].mag) ++j;
      if (tempstar.mag < stars[j].mag) {
	stars[i] = stars[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    stars[i] = tempstar;
  }
}


# if (0)

  /***  this is tailored for LONEOS ***/
  Nccd = -1;
  fits_scan (&header, "NCCD", "%d", 1, &Nccd);
  if (Nccd == -1) {  /* no ccd info in header, not loneos */
    coords[0].crpix1 = 0.5*(*NX);
    coords[0].crpix2 = 0.5*(*NY);
  }
  if (Nccd == 0) {  /* chip 0 (a) *** might be wrong *** */
    coords[0].crpix1 = 0;
    coords[0].crpix2 = 0.5*(*NY);
  }
  if (Nccd == 1) {  /* chip 1 (b) */
    coords[0].crpix1 = (*NX);
    coords[0].crpix2 = 0.5*(*NY);
  }

# endif
