/*  readphotom.c  */
/*  (c)  Jure Skvarc  May 2000  */
/* A program which reads stars from the phot catalog, which are  */
/* contained in image with center coordinates */
/*  (ra, dec), rotation angle fi  and dimensions (width, height)     */
/*  and outputs their pixel coordinates and celestial coordinates        */
/*  TAN projection is assumed                    */
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <fitsio.h>
#include <getopt.h>

#include "consts.h"
#include "structs.h"
#include "functs.h"
/*  Read stars from file name which have declination between */
/*  declow and dechigh  */
int
read_pht_region(char *name, char *accname, double declow, double dechigh, STAR **stars, int *n, char **error)

{
  FILE *accfp, *fp;
  PHT_ACCEL acc[19];
  int i, k, nacc;
  int lok = 0, hok = 0;
  int indlow = 0, indhigh;

  /*  Read the accelerator file  */
  accfp = fopen(accname, "r");
  if (accfp == NULL) {
    asprintf(error, "Can't open the PHT accelerator file %s!", accname);
    return FILE_ERROR;
  }
  nacc = 0;
  while (!feof(accfp)) {
    k = fscanf(accfp, "%f %d %lf", &acc[nacc].limit, &acc[nacc].num, &acc[nacc].reallimit);
    if (k == 3) nacc++;
  }
  fclose(accfp);
  indhigh = acc[nacc - 1].num;
  /*  Find the region which has to be read in  */
  for (i = 1; i < nacc && !(lok && hok); i++) {
    if (!lok && declow >= acc[i - 1].limit && declow < acc[i].limit) {
      indlow = acc[i - 1].num;
      lok = 1;
    }
    if (!hok && dechigh >= acc[i - 1].limit && dechigh < acc[i].limit) {
      indhigh = acc[i].num;
      hok = 1;
    }
  }
  *n = indhigh - indlow;
  /*  Open the catalog file  */
  fp = fopen(name, "r");
  if (fp == NULL) {
    asprintf(error, "Can't open the PHT catalog file %s!", name);
    return FILE_ERROR;
  }
  /*  Reserve space for stars  */
  *stars = (STAR *) malloc(sizeof(STAR) * *n);
  if (*stars == NULL) {
    asprintf(error, "Can't reserve space for stars!");
    return MEMORY_ERROR;
  }
  /*  Read all of the stars between the two indexes  */
  if (fseek(fp, indlow * sizeof(STAR), SEEK_SET)) {
    asprintf(error, "Error in fseek for file %s", name);
    return IO_ERROR;
  }
  k = fread(*stars, sizeof(STAR), *n, fp);
  if (k != *n) {
    asprintf(error, "Error in reading star data from file %s!", name);
    return IO_ERROR;
  }
  fclose(fp);
  return 0;
}

/*  transforms s list of stars in STAR (PHT) format into a list in W_STAR format */
/*  it copies only those which fall within image limits  */
/*    */
int
extract_pht_stars(char *name, char *accname, double declow, double dechigh, 
		  W_STAR *starlist, MAP_IMAGE *image, double jd)

{
  int i, nstar, err, n;
  double xpix, ypix, xinc, yinc;
  STAR *star;
  int wmin, hmin, wmax, hmax;
  double dt;
  char *error;

  /*Always zero, because proper motion information is not in STAR structure */
#if 0
  int proper_motion = 0; 
#endif  
  nstar = image->nm;
  xinc = image->width / image->w;
  yinc = image->height / image->h;
  
  wmin = -image->extarea * image->w;
  hmin = -image->extarea * image->h;
  wmax =  (1 + image->extarea) * image->w;
  hmax =  (1 + image->extarea) * image->h;
  if ((err = read_pht_region(name, accname, declow, dechigh, &star, &n, &error))) {
    return err;
  }
  for (i = 0; i < n; i++) {
    /*  Check if the star falls into the image  */
    err = 0;
    fits_world_to_pix(star[i].ra, star[i].dec, image->ra, image->dec, 
		0.5 * image->w, 0.5 * image->h, 
		xinc, yinc, image->fi, "-TAN", 
		&xpix, &ypix, &err);
    if (!err && xpix >= wmin && xpix < wmax && ypix >= hmin && ypix < hmax) {
      dt = (jd - 2451545) / 365.25;
      strcpy(starlist[nstar].name, star[i].name);
#if 0
      if (proper_motion) {
	/*  Make correction for the proper motion  */
	starlist[nstar].ra = star[i].ra + star[i].pmra * dt / 240.0;
	if (starlist[nstar].ra > 360) starlist[nstar].ra -= 360;
	if (starlist[nstar].ra < 0) starlist[nstar].ra += 360;
	/*  Make correction for the proper motion  */
	starlist[nstar].dec = star[i].dec + star[i].pmdec * dt / 3600.0;
      }
#endif
      starlist[nstar].ra = star[i].ra;
      starlist[nstar].dec = star[i].dec;
      starlist[nstar].sra = star[i].sra;
      starlist[nstar].sdec = star[i].sdec;
      starlist[nstar].u = star[i].u;
      starlist[nstar].su = star[i].su;
      starlist[nstar].b = star[i].b;
      starlist[nstar].sb = star[i].sb;
      starlist[nstar].v = star[i].v;
      starlist[nstar].sv = star[i].sv;
      starlist[nstar].r = star[i].r;
      starlist[nstar].sr = star[i].sr;
      starlist[nstar].i = star[i].i;
      starlist[nstar].si = star[i].si;
      starlist[nstar].x = xpix - 1;
      starlist[nstar].y = ypix - 1;
      nstar++;
      if (nstar == MAXSTARNUM) {
	image->nm = nstar;
	return TOO_MANY_STARS;
      }
    }
  }
  free(star);
  image->nm = nstar;
  return 0;
}


int
get_pht_stars(MAP_IMAGE *image, W_STAR *starlist, char *path, double jd) 

{
  STAR corner[4];
  double x, w;
  double declow, dechigh;
  double ralow[2], rahigh[2];
  int rarange;
  double xpix, ypix;
  int i, j;
  int err;
  char *fname, *accname;
  int indlow, indhigh;
  double xinc = image->width / image->w;
  double yinc = image->height / image->h;

  /*  Reserve space for the file name  */
  if ((fname = malloc(sizeof(char) * (strlen(path) + 50))) == NULL) {
    return MEMORY_ERROR;
  }
  if ((accname = malloc(sizeof(char) * (strlen(path) + 50))) == NULL) {
    return MEMORY_ERROR;
  }
  if ((err = calculate_corners(image, corner))) return err;
  /*  Find the highest and the lowest declination value  */
  declow = dechigh = corner[0].dec;
  ralow[0] = rahigh[0] = corner[0].ra;
  for (i = 1; i < 4; i++) {
    if (corner[i].dec < declow) declow = corner[i].dec;
    if (corner[i].dec > dechigh) dechigh = corner[i].dec;
    if (corner[i].ra < ralow[0]) ralow[0] = corner[i].ra;
    if (corner[i].ra > rahigh[0]) rahigh[0] = corner[i].ra;
  }
  w = max(image->width, image->height);
  x = declow - image->extarea * w;
  if (x > -90.0) declow = x;
  else declow = -90;
  x = dechigh + image->extarea * w;
  if (x < 90) dechigh = x;
  else dechigh = 90;

  x = ralow[0] - image->extarea * w;
  if (x < 0) {
    rahigh[0] = x + 360.0;
    ralow[0] = rahigh[0] + image->extarea * w;
  }
  else {
    ralow[0] = x;
    x = rahigh[0] + image->extarea * w;
    if (x > 360) {
      rahigh[0] = ralow[0];
      ralow[0] = x - 360.0;
    }
  }

  /* Check if the north pole falls inside the rectangle  */
  err = 0;
  fits_world_to_pix(0.0, 90.0, image->ra, image->dec, 
		    0.5 * image->w, 0.5 * image->h, 
	      xinc, yinc, image->fi, "-TAN", 
	      &xpix, &ypix, &err);
  if (!err && (xpix >= 0.0 && xpix <= image->w && ypix >= 0.0 && ypix <= image->h)) {
      ralow[0] = 0;
      rahigh[0] = 360;
      rarange = 1;
      dechigh = 90.0;
  }
  else {
    /* Check if the south pole falls inside the rectangle  */
    err = 0;
    fits_world_to_pix(0.0, -90.0, image->ra, image->dec,
		      0.5 * image->w, 0.5 * image->h, 
		      xinc, yinc, image->fi, "-TAN", 
		      &xpix, &ypix, &err);
    if (!err && (xpix >= 0.0 && xpix <= image->w && ypix >= 0.0 && ypix <= image->h)) {
      ralow[0] = 0;
      rahigh[0] = 360;
      rarange = 1;
      declow = -90.0;
    }
    else {
      if (rahigh[0] - ralow[0] >= 180) {
	ralow[1] = rahigh[0];
	rahigh[1] = 360;
	rahigh[0] = ralow[0];
	ralow[0] = 0;
	rarange = 2;
      }
      else {
	rarange = 1;
      }
    }
  }
  for (i = 0; i < rarange; i++) {
    /* indexes  of files which contain stars of interest  */
    indlow = floor(ralow[i] / PHT_RAWIDTH);
    indhigh = floor((rahigh[i] - 1e-9) / PHT_RAWIDTH);
    /*  Visit all files which contain stars  */
    for (j = indlow; j <= indhigh; j++) {
      /*  Construct file names  */
      sprintf(fname, "%s/phot%02d%02d.bin", path, j / 2, (j % 2) * 30);
      sprintf(accname, "%s/phot%02d%02d.acc", path, j / 2, (j % 2) * 30);
      err = extract_pht_stars(fname, accname, declow, dechigh, starlist, image, jd);
      if (err) return err;
    }
  }
  free(accname);
  free(fname);
  return 0;
}