# include "astro.h"

main (argc, argv)
int argc;
char **argv;
{

  int status;
  Header header;
  char  filename[1024], reffile[1024];
  double X_O, X_X, X_Y, Y_O, Y_X, Y_Y;
  double RA_O, RA_X, RA_Y, DEC_O, DEC_X, DEC_Y;
  double ra_o, ra_x, ra_y, dec_o, dec_x, dec_y;
  double dRA, dDEC, dX, dY, dra, ddec;

  while (fscanf (stdin, "%s", filename) != EOF) {

    status = fits_read_header (filename, &header);
    if (!status) {
      fprintf (stderr, "error opening file %s", filename);
      continue;
    }
    status = fits_scan (&header, "rREF", "%s", 1, reffile);
    status &= fits_scan (&header, "X_O", "%lf", 1, &X_O);
    status &= fits_scan (&header, "X_X", "%lf", 1, &X_X);
    status &= fits_scan (&header, "X_Y", "%lf", 1, &X_Y);
    status &= fits_scan (&header, "Y_O", "%lf", 1, &Y_O);
    status &= fits_scan (&header, "Y_X", "%lf", 1, &Y_X);
    status &= fits_scan (&header, "Y_Y", "%lf", 1, &Y_Y);
    status &= fits_scan (&header, "dX", "%lf", 1, &dX);
    status &= fits_scan (&header, "dY", "%lf", 1, &dY);
    if (!status) {
      fprintf (stderr, "file missing rastro info: %s", filename);
      fits_free_header (&header);
      continue;
    }

    status  = fits_read_header (reffilename, &refhead);
    if (!status) {
      fprintf (stderr, "error opening file %s", reffilename);
      continue;
    }
    status &= fits_scan (&refhead, "RA_O", "%lf", 1, &RA_O);
    status &= fits_scan (&refhead, "RA_X", "%lf", 1, &RA_X);
    status &= fits_scan (&refhead, "RA_Y", "%lf", 1, &RA_Y);
    status &= fits_scan (&refhead, "DEC_O", "%lf", 1, &DEC_O);
    status &= fits_scan (&refhead, "DEC_X", "%lf", 1, &DEC_X);
    status &= fits_scan (&refhead, "DEC_Y", "%lf", 1, &DEC_Y);
    status &= fits_scan (&refhead, "dRA",   "%lf", 1, &dRA);
    status &= fits_scan (&refhead, "dDEC",  "%lf", 1, &dDEC);
    if (!status) {
      fprintf (stderr, "file missing rastro info: %s", reffile);
      fits_free_header (&header);
      fits_free_header (&refhead);
      continue;
    }

    
    ra_o = RA_X*X_O + RA_Y*Y_O + RA_O;
    dec_o = DEC_X*X_O + DEC_Y*Y_O + DEC_O;
    
    ra_x = RA_X*X_X + RA_Y*Y_X;
    ra_y = RA_X*X_Y + RA_Y*Y_Y;
    dec_x = DEC_X*X_X + DEC_Y*Y_X;
    dec_y = DEC_X*X_Y + DEC_Y*Y_Y;
    
    dra = sqrt(dRA*dRA + SQ(3600*dX*RA_X) + SQ(3600*dY*RA_Y));
    ddec = sqrt(dDEC*dDEC + SQ(3600*dX*DEC_X) + SQ(3600*dY*DEC_Y));
			   
    status  = fits_modify (&header, "rREF", "%s", 1, reffile);
    status &= fits_modify (&header, "RA_O", "%lf", 1, RA_O);
    status &= fits_modify (&header, "RA_X", "%lf", 1, RA_X);
    status &= fits_modify (&header, "RA_Y", "%lf", 1, RA_Y);
    status &= fits_modify (&header, "DEC_O", "%lf", 1, DEC_O);
    status &= fits_modify (&header, "DEC_X", "%lf", 1, DEC_X);
    status &= fits_modify (&header, "DEC_Y", "%lf", 1, DEC_Y);
    status &= fits_modify (&header, "dRA",   "%lf", 1, dRA);
    status &= fits_modify (&header, "dDEC",  "%lf", 1, dDEC);
    
    fits_free_header (&header);
    fits_free_header (&refhead);
    
  }

}

