# include "astro.h"
# define NBYTES 160000

void extract_catalog (filename, ra, dec, mag, Nstars)
char *filename;
double **ra;
double **dec;
double **mag;
int     *Nstars;
{
  
  Header  header;
  Table   table;
  char **temp1;
  char origin[300];
  int i;
  double *temp2, R, D, minM, maxM;
  double minR, maxR, minD, maxD;

  if (*Nstars == 0) {
    ALLOCATE (*ra, double, 1);
    ALLOCATE (*dec, double, 1);
    ALLOCATE (*mag, double, 1);
  }
  
  fits_read_header (filename, &header);
  fits_scan (&header, "ORIGIN", "%s", 1, origin);
  fits_free_header (&header);

  minM = 100;
  maxM = 0;
  minR = 390;
  maxR = -10;
  minD = 100;
  maxD = -100;
  if (!strcmp (origin, "ST ScI")) {   /* Assuming it is guide star catalog.. */
    fits_read_table   (filename, &table);  
    fits_table_column (&table, "RA_DEG", "%lf", &temp2);
    REALLOCATE (*ra, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*ra)[*Nstars+i] = temp2[i];
      minR = MIN (temp2[i], minR);
      maxR = MAX (temp2[i], maxR);
    }
    free (temp2);
    fits_table_column (&table, "DEC_DEG", "%lf", &temp2);
    REALLOCATE (*dec, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*dec)[*Nstars+i] = temp2[i];
      minD = MIN (temp2[i], minD);
      maxD = MAX (temp2[i], maxD);
    }
    free (temp2);
    fits_table_column (&table, "MAG", "%lf", &temp2);
    REALLOCATE (*mag, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*mag)[*Nstars+i] = temp2[i];
      minM = MIN (temp2[i], minM);
      maxM = MAX (temp2[i], maxM);
    }
    *Nstars += table.Naxis[1];
    free (temp2);
    fprintf (stderr, "mag range: %f -- %f\n", minM, maxM);
    fprintf (stderr, "RA  range: %f -- %f, \n", minR, maxR);
    fprintf (stderr, "DEC range: %f -- %f, \n", minD, maxD);
    return;
  }

  if (!strncmp (origin, "U.Minn. PDP11/60", strlen("U.Minn. PDP11/60"))) {  /* Assuming it is Berkhuijsen etal catalog.. */
    fits_read_table   (filename, &table);  
    fits_table_column (&table, "RA(1950)", "%s", &temp1);
    REALLOCATE (*ra, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*ra)[*Nstars+i] = 15*dms_to_deg (temp1[i], ':');
      minR = MIN ((*ra)[*Nstars+i], minR);
      maxR = MAX ((*ra)[*Nstars+i], maxR);
      free (temp1[i]);
    }
    free (temp1);
    fits_table_column (&table, "DEC(1950)", "%s", &temp1);
    REALLOCATE (*dec, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*dec)[*Nstars+i] = dms_to_deg (temp1[i], ':');
      minD = MIN ((*dec)[*Nstars+i], minD);
      maxD = MAX ((*dec)[*Nstars+i], maxD);
      free (temp1[i]);
    }
    free (temp1);
    fits_table_column (&table, "V-MAG", "%lf", &temp2);
    REALLOCATE (*mag, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*mag)[*Nstars+i] = temp2[i];
      minM = MIN (temp2[i], minM);
      maxM = MAX (temp2[i], maxM);
    }
    free (temp2);
    for (i = *Nstars; i < *Nstars + table.Naxis[1]; i++) {
      if (i < 20) 
	fprintf (stderr, "%lf %lf -> ", (*ra)[i], (*dec)[i]);
      /* see Magnier et al. 1993 A&A 272 695 */
      R = 0.99999848*(*ra)[i]  - 0.00152332*(*dec)[i] + 0.0614972;
      D = 0.99999848*(*dec)[i] + 0.0009958* (*ra)[i]  - 0.0093577;
      precess (&R, &D, 1949.999790, 2000.0);  /* B1950.0 -> J2000.0 */
      (*ra)[i] = R;
      (*dec)[i] = D;
      if (i < 20) 
	fprintf (stderr, "%lf %lf\n", (*ra)[i], (*dec)[i]);
    }
    *Nstars += table.Naxis[1];
    fprintf (stderr, "mag range: %f -- %f\n", minM, maxM);
    return;
  }

  if (!strncmp (origin, "MIT-CSR", strlen("MIT-CSR"))) {  /* Assuming it is Magnier etal catalog.. */
    fits_read_table   (filename, &table);  
    fits_table_column (&table, "RA", "%lf", &temp2);
    REALLOCATE (*ra, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*ra)[*Nstars+i] = temp2[i];
      minR = MIN ((*ra)[*Nstars+i], minR);
      maxR = MAX ((*ra)[*Nstars+i], maxR);
    }
    free (temp2);
    fits_table_column (&table, "DEC", "%lf", &temp2);
    REALLOCATE (*dec, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*dec)[*Nstars+i] = temp2[i];
      minD = MIN ((*dec)[*Nstars+i], minD);
      maxD = MAX ((*dec)[*Nstars+i], maxD);
    }
    free (temp2);
    fits_table_column (&table, "V-MAG", "%lf", &temp2);
    REALLOCATE (*mag, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*mag)[*Nstars+i] = temp2[i];
      minM = MIN (temp2[i], minM);
      maxM = MAX (temp2[i], maxM);
    }
    free (temp2);
    *Nstars += table.Naxis[1];
    fprintf (stderr, "mag range: %f -- %f\n", minM, maxM);
    return;
  }

  if (!strncmp (origin, "RELPHOT", strlen("RELPHOT"))) {  /* Assuming it is from relphot.c */
    fits_read_table   (filename, &table);  
    fits_table_column (&table, "RA", "%lf", &temp2);
    REALLOCATE (*ra, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*ra)[*Nstars+i] = temp2[i];
      minR = MIN ((*ra)[*Nstars+i], minR);
      maxR = MAX ((*ra)[*Nstars+i], maxR);
    }
    free (temp2);
    fits_table_column (&table, "DEC", "%lf", &temp2);
    REALLOCATE (*dec, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*dec)[*Nstars+i] = temp2[i];
      minD = MIN ((*dec)[*Nstars+i], minD);
      maxD = MAX ((*dec)[*Nstars+i], maxD);
    }
    free (temp2);
    fits_table_column (&table, "MAG", "%lf", &temp2);
    REALLOCATE (*mag, double, *Nstars + table.Naxis[1]);
    for (i = 0; i < table.Naxis[1]; i++) {
      (*mag)[*Nstars+i] = temp2[i];
      minM = MIN (temp2[i], minM);
      maxM = MAX (temp2[i], maxM);
    }
    free (temp2);
    *Nstars += table.Naxis[1];
    fprintf (stderr, "mag range: %f -- %f\n", minM, maxM);
    return;
  }

  fprintf (stderr, "unknown catalog format, ORIGIN: %s\n", origin);
  exit (0);

}
