359 lines
8.3 KiB
C
359 lines
8.3 KiB
C
|
#include <stdio.h>
|
||
|
#include <string.h>
|
||
|
#include <stdlib.h>
|
||
|
#include <math.h>
|
||
|
#include <wcslib/cel.h>
|
||
|
#include <cpgplot.h>
|
||
|
#include "qfits.h"
|
||
|
#include <gsl/gsl_multifit.h>
|
||
|
|
||
|
#define LIM 256
|
||
|
#define D2R M_PI/180.0
|
||
|
#define R2D 180.0/M_PI
|
||
|
#define NMAX 4096
|
||
|
|
||
|
struct catalog {
|
||
|
int n;
|
||
|
float x[NMAX],y[NMAX];
|
||
|
double ra[NMAX],de[NMAX];
|
||
|
double rx[NMAX],ry[NMAX];
|
||
|
float xres[NMAX],yres[NMAX],res[NMAX];
|
||
|
float xrms,yrms,rms;
|
||
|
int usage[NMAX];
|
||
|
};
|
||
|
struct image {
|
||
|
int 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 ra0,de0;
|
||
|
float a[3],b[3];
|
||
|
float x0,y0;
|
||
|
};
|
||
|
int fgetline(FILE *file,char *s,int lim);
|
||
|
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);
|
||
|
struct catalog read_catalog(char *filename);
|
||
|
void lfit2d(float *x,float *y,float *z,int n,float *a);
|
||
|
void add_fits_keywords(struct transformation t,char *filename);
|
||
|
struct image read_fits(char *filename);
|
||
|
|
||
|
// 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);
|
||
|
|
||
|
return;
|
||
|
}
|
||
|
|
||
|
|
||
|
int main(int argc,char *argv[])
|
||
|
{
|
||
|
int i,j,k,l,m;
|
||
|
struct catalog c;
|
||
|
struct transformation t;
|
||
|
double ra0,de0;
|
||
|
float rmsmin;
|
||
|
float x[NMAX],y[NMAX],rx[NMAX],ry[NMAX];
|
||
|
struct image img;
|
||
|
char filename[128];
|
||
|
|
||
|
if (argc==1)
|
||
|
strcpy(filename,"test.fits");
|
||
|
else if (argc==2)
|
||
|
strcpy(filename,argv[1]);
|
||
|
|
||
|
img=read_fits(filename);
|
||
|
printf("files read\n");
|
||
|
c=read_catalog("out.dat");
|
||
|
printf("files read\n");
|
||
|
|
||
|
// Initial fit
|
||
|
t.ra0=c.ra[0];
|
||
|
t.de0=c.de[0];
|
||
|
t.x0=(float) img.naxis1/2.0;
|
||
|
t.y0=(float) img.naxis2/2.0;
|
||
|
|
||
|
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]=c.rx[i];
|
||
|
ry[k]=c.ry[i];
|
||
|
k++;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
// Fit
|
||
|
lfit2d(x,y,rx,k,t.a);
|
||
|
lfit2d(x,y,ry,k,t.b);
|
||
|
printf("%f %f %f %f %f %f %f %f\n",t.ra0,t.de0,t.a[0],t.a[1],t.a[2],t.b[0],t.b[1],t.b[2]);
|
||
|
|
||
|
// 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,c.xrms=0.0,c.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]);
|
||
|
printf("%12.4f %12.4f %12.4f %12.4f %10.4f %10.4f\n",c.x[i],c.y[i],c.rx[i],c.ry[i],c.xres[i],c.yres[i]);
|
||
|
c.res[i]=sqrt(c.xres[i]*c.xres[i]+c.yres[i]*c.yres[i]);
|
||
|
c.xrms+=c.xres[i]*c.xres[i];
|
||
|
c.yrms+=c.yres[i]*c.yres[i];
|
||
|
c.rms+=c.xres[i]*c.xres[i]+c.yres[i]*c.yres[i];
|
||
|
m++;
|
||
|
}
|
||
|
}
|
||
|
c.xrms=sqrt(c.xrms/(float) m);
|
||
|
c.yrms=sqrt(c.yrms/(float) m);
|
||
|
c.rms=sqrt(c.rms/(float) m);
|
||
|
|
||
|
// Deselect outliers
|
||
|
for (i=0;i<c.n;i++) {
|
||
|
if (c.res[i]>2*c.rms)
|
||
|
c.usage[i]=0;
|
||
|
}
|
||
|
}
|
||
|
printf("%12.8lf %10.6lf %10.6lf %8.4f %8.4f %8.4f %8.4f\n",img.mjd,t.ra0,t.de0,t.a[1],t.a[2],t.b[1],t.b[2]);
|
||
|
printf("%d/%d %f %f %f\n",m,c.n,c.xrms,c.yrms,c.rms);
|
||
|
|
||
|
// add_fits_keywords(t,"test.fits");
|
||
|
modify_fits_keywords(t,filename);
|
||
|
|
||
|
return 0;
|
||
|
}
|
||
|
|
||
|
// 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;
|
||
|
}
|
||
|
|
||
|
// Read catalog
|
||
|
struct catalog read_catalog(char *filename)
|
||
|
{
|
||
|
int i=0;
|
||
|
char line[LIM];
|
||
|
FILE *file;
|
||
|
struct catalog c;
|
||
|
|
||
|
file=fopen(filename,"r");
|
||
|
while (fgetline(file,line,LIM)>0) {
|
||
|
sscanf(line,"%f %f %lf %lf",&c.x[i],&c.y[i],&c.ra[i],&c.de[i]);
|
||
|
c.usage[i]=1;
|
||
|
|
||
|
i++;
|
||
|
}
|
||
|
fclose(file);
|
||
|
c.n=i;
|
||
|
|
||
|
return c;
|
||
|
}
|
||
|
|
||
|
// 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++;
|
||
|
}
|
||
|
}
|
||
|
}
|
||
|
|
||
|
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;
|
||
|
}
|
||
|
|
||
|
// 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.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"));
|
||
|
img.mjd=0.0;
|
||
|
|
||
|
|
||
|
return img;
|
||
|
}
|