spacecruft
/
sattools
1
0
Fork 0
Code Releases Activity
sattools/src/wcsfit.c

389 lines
9.1 KiB
C
Raw Permalink Blame History

#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;
}
View Git Blame Copy Permalink