From 7a943718d3f4409662cb05d234d8baedd1920267 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Mon, 9 Sep 2013 08:36:40 +0200 Subject: [PATCH] Adding a few more tools --- angular.c | 236 +++++++++++++++++++++ calibrate.c | 580 ++++++++++++++++++++++++++++++++++++++++++++++++++++ dec2sex.c | 57 ++++++ jpg2fits.c | 289 ++++++++++++++++++++++++++ sex2dec.c | 42 ++++ 5 files changed, 1204 insertions(+) create mode 100644 angular.c create mode 100644 calibrate.c create mode 100644 dec2sex.c create mode 100644 jpg2fits.c create mode 100644 sex2dec.c diff --git a/angular.c b/angular.c new file mode 100644 index 0000000..a00f275 --- /dev/null +++ b/angular.c @@ -0,0 +1,236 @@ +#include +#include +#include +#include +#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 \n",argv[0]); + printf(" Computes angular offset\n"); + printf("Usage: %s -d \n",argv[0]); + printf(" Applies angular offset\n"); + printf("Usage: %s -x \n",argv[0]); + printf(" Computes x-offset only\n"); + printf("Usage: %s -y \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; +} diff --git a/calibrate.c b/calibrate.c new file mode 100644 index 0000000..685613b --- /dev/null +++ b/calibrate.c @@ -0,0 +1,580 @@ +#include +#include +#include +#include +#include "qfits.h" +#include "cpgplot.h" +#include "cel.h" +#include + +#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=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;i0) { + 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;immin) + 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 && x0.0 && yn;i++) + cat->select[i]=0; + for (i=0;in;i++) + ast->select[i]=0; + + file=fopen("out.dat","w"); + for (i=0,n=0;in;i++) { + for (j=0,flag=0;jn;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 || rx[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; +} diff --git a/dec2sex.c b/dec2sex.c new file mode 100644 index 0000000..6c2ce43 --- /dev/null +++ b/dec2sex.c @@ -0,0 +1,57 @@ +/* Compute sexagesimal from decimal */ +#include +#include +#include +#include + +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] \n\nCompute sexagesimal from decimal input x.\n"); + printf("Options: -r gives hours instead of degrees\n"); + printf(" -d 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; +} diff --git a/jpg2fits.c b/jpg2fits.c new file mode 100644 index 0000000..1cc077e --- /dev/null +++ b/jpg2fits.c @@ -0,0 +1,289 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +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_scanlineifd[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=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; +} diff --git a/sex2dec.c b/sex2dec.c new file mode 100644 index 0000000..85420f4 --- /dev/null +++ b/sex2dec.c @@ -0,0 +1,42 @@ +// Compute decimal from sexagesimal +#include +#include +#include +#include + +// 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 \n -or- \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; +}