#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <wcslib/cel.h>
#include <cpgplot.h>
#include "qfits.h"
#include <gsl/gsl_multifit.h>
#include <getopt.h>

#define NMAX 8192
#define LIM 128
#define D2R M_PI/180.0
#define R2D 180.0/M_PI

struct image {
  int naxis,naxis1,naxis2,nframes;
  float *zavg,*zstd,*zmax,*znum;
  double ra0,de0;
  float x0,y0;
  float a[2],b[2];
  double mjd;
  float *dt;
};
struct transformation {
  double mjd;
  double ra0,de0;
  float a[3],b[3];
  float x0,y0;
  float xrms,yrms,rms;
};
struct star {
  double ra,de;
  float pmra,pmde;
  float mag;
};
struct catalog {
  int n;
  float x[NMAX],y[NMAX],imag[NMAX],fm[NMAX],fb[NMAX],bg[NMAX];
  double ra[NMAX],de[NMAX],vmag[NMAX];
  double rx[NMAX],ry[NMAX];
  float xres[NMAX],yres[NMAX],res[NMAX];
  int usage[NMAX];
};
struct image read_fits(char *filename);
void forward(double ra0,double de0,double ra,double de,double *x,double *y);
void reverse(double ra0,double de0,double x,double y,double *ra,double *de);
double gmst(double mjd);
double modulo(double x,double y);
int fgetline(FILE *file,char *s,int lim);
struct catalog match_catalogs(char *pixcat,char *astcat,struct transformation t,struct image img,float rmax,float mmin);
void plot_astrometric_catalog(struct transformation t,struct image img,float mmin);
void plot_pixel_catalog(char *filename);
void lfit2d(float *x,float *y,float *z,int n,float *a);
void add_fits_keywords(struct transformation t,char *filename);
void modify_fits_keywords(struct transformation t,char *filename);
void precess(double mjd0,double ra0,double de0,double mjd,double *ra,double *de);

void plot_image(struct image img,struct transformation t,struct catalog c,char *filename,float mmin)
{
  int i;
  float tr[]={-0.5,1.0,0.0,-0.5,0.0,1.0};
  float heat_l[]={0.0,0.2,0.4,0.6,1.0};
  float heat_r[]={0.0,0.5,1.0,1.0,1.0};
  float heat_g[]={0.0,0.0,0.5,1.0,1.0};
  float heat_b[]={0.0,0.0,0.0,0.3,1.0};
  float zmin,zmax,zavg,zstd;

  for (i=0,zavg=0.0;i<img.naxis1*img.naxis2;i++)
    zavg+=img.zavg[i];
  zavg/=(float) img.naxis1*img.naxis2;
  for (i=0,zstd=0.0;i<img.naxis1*img.naxis2;i++)
    zstd+=pow(img.zavg[i]-zavg,2);
  zstd=sqrt(zstd/(float) (img.naxis1*img.naxis2));
  zmin=zavg-2*zstd;
  zmax=zavg+6*zstd;

  cpgopen("1/xs");
  cpgwnad(0.0,img.naxis1,0.0,img.naxis2);
  cpgctab (heat_l,heat_r,heat_g,heat_b,5,1.0,0.5);
    
  cpgimag(img.zavg,img.naxis1,img.naxis2,1,img.naxis1,1,img.naxis2,zmin,zmax,tr);
  cpgbox("BCTSNI",0.,0,"BCTSNI",0.,0);

  cpgsci(3);
  plot_pixel_catalog(filename);

  cpgsci(4);
  plot_astrometric_catalog(t,img,mmin);
  cpgsci(2);
  for (i=0;i<c.n;i++)
    cpgpt1(c.x[i]+t.x0,c.y[i]+t.y0,24);

  cpgend();

  return;
}

void usage(float mmin,float rmin)
{
  printf("addwcs: Add/fit World Coordinate System to a FITS file\n\n");

  printf("-f <file>:    FITS file to add/fit WCS to [required]\n");
  printf("-r <file>:    FITS file with reference WCS [required]\n");
  printf("-m <float>:   Magnitude cut-off in Tycho-2 catalog [optional; default %.1f]\n",mmin);
  printf("-R <float>:   Radius cut-off for matching [optional; default %.1f pix]\n",rmin);
  printf("-p            Plot image and selected stars [optional]\n");
  printf("-a            Add WCS keywords to input file (instead of modify) [optional]\n");
  printf("-t            Track on a fixed RA/Dec (correct for field rotation)\n");
  printf("-h            Print this help\n");

  return;
}

// Get reference transformation
struct transformation reference(char *filename)
{
  struct transformation t;

  t.mjd=atof(qfits_query_hdr(filename,"MJD-OBS"));
  t.ra0=atof(qfits_query_hdr(filename,"CRVAL1"));
  t.de0=atof(qfits_query_hdr(filename,"CRVAL2"));
  t.x0=atof(qfits_query_hdr(filename,"CRPIX1"));
  t.y0=atof(qfits_query_hdr(filename,"CRPIX2"));
  t.a[0]=0.0;
  t.a[1]=3600.0*atof(qfits_query_hdr(filename,"CD1_1"));
  t.a[2]=3600.0*atof(qfits_query_hdr(filename,"CD1_2"));
  t.b[0]=0.0;
  t.b[1]=3600.0*atof(qfits_query_hdr(filename,"CD2_1"));
  t.b[2]=3600.0*atof(qfits_query_hdr(filename,"CD2_2"));

  return t;
}

void rotate(float theta,float *x,float *y)
{
  float ct,st;
  float x0,y0;
  
  ct=cos(theta*D2R);
  st=sin(theta*D2R);
  x0= *x;
  y0= *y;

  *x=ct*x0-st*y0;
  *y=st*x0+ct*y0;

  return;
}

int main(int argc,char *argv[])
{
  int i,j,k,l,m;
  struct transformation t;
  struct image img;
  char *fitsfile=NULL,*reffile=NULL,catfile[128],calfile[128];
  FILE *outfile;
  struct catalog c;
  float mmin=10.0,rmin=10.0;
  double mjd0=51544.5,ra0,de0,ra1,de1;
  float q0,q1;
  float rmsmin;
  float x[NMAX],y[NMAX],rx[NMAX],ry[NMAX];
  int arg=0,plot=0,add=0,track=0;
  char *env,starfile[128];

  // Environment variables
  env=getenv("ST_DATADIR");
  sprintf(starfile,"%s/data/tycho2.dat",env);

  // Decode options
  if (argc>1) {
    while ((arg=getopt(argc,argv,"f:r:m:R:hpnta"))!=-1) {
      switch (arg) {
	
      case 'f':
	fitsfile=optarg;
	break;
	
      case 'r':
	reffile=optarg;
	break;
	
      case 'm':
	mmin=atof(optarg);
	break;

      case 't':
	track=1;
	break;

      case 'R':
	rmin=atof(optarg);
	break;
	
      case 'p':
	plot=1;
	break;

      case 'a':
	add=1;
	break;
	
      case 'h':
	usage(mmin,rmin);
	return 0;
	
      default:
	usage(mmin,rmin);
	return 0;
      }
    } 
  } else {
    usage(mmin,rmin);
    return 0;
  }

  // Check if minimum input is provided
  if (fitsfile==NULL || reffile==NULL) {
    usage(mmin,rmin);
    return 0;
  }

  // Check this is indeed a FITS file 
  if (is_fits_file(fitsfile)!=1) {
    printf("%s is not a FITS file\n",fitsfile);
    return -1 ;
  }
    
  // Check this is indeed a FITS file 
  if (is_fits_file(reffile)!=1) {
    printf("%s is not a FITS file\n",reffile);
    return -1 ;
  }  

  // Read fits file
  img=read_fits(fitsfile);
  sprintf(catfile,"%s.cat",fitsfile);
  sprintf(calfile,"%s.cal",fitsfile);

  // Read reference transformation
  t=reference(reffile);

  // Correct astrometry for fixed or tracked setup
  if (track==0) {
    precess(mjd0,t.ra0,t.de0,t.mjd,&ra1,&de1);
    ra1=modulo(ra1+gmst(img.mjd)-gmst(t.mjd),360.0);
    precess(img.mjd,ra1,de1,mjd0,&t.ra0,&t.de0);
  }

  // Match catalog
  c=match_catalogs(catfile,starfile,t,img,rmin,mmin);

  
  // Plot
  if (plot==1)
    plot_image(img,t,c,catfile,mmin);

  // Do fit
  if (c.n>10) {
    for (l=0;l<10;l++) {
      for (j=0;j<5;j++) {
	// Transform
	for (i=0;i<c.n;i++) 
	  forward(t.ra0,t.de0,c.ra[i],c.de[i],&c.rx[i],&c.ry[i]);
      
	// Select
	for (i=0,k=0;i<c.n;i++) {
	  if (c.usage[i]==1) {
	    x[k]=c.x[i];
	    y[k]=c.y[i];
	    rx[k]=(float) c.rx[i];
	    ry[k]=(float) c.ry[i];
	    k++;
	  }
	}
	// Fit
	lfit2d(x,y,rx,k,t.a);
	lfit2d(x,y,ry,k,t.b);
	
	// Move reference point
	reverse(t.ra0,t.de0,t.a[0],t.b[0],&ra0,&de0);
	t.ra0=ra0;
	t.de0=de0;
      }
      
      // Compute and plot residuals
      for (i=0,t.xrms=0.0,t.yrms=0.0,m=0;i<c.n;i++) {
	if (c.usage[i]==1) {
	  c.xres[i]=c.rx[i]-(t.a[0]+t.a[1]*c.x[i]+t.a[2]*c.y[i]);
	  c.yres[i]=c.ry[i]-(t.b[0]+t.b[1]*c.x[i]+t.b[2]*c.y[i]);
	  
	  c.res[i]=sqrt(c.xres[i]*c.xres[i]+c.yres[i]*c.yres[i]);
	  t.xrms+=c.xres[i]*c.xres[i];
	  t.yrms+=c.yres[i]*c.yres[i];
	  t.rms+=c.xres[i]*c.xres[i]+c.yres[i]*c.yres[i];
	  m++;
	}
      }
      t.xrms=sqrt(t.xrms/(float) m);
      t.yrms=sqrt(t.yrms/(float) m);
      t.rms=sqrt(t.rms/(float) m);
      
      // Deselect outliers
      for (i=0;i<c.n;i++) {
	if (c.res[i]>2*t.rms)
	  c.usage[i]=0;
      }
    }
  } else {
    t.xrms=0.0;
    t.yrms=0.0;
    t.rms=0.0;
  }

  // Print results
  outfile=fopen(calfile,"w");
  for (i=0;i<c.n;i++) 
    if (c.usage[i]==1)
      fprintf(outfile,"%10.4f %10.4f %10.6f %10.6f %8.3f %8.3f %8.3f %8.3f %8.3f\n",c.x[i],c.y[i],c.ra[i],c.de[i],c.vmag[i],c.imag[i],c.fb[i],c.fm[i],c.bg[i]);
  fclose(outfile);

  printf("%s %8.4lf %8.4lf ",fitsfile,t.ra0,t.de0);
  printf("%3d/%3d %6.1f %6.1f %6.1f\n",m,c.n,t.xrms,t.yrms,t.rms);

  // Add keywords
  if (add==1)
    add_fits_keywords(t,fitsfile);
  else
    modify_fits_keywords(t,fitsfile);

  return 0;
}

// Read fits image
struct image read_fits(char *filename)
{
  int i,j,k,l,m;
  qfitsloader ql;
  char key[FITS_LINESZ+1];
  char val[FITS_LINESZ+1];
  struct image img;

  // Image size
  img.naxis=atoi(qfits_query_hdr(filename,"NAXIS"));
  img.naxis1=atoi(qfits_query_hdr(filename,"NAXIS1"));
  img.naxis2=atoi(qfits_query_hdr(filename,"NAXIS2"));

  // MJD
  img.mjd=(double) atof(qfits_query_hdr(filename,"MJD-OBS"));

  // Set parameters
  ql.xtnum=0;
  ql.ptype=PTYPE_FLOAT;
  ql.filename=filename;

  if (img.naxis==3) {
    // Number of frames
    img.nframes=atoi(qfits_query_hdr(filename,"NFRAMES"));
    
    // Timestamps
    img.dt=(float *) malloc(sizeof(float)*img.nframes);
    for (i=0;i<img.nframes;i++) {
      sprintf(key,"DT%04d",i);
      strcpy(val,qfits_query_hdr(filename,key));
      sscanf(val+1,"%f",&img.dt[i]);
      //    img.dt[i]=atof(qfits_query_hdr(filename,key));
    }
    
    // Allocate image memory
    img.zavg=(float *) malloc(sizeof(float)*img.naxis1*img.naxis2);
    img.zstd=(float *) malloc(sizeof(float)*img.naxis1*img.naxis2);
    img.zmax=(float *) malloc(sizeof(float)*img.naxis1*img.naxis2);
    img.znum=(float *) malloc(sizeof(float)*img.naxis1*img.naxis2);
    
    
    // Loop over planes
    for (k=0;k<4;k++) {
      ql.pnum=k;;
      
      // Initialize load
      if (qfitsloader_init(&ql) != 0) 
	printf("Error initializing data loading\n");
      
      // Test load
      if (qfits_loadpix(&ql) != 0) 
	printf("Error loading actual data\n");
      
      // Fill z array
      for (i=0,l=0;i<img.naxis1;i++) {
	for (j=0;j<img.naxis2;j++) {
	  if (k==0) img.zavg[l]=ql.fbuf[l];
	  if (k==1) img.zstd[l]=ql.fbuf[l];
	  if (k==2) img.zmax[l]=ql.fbuf[l];
	  if (k==3) img.znum[l]=ql.fbuf[l];
	  l++;
	}
      }
    }
  } else {
    // Allocate image memory
    img.zavg=(float *) malloc(sizeof(float)*img.naxis1*img.naxis2);
    
    ql.pnum=0;
    
    // Initialize load
    if (qfitsloader_init(&ql) != 0) 
      printf("Error initializing data loading\n");
    
    // Test load
    if (qfits_loadpix(&ql) != 0) 
      printf("Error loading actual data\n");
      
    // Fill z array
    for (i=0,l=0;i<img.naxis1;i++) {
      for (j=0;j<img.naxis2;j++) {
	img.zavg[l]=ql.fbuf[l];
	l++;
      }
    }
  }

  return img;
}


// Greenwich Mean Sidereal Time
double gmst(double mjd)
{
  double t,gmst;

  t=(mjd-51544.5)/36525.0;

  gmst=modulo(280.46061837+360.98564736629*(mjd-51544.5)+t*t*(0.000387933-t/38710000),360.0);

  return gmst;
}

// Return x modulo y [0,y)
double modulo(double x,double y)
{
  x=fmod(x,y);
  if (x<0.0) x+=y;

  return x;
}

// Read a line of maximum length int lim from file FILE into string s
int fgetline(FILE *file,char *s,int lim)
{
  int c,i=0;
 
  while (--lim > 0 && (c=fgetc(file)) != EOF && c != '\n')
    s[i++] = c;
  if (c == '\n')
    s[i++] = c;
  s[i] = '\0';
  return i;
}

// Match catalogs
struct catalog match_catalogs(char *pixcat,char *astcat,struct transformation t,struct image img,float rmax,float mmin)
{
  int i=0,imin,j,k,np;
  FILE *file;
  char line[LIM];
  struct star s;
  double rx,ry,d,dx,dy;
  int usage[NMAX];
  float xp[NMAX],yp[NMAX],mp[NMAX],x,y,fb[NMAX],fm[NMAX],bg[NMAX];
  struct catalog c;
  float r,rmin;

  // Read pixel catalog
  file=fopen(pixcat,"r");
  if (file==NULL) {
    printf("pixel catalog not found\n");
    exit(-1);
  }
  while (fgetline(file,line,LIM)>0) {
    if (strstr(line,"#")!=NULL) 
      continue;
    sscanf(line,"%f %f %f %f %f %f",&xp[i],&yp[i],&mp[i],&fb[i],&fm[i],&bg[i]);
    usage[i]=1;
    i++;
  }
  fclose(file);
  np=i;

  // Denominator
  d=t.a[1]*t.b[2]-t.a[2]*t.b[1];

  // Read astrometric catalog
  file=fopen(astcat,"rb");
  if (file==NULL) {
    printf("astrometric catalog not found\n");
    exit(-1);
  }
  j=0;
  while (!feof(file)) {
    fread(&s,sizeof(struct star),1,file);
    if (s.mag>mmin)
      continue;
    r=acos(sin(t.de0*D2R)*sin(s.de*D2R)+cos(t.de0*D2R)*cos(s.de*D2R)*cos((t.ra0-s.ra)*D2R))*R2D;
    if (r>90.0)
      continue;
    forward(t.ra0,t.de0,s.ra,s.de,&rx,&ry);

    dx=rx-t.a[0];
    dy=ry-t.b[0];
    x=(t.b[2]*dx-t.a[2]*dy)/d+t.x0;
    y=(t.a[1]*dy-t.b[1]*dx)/d+t.y0;

    // On image
    if (x>0.0 && x<img.naxis1 && y>0.0 && y<img.naxis2) {
      // Loop over pixel catalog
      for (i=0;i<np;i++) {
	r=sqrt(pow(xp[i]-x,2)+pow(yp[i]-y,2));
	if (i==0 || r<rmin) {
	  rmin=r;
	  imin=i;
	}
      }

      // Select
      if (rmin<rmax && usage[imin]==1) {
	c.x[j]=xp[imin]-t.x0;
	c.y[j]=yp[imin]-t.y0;
	c.imag[j]=mp[imin];
	c.fb[j]=fb[imin];
	c.fm[j]=fm[imin];
	c.bg[j]=bg[imin];
	c.ra[j]=s.ra;
	c.de[j]=s.de;
	c.vmag[j]=s.mag;
	c.usage[j]=1;
	usage[imin]=0;
	j++;
      }
    }
  }
  fclose(file);
  c.n=j;

  return c;
}

// Plot astrometric catalog
void plot_astrometric_catalog(struct transformation t,struct image img,float mmin)
{
  int i=0;
  FILE *file;
  struct star s;
  double rx,ry,d,r;
  double ra,de;
  float x,y;
  char *env,starfile[128];

  // Environment variables
  env=getenv("ST_DATADIR");
  sprintf(starfile,"%s/data/tycho2.dat",env);

  d=t.a[1]*t.b[2]-t.a[2]*t.b[1];

  file=fopen(starfile,"rb");
  while (!feof(file)) {
    fread(&s,sizeof(struct star),1,file);
    if (s.mag>mmin)
      continue;
    r=acos(sin(t.de0*D2R)*sin(s.de*D2R)+cos(t.de0*D2R)*cos(s.de*D2R)*cos((t.ra0-s.ra)*D2R))*R2D;
    if (r>90.0)
      continue;
    forward(t.ra0,t.de0,s.ra,s.de,&rx,&ry);
    x=(t.b[2]*rx-t.a[2]*ry)/d+t.x0;
    y=(t.a[1]*ry-t.b[1]*rx)/d+t.y0;
    if (x>0.0 && x<img.naxis1 && y>0.0 && y<img.naxis2)
      cpgpt1(x,y,24);
  }
  fclose(file);

  return;
}

// Plot pixel catalog
void plot_pixel_catalog(char *filename)
{
  int i=0;
  FILE *file;
  char line[LIM];
  float x,y,mag;

  // Read catalog
  file=fopen(filename,"r");
  while (fgetline(file,line,LIM)>0) {
    if (strstr(line,"#")!=NULL) 
      continue;
    sscanf(line,"%f %f %f",&x,&y,&mag);

    cpgpt1(x,y,4);

    i++;
  }
  fclose(file);

  return;
}

// Linear 2D fit
void lfit2d(float *x,float *y,float *z,int n,float *a)
{
  int i,j,m;
  double chisq;
  gsl_matrix *X,*cov;
  gsl_vector *yy,*w,*c;

  X=gsl_matrix_alloc(n,3);
  yy=gsl_vector_alloc(n);
  w=gsl_vector_alloc(n);

  c=gsl_vector_alloc(3);
  cov=gsl_matrix_alloc(3,3);

  // Fill matrices
  for(i=0;i<n;i++) {
    gsl_matrix_set(X,i,0,1.0);
    gsl_matrix_set(X,i,1,x[i]);
    gsl_matrix_set(X,i,2,y[i]);
    
    gsl_vector_set(yy,i,z[i]);
    gsl_vector_set(w,i,1.0);
  }

  // Do fit
  gsl_multifit_linear_workspace *work=gsl_multifit_linear_alloc(n,3);
  gsl_multifit_wlinear(X,w,yy,c,cov,&chisq,work);
  gsl_multifit_linear_free(work);

  // Save parameters
  for (i=0;i<3;i++)
    a[i]=gsl_vector_get(c,(i));

  gsl_matrix_free(X);
  gsl_vector_free(yy);
  gsl_vector_free(w);
  gsl_vector_free(c);
  gsl_matrix_free(cov);

  return;
}

// Add FITS keywords
void add_fits_keywords(struct transformation t,char *filename) 
{
  int i,j,k,l,m;
  int naxis1,naxis2,naxis3;
  qfits_header *qh;
  qfitsdumper qd;
  qfitsloader ql;
  char key[FITS_LINESZ+1];
  char val[FITS_LINESZ+1];
  char com[FITS_LINESZ+1];
  char lin[FITS_LINESZ+1];
  FILE *file;
  float *fbuf;

  naxis1=atoi(qfits_query_hdr(filename,"NAXIS1"));
  naxis2=atoi(qfits_query_hdr(filename,"NAXIS2"));
  naxis3=atoi(qfits_query_hdr(filename,"NAXIS3"));

  fbuf=malloc(sizeof(float)*naxis1*naxis2*naxis3);

  // Read header
  qh=qfits_header_read(filename);

  ql.xtnum=0;
  ql.ptype=PTYPE_FLOAT;
  ql.filename=filename;
  for (k=0,l=0;k<naxis3;k++) {
    ql.pnum=k;
    // Initialize load
    if (qfitsloader_init(&ql) != 0) 
      printf("Error initializing data loading\n");

    // Test load
    if (qfits_loadpix(&ql) != 0) 
      printf("Error loading actual data\n");

    for (i=0,m=0;i<naxis1;i++) {
      for (j=0;j<naxis2;j++) {
	fbuf[l]=ql.fbuf[m];
	l++;
	m++;
      }
    }
  }

  sprintf(val,"%e",t.yrms/3600.0);
  qfits_header_add_after(qh,"MJD-OBS","CRRES2",val," ",NULL);
  sprintf(val,"%e",t.xrms/3600.0);
  qfits_header_add_after(qh,"MJD-OBS","CRRES1",val," ",NULL);
  qfits_header_add_after(qh,"MJD-OBS","CUNIT2","'deg'"," ",NULL);
  qfits_header_add_after(qh,"MJD-OBS","CUNIT1","'deg'"," ",NULL);
  qfits_header_add_after(qh,"MJD-OBS","CTYPE2","'DEC--TAN'"," ",NULL);
  qfits_header_add_after(qh,"MJD-OBS","CTYPE1","'RA---TAN'"," ",NULL);
  sprintf(val,"%e",t.b[2]/3600.0);
  qfits_header_add_after(qh,"MJD-OBS","CD2_2",val," ",NULL);
  sprintf(val,"%e",t.b[1]/3600.0);
  qfits_header_add_after(qh,"MJD-OBS","CD2_1",val," ",NULL);
  sprintf(val,"%e",t.a[2]/3600.0);
  qfits_header_add_after(qh,"MJD-OBS","CD1_2",val," ",NULL);
  sprintf(val,"%e",t.a[1]/3600.0);
  qfits_header_add_after(qh,"MJD-OBS","CD1_1",val," ",NULL);
  sprintf(val,"%f",t.de0);
  qfits_header_add_after(qh,"MJD-OBS","CRVAL2",val," ",NULL);
  sprintf(val,"%f",t.ra0);
  qfits_header_add_after(qh,"MJD-OBS","CRVAL1",val," ",NULL);
  sprintf(val,"%f",t.y0);
  qfits_header_add_after(qh,"MJD-OBS","CRPIX2",val," ",NULL);
  sprintf(val,"%f",t.x0);
  qfits_header_add_after(qh,"MJD-OBS","CRPIX1",val," ",NULL);

  file=fopen(filename,"w");
  qfits_header_dump(qh,file);
  fclose(file);

  qfits_header_destroy(qh);

  qd.filename=filename;
  qd.npix=naxis1*naxis2*naxis3;
  qd.ptype=PTYPE_FLOAT;
  qd.fbuf=fbuf;
  qd.out_ptype=-32;

  qfits_pixdump(&qd);
  free(fbuf);

  return;
}

// Modify FITS keywords
void modify_fits_keywords(struct transformation t,char *filename) 
{
  char card[FITS_LINESZ+1];
  char key[FITS_LINESZ+1];
  char val[FITS_LINESZ+1];
  char com[FITS_LINESZ+1];

  sprintf(val,"%f",t.x0);
  keytuple2str(card,"CRPIX1",val,"");
  qfits_replace_card(filename,"CRPIX1",card);

  sprintf(val,"%f",t.y0);
  keytuple2str(card,"CRPIX2",val,"");
  qfits_replace_card(filename,"CRPIX2",card);

  sprintf(val,"%f",t.ra0);
  keytuple2str(card,"CRVAL1",val,"");
  qfits_replace_card(filename,"CRVAL1",card);

  sprintf(val,"%f",t.de0);
  keytuple2str(card,"CRVAL2",val,"");
  qfits_replace_card(filename,"CRVAL2",card);

  sprintf(val,"%e",t.a[1]/3600.0);
  keytuple2str(card,"CD1_1",val,"");
  qfits_replace_card(filename,"CD1_1",card);

  sprintf(val,"%e",t.a[2]/3600.0);
  keytuple2str(card,"CD1_2",val,"");
  qfits_replace_card(filename,"CD1_2",card);

  sprintf(val,"%e",t.b[1]/3600.0);
  keytuple2str(card,"CD2_1",val,"");
  qfits_replace_card(filename,"CD2_1",card);

  sprintf(val,"%e",t.b[2]/3600.0);
  keytuple2str(card,"CD2_2",val,"");
  qfits_replace_card(filename,"CD2_2",card);

  sprintf(val,"%e",t.xrms/3600.0);
  keytuple2str(card,"CRRES1",val,"");
  qfits_replace_card(filename,"CRRES1",card);

  sprintf(val,"%e",t.yrms/3600.0);
  keytuple2str(card,"CRRES2",val,"");
  qfits_replace_card(filename,"CRRES2",card);


  return;
}

// Precess a celestial position
void precess(double mjd0,double ra0,double de0,double mjd,double *ra,double *de)
{
  double t0,t;
  double zeta,z,theta;
  double a,b,c;

  // Angles in radians
  ra0*=D2R;
  de0*=D2R;

  // Time in centuries
  t0=(mjd0-51544.5)/36525.0;
  t=(mjd-mjd0)/36525.0;

  // Precession angles
  zeta=(2306.2181+1.39656*t0-0.000139*t0*t0)*t;
  zeta+=(0.30188-0.000344*t0)*t*t+0.017998*t*t*t;
  zeta*=D2R/3600.0;
  z=(2306.2181+1.39656*t0-0.000139*t0*t0)*t;
  z+=(1.09468+0.000066*t0)*t*t+0.018203*t*t*t;
  z*=D2R/3600.0;
  theta=(2004.3109-0.85330*t0-0.000217*t0*t0)*t;
  theta+=-(0.42665+0.000217*t0)*t*t-0.041833*t*t*t;
  theta*=D2R/3600.0;
  
  a=cos(de0)*sin(ra0+zeta);
  b=cos(theta)*cos(de0)*cos(ra0+zeta)-sin(theta)*sin(de0);
  c=sin(theta)*cos(de0)*cos(ra0+zeta)+cos(theta)*sin(de0);

  *ra=(atan2(a,b)+z)*R2D;
  *de=asin(c)*R2D;

  if (*ra<360.0)
    *ra+=360.0;
  if (*ra>360.0)
    *ra-=360.0;

  return;
}