Adding a few more tools
Cees Bassa
parent
9872a7e271
commit
7a943718d3
|
|
@ -0,0 +1,236 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include "cel.h"
|
||||
|
||||
#define LIM 80
|
||||
|
||||
void dec2sex(double,char *,int,int);
|
||||
double sex2dec(char *);
|
||||
void reverse(double,double,double,double,double *,double *);
|
||||
void forward(double,double,double,double,double *,double *);
|
||||
|
||||
int main(int argc,char *argv[])
|
||||
{
|
||||
int i;
|
||||
double ra1,de1,ra2,de2;
|
||||
double rx,ry;
|
||||
char sra[15],sde[15];
|
||||
|
||||
if (argc==1) {
|
||||
printf("Usage: %s <ra1> <de1> <ra2> <de2>\n",argv[0]);
|
||||
printf(" Computes angular offset\n");
|
||||
printf("Usage: %s -d <ra1> <de1> <dra> <dde>\n",argv[0]);
|
||||
printf(" Applies angular offset\n");
|
||||
printf("Usage: %s -x <ra1> <de1> <ra2> <de2>\n",argv[0]);
|
||||
printf(" Computes x-offset only\n");
|
||||
printf("Usage: %s -y <ra1> <de1> <ra2> <de2>\n",argv[0]);
|
||||
printf(" Computes y-offset only\n");
|
||||
return -1;
|
||||
}
|
||||
|
||||
if (strcmp(argv[1],"-d")==0) {
|
||||
if (strchr(argv[2],':')!=NULL)
|
||||
ra1=15.*sex2dec(argv[2]);
|
||||
else
|
||||
ra1=atof(argv[2]);
|
||||
if (strchr(argv[3],':')!=NULL)
|
||||
de1=sex2dec(argv[3]);
|
||||
else
|
||||
de1=atof(argv[3]);
|
||||
|
||||
rx=(double) atof(argv[4]);
|
||||
ry=(double) atof(argv[5]);
|
||||
|
||||
reverse(ra1,de1,rx,ry,&ra2,&de2);
|
||||
|
||||
dec2sex(ra2/15.0,sra,0,7);
|
||||
dec2sex(de2,sde,0,6);
|
||||
|
||||
printf("%s %s\n",sra,sde);
|
||||
} else if (strcmp(argv[1],"-x")==0) {
|
||||
if (strchr(argv[2],':')!=NULL)
|
||||
ra1=15.*sex2dec(argv[2]);
|
||||
else
|
||||
ra1=atof(argv[2]);
|
||||
if (strchr(argv[3],':')!=NULL)
|
||||
de1=sex2dec(argv[3]);
|
||||
else
|
||||
de1=atof(argv[3]);
|
||||
if (strchr(argv[4],':')!=NULL)
|
||||
ra2=15.*sex2dec(argv[4]);
|
||||
else
|
||||
ra2=atof(argv[4]);
|
||||
if (strchr(argv[5],':')!=NULL)
|
||||
de2=sex2dec(argv[5]);
|
||||
else
|
||||
de2=atof(argv[5]);
|
||||
|
||||
forward(ra1,de1,ra2,de2,&rx,&ry);
|
||||
|
||||
printf("%8.3f\n",rx);
|
||||
} else if (strcmp(argv[1],"-y")==0) {
|
||||
if (strchr(argv[2],':')!=NULL)
|
||||
ra1=15.*sex2dec(argv[2]);
|
||||
else
|
||||
ra1=atof(argv[2]);
|
||||
if (strchr(argv[3],':')!=NULL)
|
||||
de1=sex2dec(argv[3]);
|
||||
else
|
||||
de1=atof(argv[3]);
|
||||
if (strchr(argv[4],':')!=NULL)
|
||||
ra2=15.*sex2dec(argv[4]);
|
||||
else
|
||||
ra2=atof(argv[4]);
|
||||
if (strchr(argv[5],':')!=NULL)
|
||||
de2=sex2dec(argv[5]);
|
||||
else
|
||||
de2=atof(argv[5]);
|
||||
|
||||
forward(ra1,de1,ra2,de2,&rx,&ry);
|
||||
|
||||
printf("%8.3f\n",ry);
|
||||
} else {
|
||||
if (strchr(argv[1],':')!=NULL)
|
||||
ra1=15.*sex2dec(argv[1]);
|
||||
else
|
||||
ra1=atof(argv[1]);
|
||||
if (strchr(argv[2],':')!=NULL)
|
||||
de1=sex2dec(argv[2]);
|
||||
else
|
||||
de1=atof(argv[2]);
|
||||
if (strchr(argv[3],':')!=NULL)
|
||||
ra2=15.*sex2dec(argv[3]);
|
||||
else
|
||||
ra2=atof(argv[3]);
|
||||
if (strchr(argv[4],':')!=NULL)
|
||||
de2=sex2dec(argv[4]);
|
||||
else
|
||||
de2=atof(argv[4]);
|
||||
|
||||
forward(ra1,de1,ra2,de2,&rx,&ry);
|
||||
|
||||
printf("dRA:%8.3f dDE:%8.3f d:%8.3f\n",rx,ry,sqrt(rx*rx+ry*ry));
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Convert Decimal into Sexagesimal
|
||||
void dec2sex(double x,char *s,int f,int len)
|
||||
{
|
||||
int i;
|
||||
double sec,deg,min;
|
||||
char sign;
|
||||
char *form[]={":: ",",, ","hms"," "};
|
||||
|
||||
sign=(x<0 ? '-' : ' ');
|
||||
x=3600.*fabs(x);
|
||||
|
||||
sec=fmod(x,60.);
|
||||
x=(x-sec)/60.;
|
||||
min=fmod(x,60.);
|
||||
x=(x-min)/60.;
|
||||
// deg=fmod(x,60.);
|
||||
deg=x;
|
||||
|
||||
if (len==7) sprintf(s,"%c%02i%c%02i%c%07.4f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
|
||||
if (len==6) sprintf(s,"%c%02i%c%02i%c%06.3f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
|
||||
if (len==5) sprintf(s,"%c%02i%c%02i%c%05.2f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
|
||||
if (len==4) sprintf(s,"%c%02i%c%02i%c%04.1f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
|
||||
if (len==2) sprintf(s,"%c%02i%c%02i%c%02i%c",sign,(int) deg,form[f][0],(int) min,form[f][1],(int) floor(sec),form[f][2]);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
// Convert Sexagesimal into Decimal
|
||||
double sex2dec(char *s)
|
||||
{
|
||||
double x;
|
||||
float deg,min,sec;
|
||||
char t[LIM];
|
||||
|
||||
strcpy(t,s);
|
||||
|
||||
deg=fabs(atof(strtok(t," :")));
|
||||
min=fabs(atof(strtok(NULL," :")));
|
||||
sec=fabs(atof(strtok(NULL," :")));
|
||||
|
||||
x=(double) deg+(double) min/60.+(double) sec/3600.;
|
||||
if (s[0]=='-') x= -x;
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
// Get a RA and Decl from x and y
|
||||
void reverse(double ra0,double de0,double x,double y,double *ra,double *de)
|
||||
{
|
||||
int i;
|
||||
char pcode[4]="TAN";
|
||||
double phi,theta;
|
||||
struct celprm cel;
|
||||
struct prjprm prj;
|
||||
|
||||
x/=3600.;
|
||||
y/=3600.;
|
||||
|
||||
// Initialize Projection Parameters
|
||||
prj.flag=0;
|
||||
prj.r0=0.;
|
||||
for (i=0;i<10;prj.p[i++]=0.);
|
||||
|
||||
// Initialize Reference Angles
|
||||
cel.ref[0]=ra0;
|
||||
cel.ref[1]=de0;
|
||||
cel.ref[2]=999.;
|
||||
cel.ref[3]=999.;
|
||||
cel.flag=0.;
|
||||
|
||||
if (celset(pcode,&cel,&prj)) {
|
||||
printf("Error in Projection (celset)\n");
|
||||
return;
|
||||
} else {
|
||||
if (celrev(pcode,x,y,&prj,&phi,&theta,&cel,ra,de)) {
|
||||
printf("Error in Projection (celrev)\n");
|
||||
return;
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
// Get a x and y from a RA and Decl
|
||||
void forward(double ra0,double de0,double ra,double de,double *x,double *y)
|
||||
{
|
||||
int i;
|
||||
char pcode[4]="TAN";
|
||||
double phi,theta;
|
||||
struct celprm cel;
|
||||
struct prjprm prj;
|
||||
|
||||
// Initialize Projection Parameters
|
||||
prj.flag=0;
|
||||
prj.r0=0.;
|
||||
for (i=0;i<10;prj.p[i++]=0.);
|
||||
|
||||
// Initialize Reference Angles
|
||||
cel.ref[0]=ra0;
|
||||
cel.ref[1]=de0;
|
||||
cel.ref[2]=999.;
|
||||
cel.ref[3]=999.;
|
||||
cel.flag=0.;
|
||||
|
||||
if (celset(pcode,&cel,&prj)) {
|
||||
printf("Error in Projection (celset)\n");
|
||||
return;
|
||||
} else {
|
||||
if (celfwd(pcode,ra,de,&cel,&phi,&theta,&prj,x,y)) {
|
||||
printf("Error in Projection (celfwd)\n");
|
||||
return;
|
||||
}
|
||||
}
|
||||
*x *=3600.;
|
||||
*y *=3600.;
|
||||
|
||||
return;
|
||||
}
|
||||
|
|
@ -0,0 +1,580 @@
|
|||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include "qfits.h"
|
||||
#include "cpgplot.h"
|
||||
#include "cel.h"
|
||||
#include <gsl/gsl_multifit.h>
|
||||
|
||||
#define LIM 256
|
||||
#define NMAX 8192
|
||||
#define D2R M_PI/180.0
|
||||
#define R2D 180.0/M_PI
|
||||
|
||||
struct star {
|
||||
double ra,de;
|
||||
float pmra,pmde;
|
||||
float mag;
|
||||
};
|
||||
struct image {
|
||||
int naxis1,naxis2,naxis3;
|
||||
float *z;
|
||||
float zmin,zmax;
|
||||
double ra0,de0;
|
||||
float x0,y0;
|
||||
float a[3],b[3];
|
||||
double mjd;
|
||||
} img;
|
||||
struct catalog {
|
||||
int n;
|
||||
float x[NMAX],y[NMAX],mag[NMAX];
|
||||
double ra[NMAX],de[NMAX],rx[NMAX],ry[NMAX];
|
||||
int select[NMAX];
|
||||
};
|
||||
struct image read_fits(char *filename,int pnum);
|
||||
struct catalog read_pixel_catalog(char *filename);
|
||||
int fgetline(FILE *file,char *s,int lim);
|
||||
struct catalog read_astrometric_catalog(char *filename,float mmin,float sx,float sy,float angle);
|
||||
void forward(double ra0,double de0,double ra,double de,double *x,double *y);
|
||||
int select_nearest(struct catalog c,float x,float y);
|
||||
void lfit2d(float *x,float *y,float *z,int n,float *a);
|
||||
void fit_transformation(struct catalog cat,struct catalog ast,int nselect);
|
||||
struct catalog reread_astrometric_catalog(char *filename,float mmin);
|
||||
int match_catalogs(struct catalog *cat,struct catalog *ast,float rmax);
|
||||
|
||||
int main(int argc,char *argv[])
|
||||
{
|
||||
int i;
|
||||
float xmin,xmax,ymin,ymax,zmin,zmax;
|
||||
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};
|
||||
int redraw=1,plotcat=0,click=0,nselect=0;
|
||||
float x,y,width;
|
||||
char c,pixcat[LIM];
|
||||
struct catalog cat,ast;
|
||||
float sx,sy,q;
|
||||
char *env,starfile[128];
|
||||
float r,rmin=1.0,rmax=10.0,mmin=1.0,mmax=8.0,mag=8.0;
|
||||
|
||||
// Environment variables
|
||||
env=getenv("ST_DATADIR");
|
||||
sprintf(starfile,"%s/data/tycho2.dat",env);
|
||||
|
||||
// Read image
|
||||
img=read_fits(argv[1],0);
|
||||
|
||||
// Hard coded
|
||||
img.ra0=301.53;
|
||||
img.de0=8.07;
|
||||
sx=-10.0;
|
||||
sy=10.0;
|
||||
q=10.0;
|
||||
img.x0=0.5*(float) img.naxis1;
|
||||
img.y0=0.5*(float) img.naxis2;
|
||||
|
||||
// Read pixel catalog
|
||||
sprintf(pixcat,"%s.cat",argv[1]);
|
||||
cat=read_pixel_catalog(pixcat);
|
||||
|
||||
// Read astrometric catalog
|
||||
ast=read_astrometric_catalog(starfile,mag,sx,sy,q);
|
||||
|
||||
// Open server
|
||||
cpgopen("/xs");
|
||||
// cpgpap(0.,1.0);
|
||||
cpgask(0);
|
||||
cpgsch(0.8);
|
||||
|
||||
// Default limits
|
||||
width=(img.naxis1>img.naxis2) ? img.naxis1 : img.naxis2;
|
||||
xmin=0.5*(img.naxis1-width);
|
||||
xmax=0.5*(img.naxis1+width);
|
||||
ymin=0.5*(img.naxis2-width);
|
||||
ymax=0.5*(img.naxis2+width);
|
||||
zmin=img.zmin;
|
||||
zmax=img.zmax;
|
||||
|
||||
// Forever loop
|
||||
for (;;) {
|
||||
if (redraw==1) {
|
||||
cpgeras();
|
||||
|
||||
cpgsvp(0.1,0.95,0.1,0.95);
|
||||
cpgwnad(xmin,xmax,ymin,ymax);
|
||||
cpglab("x (pix)","y (pix)"," ");
|
||||
cpgsfs(2);
|
||||
cpgctab (heat_l,heat_r,heat_g,heat_b,5,1.0,0.5);
|
||||
|
||||
cpgimag(img.z,img.naxis1,img.naxis2,1,img.naxis1,1,img.naxis2,zmin,zmax,tr);
|
||||
cpgbox("BCTSNI",0.,0,"BCTSNI",0.,0);
|
||||
|
||||
// Plot catalog
|
||||
if (plotcat==1) {
|
||||
cpgsci(3);
|
||||
for (i=0;i<cat.n;i++) {
|
||||
if (cat.select[i]!=0)
|
||||
cpgpt1(cat.x[i],cat.y[i],6);
|
||||
else
|
||||
cpgpt1(cat.x[i],cat.y[i],4);
|
||||
}
|
||||
cpgsci(1);
|
||||
}
|
||||
|
||||
cpgsci(4);
|
||||
for (i=0;i<ast.n;i++) {
|
||||
r=rmax-(rmax-rmin)*(ast.mag[i]-mmin)/(mmax-mmin);
|
||||
|
||||
// Upscale for image size
|
||||
r*=img.naxis1/752.0;
|
||||
if (ast.select[i]!=0)
|
||||
cpgpt1(ast.x[i],ast.y[i],6);
|
||||
cpgcirc(ast.x[i],ast.y[i],r);
|
||||
}
|
||||
cpgsci(1);
|
||||
redraw=0;
|
||||
}
|
||||
|
||||
// Get cursor
|
||||
cpgcurs(&x,&y,&c);
|
||||
|
||||
// Quit
|
||||
if (c=='q')
|
||||
break;
|
||||
|
||||
// Select pixel catalog
|
||||
if (c=='a' && click==0) {
|
||||
i=select_nearest(cat,x,y);
|
||||
cat.select[i]=nselect+1;
|
||||
redraw=1;
|
||||
click=1;
|
||||
}
|
||||
|
||||
// Select catalog
|
||||
if (c=='b' && click==1) {
|
||||
i=select_nearest(ast,x,y);
|
||||
ast.select[i]=nselect+1;
|
||||
redraw=1;
|
||||
click=0;
|
||||
nselect++;
|
||||
}
|
||||
|
||||
// Center
|
||||
if (c=='c') {
|
||||
xmin=x-0.5*width;
|
||||
xmax=x+0.5*width;
|
||||
ymin=y-0.5*width;
|
||||
ymax=y+0.5*width;
|
||||
redraw=1;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Fit
|
||||
if (c=='f' && nselect>=3) {
|
||||
fit_transformation(cat,ast,nselect);
|
||||
ast=reread_astrometric_catalog(starfile,mag+1);
|
||||
redraw=1;
|
||||
}
|
||||
|
||||
// Zoom
|
||||
if (c=='z' || c=='+') {
|
||||
width/=1.25;
|
||||
xmin=x-0.5*width;
|
||||
xmax=x+0.5*width;
|
||||
ymin=y-0.5*width;
|
||||
ymax=y+0.5*width;
|
||||
redraw=1;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Match catalogs
|
||||
if (c=='m') {
|
||||
nselect=match_catalogs(&cat,&ast,10.0);
|
||||
redraw=1;
|
||||
}
|
||||
|
||||
// Unzoom
|
||||
if (c=='x' || c=='+' || c=='=') {
|
||||
width*=1.25;
|
||||
xmin=x-0.5*width;
|
||||
xmax=x+0.5*width;
|
||||
ymin=y-0.5*width;
|
||||
ymax=y+0.5*width;
|
||||
redraw=1;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Plot catalog
|
||||
if (c=='p') {
|
||||
plotcat = (plotcat==1) ? 0 : 1;
|
||||
redraw=1;
|
||||
}
|
||||
|
||||
// Reset
|
||||
if (c=='r') {
|
||||
width=(img.naxis1>img.naxis2) ? img.naxis1 : img.naxis2;
|
||||
xmin=0.5*(img.naxis1-width);
|
||||
xmax=0.5*(img.naxis1+width);
|
||||
ymin=0.5*(img.naxis2-width);
|
||||
ymax=0.5*(img.naxis2+width);
|
||||
redraw=1;
|
||||
continue;
|
||||
}
|
||||
}
|
||||
|
||||
cpgend();
|
||||
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Read fits image
|
||||
struct image read_fits(char *filename,int pnum)
|
||||
{
|
||||
int i,j,k,l,m;
|
||||
qfitsloader ql;
|
||||
char key[FITS_LINESZ+1] ;
|
||||
struct image img;
|
||||
float s1,s2,avg,std;
|
||||
|
||||
// Set plane
|
||||
ql.xtnum = 0;
|
||||
ql.pnum = pnum;
|
||||
|
||||
// Set loadtype
|
||||
ql.ptype = PTYPE_FLOAT;
|
||||
|
||||
// Set filename
|
||||
ql.filename=filename;
|
||||
|
||||
// Image size
|
||||
img.naxis1=atoi(qfits_query_hdr(filename,"NAXIS1"));
|
||||
img.naxis2=atoi(qfits_query_hdr(filename,"NAXIS2"));
|
||||
// img.mjd=atof(qfits_query_hdr(filename,"MJD-OBS"));
|
||||
|
||||
// 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");
|
||||
|
||||
// Allocate image memory
|
||||
img.z=(float *) malloc(sizeof(float) * img.naxis1*img.naxis2);
|
||||
|
||||
// Fill z array
|
||||
for (i=0,l=0,m=0;i<img.naxis1;i++) {
|
||||
for (j=0;j<img.naxis2;j++) {
|
||||
img.z[l]=ql.fbuf[l];
|
||||
l++;
|
||||
}
|
||||
}
|
||||
|
||||
// Get levels
|
||||
for (i=0,s1=0.0,s2=0.0;i<img.naxis1*img.naxis2;i++) {
|
||||
s1+=img.z[i];
|
||||
s2+=img.z[i]*img.z[i];
|
||||
}
|
||||
avg=s1/(float) (img.naxis1*img.naxis2);
|
||||
std=sqrt(s2/(float) (img.naxis1*img.naxis2)-avg*avg);
|
||||
img.zmin=avg-4.0*std;
|
||||
img.zmax=avg+16.0*std;
|
||||
|
||||
return img;
|
||||
}
|
||||
|
||||
// Read pixel catalog
|
||||
struct catalog read_pixel_catalog(char *filename)
|
||||
{
|
||||
int i=0;
|
||||
FILE *file;
|
||||
char line[LIM];
|
||||
struct catalog c;
|
||||
|
||||
// Read catalog
|
||||
file=fopen(filename,"r");
|
||||
if (file==NULL) {
|
||||
fprintf(stderr,"%s not found!\n",filename);
|
||||
exit(0);
|
||||
}
|
||||
while (fgetline(file,line,LIM)>0) {
|
||||
if (i>=NMAX)
|
||||
break;
|
||||
if (strstr(line,"#")!=NULL)
|
||||
continue;
|
||||
sscanf(line,"%f %f %f",&c.x[i],&c.y[i],&c.mag[i]);
|
||||
c.select[i]=0;
|
||||
i++;
|
||||
}
|
||||
fclose(file);
|
||||
c.n=i;
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
struct catalog read_astrometric_catalog(char *filename,float mmin,float sx,float sy,float angle)
|
||||
{
|
||||
int i=0;
|
||||
FILE *file;
|
||||
char line[LIM];
|
||||
struct catalog c;
|
||||
double rx,ry,x,y,ra,de;
|
||||
struct star s;
|
||||
double d,dx,dy;
|
||||
double mjd0=51544.5;
|
||||
|
||||
file=fopen(filename,"rb");
|
||||
if (file==NULL) {
|
||||
fprintf(stderr,"%s not found!\n",filename);
|
||||
exit(0);
|
||||
}
|
||||
while (!feof(file)) {
|
||||
fread(&s,sizeof(struct star),1,file);
|
||||
if (s.mag>mmin)
|
||||
continue;
|
||||
forward(img.ra0,img.de0,s.ra,s.de,&rx,&ry);
|
||||
x=img.x0+1.0/sx*(cos(angle*D2R)*rx+sin(angle*D2R)*ry);
|
||||
y=img.y0+1.0/sy*(-sin(angle*D2R)*rx+cos(angle*D2R)*ry);
|
||||
if (x>0.0 && x<(float) img.naxis1 && y>0.0 && y<(float) img.naxis2) {
|
||||
c.x[i]=x;
|
||||
c.y[i]=y;
|
||||
c.rx[i]=rx;
|
||||
c.ry[i]=ry;
|
||||
c.ra[i]=s.ra;
|
||||
c.de[i]=s.de;
|
||||
c.mag[i]=s.mag;
|
||||
c.select[i]=0;
|
||||
|
||||
if (i>=NMAX)
|
||||
break;
|
||||
i++;
|
||||
}
|
||||
}
|
||||
fclose(file);
|
||||
c.n=i;
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
// Get a x and y from a RA and Decl
|
||||
void forward(double ra0,double de0,double ra,double de,double *x,double *y)
|
||||
{
|
||||
int i;
|
||||
char pcode[4]="STG";
|
||||
double phi,theta;
|
||||
struct celprm cel;
|
||||
struct prjprm prj;
|
||||
|
||||
// Initialize Projection Parameters
|
||||
prj.flag=0;
|
||||
prj.r0=0.;
|
||||
for (i=0;i<10;prj.p[i++]=0.);
|
||||
|
||||
// Initialize Reference Angles
|
||||
cel.ref[0]=ra0;
|
||||
cel.ref[1]=de0;
|
||||
cel.ref[2]=999.;
|
||||
cel.ref[3]=999.;
|
||||
cel.flag=0.;
|
||||
|
||||
if (celset(pcode,&cel,&prj)) {
|
||||
printf("Error in Projection (celset)\n");
|
||||
return;
|
||||
} else {
|
||||
if (celfwd(pcode,ra,de,&cel,&phi,&theta,&prj,x,y)) {
|
||||
printf("Error in Projection (celfwd)\n");
|
||||
return;
|
||||
}
|
||||
}
|
||||
*x *=3600.;
|
||||
*y *=3600.;
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
// Select nearest object
|
||||
int select_nearest(struct catalog c,float x,float y)
|
||||
{
|
||||
int i,imin;
|
||||
float r,rmin;
|
||||
|
||||
for (i=0;i<c.n;i++) {
|
||||
r=sqrt(pow(x-c.x[i],2)+pow(y-c.y[i],2));
|
||||
if (i==0 || r<rmin) {
|
||||
imin=i;
|
||||
rmin=r;
|
||||
}
|
||||
}
|
||||
|
||||
return imin;
|
||||
}
|
||||
|
||||
// Fit transformation
|
||||
void fit_transformation(struct catalog cat,struct catalog ast,int nselect)
|
||||
{
|
||||
int i,j;
|
||||
float *x,*y,*rx,*ry;
|
||||
|
||||
x=(float *) malloc(sizeof(float)*nselect);
|
||||
y=(float *) malloc(sizeof(float)*nselect);
|
||||
rx=(float *) malloc(sizeof(float)*nselect);
|
||||
ry=(float *) malloc(sizeof(float)*nselect);
|
||||
|
||||
for (i=0;i<nselect;i++) {
|
||||
for (j=0;j<cat.n;j++) {
|
||||
if (cat.select[j]==i+1) {
|
||||
x[i]=cat.x[j]-img.x0;
|
||||
y[i]=cat.y[j]-img.y0;
|
||||
}
|
||||
}
|
||||
for (j=0;j<ast.n;j++) {
|
||||
if (ast.select[j]==i+1) {
|
||||
rx[i]=ast.rx[j];
|
||||
ry[i]=ast.ry[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
lfit2d(x,y,rx,nselect,img.a);
|
||||
lfit2d(x,y,ry,nselect,img.b);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
// Linear 2D fit
|
||||
void lfit2d(float *x,float *y,float *z,int n,float *a)
|
||||
{
|
||||
int i;
|
||||
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;
|
||||
}
|
||||
|
||||
// Read astrometric catalog
|
||||
struct catalog reread_astrometric_catalog(char *filename,float mmin)
|
||||
{
|
||||
int i=0;
|
||||
FILE *file;
|
||||
char line[LIM];
|
||||
struct catalog c;
|
||||
double rx,ry,x,y;
|
||||
struct star s;
|
||||
double d,dx,dy,ra,de;
|
||||
double mjd0=51544.5;
|
||||
|
||||
file=fopen(filename,"rb");
|
||||
while (!feof(file)) {
|
||||
fread(&s,sizeof(struct star),1,file);
|
||||
if (s.mag>mmin)
|
||||
continue;
|
||||
forward(img.ra0,img.de0,s.ra,s.de,&rx,&ry);
|
||||
dx=rx-img.a[0];
|
||||
dy=ry-img.b[0];
|
||||
d=img.a[1]*img.b[2]-img.a[2]*img.b[1];
|
||||
x=(img.b[2]*dx-img.a[2]*dy)/d+img.x0;
|
||||
y=(img.a[1]*dy-img.b[1]*dx)/d+img.y0;
|
||||
if (x>0.0 && x<img.naxis1 && y>0.0 && y<img.naxis2) {
|
||||
c.x[i]=x;
|
||||
c.y[i]=y;
|
||||
c.rx[i]=rx;
|
||||
c.ry[i]=ry;
|
||||
c.ra[i]=s.ra;
|
||||
c.de[i]=s.de;
|
||||
c.mag[i]=s.mag;
|
||||
c.select[i]=0;
|
||||
i++;
|
||||
}
|
||||
}
|
||||
fclose(file);
|
||||
c.n=i;
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
int match_catalogs(struct catalog *cat,struct catalog *ast,float rmax)
|
||||
{
|
||||
int i,j,jmin,n,flag=0;
|
||||
float r,rmin;
|
||||
FILE *file;
|
||||
|
||||
// Reset
|
||||
for (i=0;i<cat->n;i++)
|
||||
cat->select[i]=0;
|
||||
for (i=0;i<ast->n;i++)
|
||||
ast->select[i]=0;
|
||||
|
||||
file=fopen("out.dat","w");
|
||||
for (i=0,n=0;i<cat->n;i++) {
|
||||
for (j=0,flag=0;j<ast->n;j++) {
|
||||
if (ast->select[j]!=0)
|
||||
continue;
|
||||
r=sqrt(pow(cat->x[i]-ast->x[j],2)+pow(cat->y[i]-ast->y[j],2));
|
||||
if (flag==0 || r<rmin) {
|
||||
rmin=r;
|
||||
jmin=j;
|
||||
flag=1;
|
||||
}
|
||||
}
|
||||
if (rmin<rmax) {
|
||||
fprintf(file,"%10.4f %10.4f %10.6f %10.6f\n",cat->x[i]-img.x0,cat->y[i]-img.y0,ast->ra[jmin],ast->de[jmin]);
|
||||
cat->select[i]=n+1;
|
||||
ast->select[jmin]=n+1;
|
||||
n++;
|
||||
}
|
||||
}
|
||||
fclose(file);
|
||||
|
||||
printf("%d stars matched\n",n);
|
||||
return n;
|
||||
}
|
||||
|
|
@ -0,0 +1,57 @@
|
|||
/* Compute sexagesimal from decimal */
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
|
||||
int main(int argc,char *argv[])
|
||||
{
|
||||
int h=0;
|
||||
double sec,deg,min;
|
||||
char c,sign,d=' ';
|
||||
double x;
|
||||
char s[14];
|
||||
|
||||
if (argc==1) {
|
||||
printf("Usage: dec2sex [options] <x>\n\nCompute sexagesimal from decimal input x.\n");
|
||||
printf("Options: -r gives hours instead of degrees\n");
|
||||
printf(" -d<a> uses character a as delimiter\n");
|
||||
printf(" -h uses hms as delimiters\n");
|
||||
printf(" -s always prints sign\n\n");
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
// Get Decimal Value
|
||||
x=(double) atof(argv[--argc]);
|
||||
sign=(x<0 ? '-' : ' ');
|
||||
|
||||
// Get Options
|
||||
while (--argc > 0 && (*++argv)[0] == '-') {
|
||||
while (c = *++argv[0]) {
|
||||
if (c == 'r')
|
||||
x /= 15.;
|
||||
if (c == 's')
|
||||
sign=(x<0 ? '-' : '+');
|
||||
if (c == 'd')
|
||||
if (strlen(argv[0])!=1)
|
||||
d=(*argv)[1];
|
||||
if (c == 'h')
|
||||
h++;
|
||||
}
|
||||
}
|
||||
|
||||
x=3600.*fabs(x);
|
||||
sec=fmod(x,60.);
|
||||
x=(x-sec)/60.;
|
||||
min=fmod(x,60.);
|
||||
x=(x-min)/60.;
|
||||
deg=x;
|
||||
|
||||
if (h==0)
|
||||
printf("%c%02i%c%02i%c%06.3f\n",sign,(int) deg,d,(int) min,d,(float) sec);
|
||||
else
|
||||
printf("%c%02ih%02im%06.3fs\n",sign,(int) deg,(int) min,(float) sec);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,289 @@
|
|||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <qfits.h>
|
||||
#include <jpeglib.h>
|
||||
#include <libexif/exif-data.h>
|
||||
#include <getopt.h>
|
||||
|
||||
struct image {
|
||||
int nx,ny,nz;
|
||||
float *z;
|
||||
float xmin,xmax,ymin,ymax;
|
||||
double tr[6];
|
||||
double mjd;
|
||||
char nfd[32];
|
||||
};
|
||||
struct image read_jpg(char *filename);
|
||||
void write_fits(struct image img,char *filename);
|
||||
double date2mjd(int year,int month,double day);
|
||||
double nfd2mjd(char *date);
|
||||
|
||||
int main(int argc,char *argv[])
|
||||
{
|
||||
int arg;
|
||||
struct image img;
|
||||
char infile[64],outfile[64],nfd[32];
|
||||
double mjd;
|
||||
|
||||
// Decode options
|
||||
while ((arg=getopt(argc,argv,"i:t:o:"))!=-1) {
|
||||
switch(arg) {
|
||||
|
||||
case 'i':
|
||||
strcpy(infile,optarg);
|
||||
break;
|
||||
|
||||
case 'o':
|
||||
strcpy(outfile,optarg);
|
||||
break;
|
||||
|
||||
case 't':
|
||||
strcpy(nfd,optarg);
|
||||
mjd=nfd2mjd(nfd);
|
||||
break;
|
||||
|
||||
default:
|
||||
return 0;
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if (infile!=NULL)
|
||||
img=read_jpg(infile);
|
||||
|
||||
if (nfd!=NULL) {
|
||||
strcpy(img.nfd,nfd);
|
||||
img.mjd=mjd;
|
||||
}
|
||||
|
||||
if (outfile!=NULL)
|
||||
write_fits(img,outfile);
|
||||
|
||||
free(img.z);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
struct image read_jpg(char *filename)
|
||||
{
|
||||
int i=0,j,k,l,m;
|
||||
unsigned long location=0;
|
||||
struct image img;
|
||||
struct jpeg_decompress_struct cinfo;
|
||||
struct jpeg_error_mgr jerr;
|
||||
JSAMPROW row_pointer[1];
|
||||
unsigned char *raw_image=NULL;
|
||||
FILE *file;
|
||||
ExifData *ed;
|
||||
ExifEntry *entry;
|
||||
|
||||
// Open file
|
||||
file=fopen(filename,"rb");
|
||||
if (!file)
|
||||
perror("Error opening file");
|
||||
|
||||
// Get header info, decompress
|
||||
cinfo.err=jpeg_std_error(&jerr);
|
||||
jpeg_create_decompress(&cinfo);
|
||||
jpeg_stdio_src(&cinfo,file);
|
||||
jpeg_read_header(&cinfo,TRUE);
|
||||
jpeg_start_decompress(&cinfo);
|
||||
|
||||
// Allocate memory
|
||||
raw_image=(unsigned char *) malloc(cinfo.output_width*cinfo.output_height*cinfo.num_components);
|
||||
|
||||
// Read image, one scan at a time
|
||||
row_pointer[0]=(unsigned char *) malloc(cinfo.output_width*cinfo.num_components);
|
||||
while(cinfo.output_scanline<cinfo.image_height) {
|
||||
jpeg_read_scanlines(&cinfo,row_pointer,1);
|
||||
for(i=0;i<cinfo.image_width*cinfo.num_components;i++)
|
||||
raw_image[location++]=row_pointer[0][i];
|
||||
}
|
||||
// wrap up decompression, destroy objects, free pointers and close open files
|
||||
jpeg_finish_decompress(&cinfo);
|
||||
jpeg_destroy_decompress(&cinfo);
|
||||
|
||||
// Copy image to image struct
|
||||
img.nx=cinfo.image_width;
|
||||
img.ny=cinfo.image_height;
|
||||
img.nz=cinfo.num_components;
|
||||
img.z=(float *) malloc(sizeof(float)*img.nx*img.ny*img.nz);
|
||||
|
||||
// Fill image
|
||||
for (i=0;i<img.nx;i++) {
|
||||
for (j=0;j<img.ny;j++) {
|
||||
k=i+(img.ny-j-1)*img.nx;
|
||||
img.z[k]=0.0;
|
||||
for (l=0;l<img.nz;l++) {
|
||||
m=img.nz*(i+img.nx*j)+l;
|
||||
img.z[k]+=(float) raw_image[m];
|
||||
}
|
||||
img.z[k]/=3.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Limits
|
||||
img.xmin=0.0;
|
||||
img.xmax=(float) img.nx;
|
||||
img.ymin=0.0;
|
||||
img.ymax=(float) img.ny;
|
||||
|
||||
// Translation matrix
|
||||
img.tr[1]=(img.xmax-img.xmin)/(float) img.nx;
|
||||
img.tr[2]=0.0;
|
||||
img.tr[0]=img.xmin-0.5*img.tr[1];
|
||||
img.tr[4]=0.0;
|
||||
img.tr[5]=(img.ymax-img.ymin)/(float) img.ny;
|
||||
img.tr[3]=img.ymin-0.5*img.tr[5];
|
||||
|
||||
// Free allocated memory
|
||||
free(row_pointer[0]);
|
||||
free(raw_image);
|
||||
|
||||
// Close file
|
||||
fclose(file);
|
||||
|
||||
// Get exif info
|
||||
ed=exif_data_new_from_file(filename);
|
||||
if (!ed) {
|
||||
printf("File not readable or no EXIF data in file %s\n",filename);
|
||||
} else {
|
||||
entry=exif_content_get_entry(ed->ifd[0],EXIF_TAG_DATE_TIME);
|
||||
exif_entry_get_value(entry,img.nfd, sizeof(img.nfd));
|
||||
img.nfd[4]='-';
|
||||
img.nfd[7]='-';
|
||||
img.nfd[10]='T';
|
||||
img.nfd[20]='\0';
|
||||
img.mjd=nfd2mjd(img.nfd);
|
||||
}
|
||||
|
||||
return img;
|
||||
}
|
||||
|
||||
// Write fits file
|
||||
void write_fits(struct image img,char *filename)
|
||||
{
|
||||
int i,j,k,l;
|
||||
int *ibuf;
|
||||
qfitsdumper qd;
|
||||
qfits_header *qh;
|
||||
char key[FITS_LINESZ+1] ;
|
||||
char val[FITS_LINESZ+1] ;
|
||||
char com[FITS_LINESZ+1] ;
|
||||
char lin[FITS_LINESZ+1] ;
|
||||
FILE *file;
|
||||
|
||||
// Create FITS header
|
||||
qh=qfits_header_default();
|
||||
|
||||
// Add stuff
|
||||
qfits_header_add(qh,"BITPIX","16"," ",NULL);
|
||||
qfits_header_add(qh,"NAXIS","2"," ",NULL);
|
||||
sprintf(val,"%i",img.nx);
|
||||
qfits_header_add(qh,"NAXIS1",val," ",NULL);
|
||||
sprintf(val,"%i",img.ny);
|
||||
qfits_header_add(qh,"NAXIS2",val," ",NULL);
|
||||
qfits_header_add(qh,"BSCALE","1.0"," ",NULL);
|
||||
qfits_header_add(qh,"BZERO","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"DATAMAX","255.0"," ",NULL);
|
||||
qfits_header_add(qh,"DATAMIN","0.0"," ",NULL);
|
||||
|
||||
// Astrometry keywors
|
||||
sprintf(val,"%f",img.nx/2.0);
|
||||
qfits_header_add(qh,"CRPIX1",val," ",NULL);
|
||||
sprintf(val,"%f",img.ny/2.0);
|
||||
qfits_header_add(qh,"CRPIX2",val," ",NULL);
|
||||
qfits_header_add(qh,"CRVAL1","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CRVAL2","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CD1_1","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CD1_2","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CD2_1","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CD2_2","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CTYPE1","'RA---TAN'"," ",NULL);
|
||||
qfits_header_add(qh,"CTYPE2","'DEC--TAN'"," ",NULL);
|
||||
qfits_header_add(qh,"CUNIT1","'deg'"," ",NULL);
|
||||
qfits_header_add(qh,"CUNIT2","'deg'"," ",NULL);
|
||||
qfits_header_add(qh,"CRRES1","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"CRRES2","0.0"," ",NULL);
|
||||
qfits_header_add(qh,"EQUINOX","2000.0"," ",NULL);
|
||||
qfits_header_add(qh,"RADECSYS","ICRS"," ",NULL);
|
||||
sprintf(val,"%s",img.nfd);
|
||||
qfits_header_add(qh,"DATE-OBS",val," ",NULL);
|
||||
sprintf(val,"%lf",img.mjd);
|
||||
qfits_header_add(qh,"MJD-OBS",val," ",NULL);
|
||||
qfits_header_add(qh,"COSPAR","4171"," ",NULL);
|
||||
qfits_header_add(qh,"EXPTIME","10.060"," ",NULL);
|
||||
qfits_header_add(qh,"OBSERVER","Cees Bassa"," ",NULL);
|
||||
|
||||
// Dump fitsheader
|
||||
// qfits_header_dump(qh,stdout);
|
||||
|
||||
// Dump to file
|
||||
file=fopen(filename,"w");
|
||||
qfits_header_dump(qh,file);
|
||||
fclose(file);
|
||||
|
||||
|
||||
// Fill buffer
|
||||
ibuf=malloc(img.nx*img.ny*sizeof(int));
|
||||
for (i=0,l=0;i<img.nx;i++) {
|
||||
for (j=img.ny-1;j>=0;j--) {
|
||||
ibuf[l]=(int) img.z[l];
|
||||
|
||||
l++;
|
||||
}
|
||||
}
|
||||
|
||||
// Set parameters
|
||||
qd.filename=filename;
|
||||
qd.npix=img.nx*img.ny;
|
||||
qd.ptype=PTYPE_INT;
|
||||
qd.ibuf=ibuf;
|
||||
qd.out_ptype=BPP_16_SIGNED;
|
||||
|
||||
// Dump
|
||||
qfits_pixdump(&qd);
|
||||
|
||||
free(ibuf);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
// Compute Julian Day from Date
|
||||
double date2mjd(int year,int month,double day)
|
||||
{
|
||||
int a,b;
|
||||
double jd;
|
||||
|
||||
if (month<3) {
|
||||
year--;
|
||||
month+=12;
|
||||
}
|
||||
|
||||
a=floor(year/100.);
|
||||
b=2.-a+floor(a/4.);
|
||||
|
||||
if (year<1582) b=0;
|
||||
if (year==1582 && month<10) b=0;
|
||||
if (year==1852 && month==10 && day<=4) b=0;
|
||||
|
||||
jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
|
||||
|
||||
return jd-2400000.5;
|
||||
}
|
||||
|
||||
// nfd2mjd
|
||||
double nfd2mjd(char *date)
|
||||
{
|
||||
int year,month,day,hour,min,sec;
|
||||
double mjd,dday;
|
||||
|
||||
sscanf(date,"%04d-%02d-%02dT%02d:%02d:%02d",&year,&month,&day,&hour,&min,&sec);
|
||||
|
||||
dday=day+hour/24.0+min/1440.0+sec/86400.0;
|
||||
mjd=date2mjd(year,month,dday);
|
||||
|
||||
return mjd;
|
||||
}
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
// Compute decimal from sexagesimal
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
|
||||
// Usage: sex2dec [options] hh:mm:ss.ss
|
||||
int main(int argc,char *argv[])
|
||||
{
|
||||
double x;
|
||||
int sign=1;
|
||||
float deg,min,sec;
|
||||
char t[20],c;
|
||||
|
||||
if (argc==1) {
|
||||
printf("Usage: sex2dec -r <hh:mm:ss.sss>\n -or- <ddd:mm:ss.sss>\nCompute sexagesimal from decimal input x.\n");
|
||||
printf("Options: -r Converts hours into degrees\n");
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
// Get Sexagesimal string
|
||||
strcpy(t,argv[--argc]);
|
||||
if (t[0]=='-') sign=-1;
|
||||
|
||||
deg=fabs(atof(strtok(t," :,")));
|
||||
min=fabs(atof(strtok(NULL," :,")));
|
||||
sec=fabs(atof(strtok(NULL," :,")));
|
||||
|
||||
x=(double) deg+(double) min/60.+(double) sec/3600.;
|
||||
|
||||
// Get Options
|
||||
while (--argc > 0 && (*++argv)[0] == '-') {
|
||||
while (c = *++argv[0]) {
|
||||
if (c == 'r')
|
||||
x *= 15.;
|
||||
}
|
||||
}
|
||||
printf("%lf\n",sign*x);
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue