1
0
Fork 0

Major update

pull/2/head
Cees Bassa 2015-05-25 23:16:03 +02:00
parent c96e25f3e3
commit 7b095cbdac
6 changed files with 157 additions and 26 deletions

View File

@ -6,6 +6,7 @@
#include "cpgplot.h"
#include "cel.h"
#include <gsl/gsl_multifit.h>
#include <getopt.h>
#define LIM 256
#define NMAX 8192
@ -33,7 +34,7 @@ struct catalog {
int select[NMAX];
};
struct image read_fits(char *filename,int pnum);
struct catalog read_pixel_catalog(char *filename);
struct catalog read_pixel_catalog(char *filename,float,float,float);
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);
@ -42,10 +43,16 @@ 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);
double sex2dec(char *s);
void usage(void)
{
return;
}
int main(int argc,char *argv[])
{
int i;
int i,j;
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};
@ -56,27 +63,92 @@ int main(int argc,char *argv[])
float x,y,width;
char c,pixcat[LIM];
struct catalog cat,ast;
float sx=-10.0,sy=10.0,q;
float sx=-10.0,sy=10.0,q=0.0;
char *env,starfile[128];
float r,rmin=1.0,rmax=10.0,mmin=2.0,mmax=8.0,mag=9.0;
float r,srmin=1.0,srmax=10.0,smmin=2.0,smmax=8.0,mag=9.0,rmax=10.0,rmask=-1;
int arg=0;
char *fitsfile=NULL;
char sra[16],sde[16];
double ra,de;
FILE *file;
// Environment variables
env=getenv("ST_DATADIR");
sprintf(starfile,"%s/data/tycho2.dat",env);
// Read image
img=read_fits(argv[1],0);
// Decode options
while ((arg=getopt(argc,argv,"f:R:D:s:q:m:r:M:"))!=-1) {
switch(arg) {
// Hard coded
img.ra0=atof(argv[2]);
img.de0=atof(argv[3]);
q=atof(argv[4]);
case 'f':
fitsfile=optarg;
break;
case 'R':
strcpy(sra,optarg);
if (strchr(sra,':')!=NULL)
ra=15.0*sex2dec(sra);
else
ra=atof(sra);
break;
case 'D':
strcpy(sde,optarg);
if (strchr(sde,':')!=NULL)
de=sex2dec(sde);
else
de=atof(sde);
break;
case 's':
if (strchr(optarg,',')!=NULL) {
sscanf(optarg,"%f,%f",&sx,&sy);
} else {
sx=-atof(optarg);
sy=-sx;
}
break;
case 'q':
q=atof(optarg);
break;
case 'm':
mag=atof(optarg);
break;
case 'M':
rmask=atof(optarg);
break;
case 'r':
rmax=atof(optarg);
break;
case 'h':
usage();
return 0;
break;
default:
usage();
return 0;
}
}
// Read image
img=read_fits(fitsfile,0);
// Set variables
img.ra0=ra;
img.de0=de;
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);
sprintf(pixcat,"%s.cat",fitsfile);
cat=read_pixel_catalog(pixcat,img.x0,img.y0,rmask);
// Read astrometric catalog
ast=read_astrometric_catalog(starfile,mag,sx,sy,q);
@ -124,7 +196,7 @@ int main(int argc,char *argv[])
cpgsci(4);
for (i=0;i<ast.n;i++) {
r=rmax-(rmax-rmin)*(ast.mag[i]-mmin)/(mmax-mmin);
r=srmax-(srmax-srmin)*(ast.mag[i]-smmin)/(smmax-smmin);
// Upscale for image size
r*=img.naxis1/752.0;
@ -190,7 +262,7 @@ int main(int argc,char *argv[])
// Match catalogs
if (c=='m') {
nselect=match_catalogs(&cat,&ast,10.0);
nselect=match_catalogs(&cat,&ast,rmax);
redraw=1;
}
@ -221,6 +293,20 @@ int main(int argc,char *argv[])
redraw=1;
continue;
}
// Dump
if (c=='d') {
file=fopen("dump.dat","w");
for (i=0;i<nselect;i++) {
for (j=0;j<cat.n;j++)
if (cat.select[j]==i+1)
fprintf(file,"%lf %lf ",cat.x[j]-img.x0,cat.y[j]-img.y0);
for (j=0;j<ast.n;j++)
if (ast.select[j]==i+1)
fprintf(file,"%lf %lf\n",ast.ra[j],ast.de[j]);
}
fclose(file);
}
}
cpgend();
@ -286,12 +372,13 @@ struct image read_fits(char *filename,int pnum)
}
// Read pixel catalog
struct catalog read_pixel_catalog(char *filename)
struct catalog read_pixel_catalog(char *filename,float x0,float y0,float rmin)
{
int i=0;
FILE *file;
char line[LIM];
struct catalog c;
float dx,dy,r;
// Read catalog
file=fopen(filename,"r");
@ -305,6 +392,11 @@ struct catalog read_pixel_catalog(char *filename)
if (strstr(line,"#")!=NULL)
continue;
sscanf(line,"%f %f %f",&c.x[i],&c.y[i],&c.mag[i]);
dx=c.x[i]-x0;
dy=c.y[i]-y0;
r=sqrt(dx*dx+dy*dy);
if (r>rmin && rmin>0.0)
continue;
c.select[i]=0;
i++;
}
@ -576,3 +668,22 @@ int match_catalogs(struct catalog *cat,struct catalog *ast,float rmax)
printf("%d stars matched\n",n);
return n;
}
// 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;
}

View File

@ -46,17 +46,17 @@ void compute_position(double mjd,xyz_t satpos,struct site s,int satno,char *desi
// Compute positions
obspos_xyz(mjd,s.lng,s.lat,s.alt,&obspos,&obsvel);
// compute difference vector
dx=satpos.x-obspos.x;
dy=satpos.y-obspos.y;
dz=satpos.z-obspos.z;
// Celestial position
r=sqrt(dx*dx+dy*dy+dz*dz);
ra=modulo(atan2(dy,dx)*R2D,360.0);
de=asin(dz/r)*R2D;
// Precess position
if (precess_flag==1) {
precess(mjd,ra,de,mjd0,&ra0,&de0);
@ -176,7 +176,7 @@ int main(int argc,char *argv[])
// fprintf(stderr,"object %d not found in %s\n",p.satno,filename);
return;
}
// Initialize
imode=init_sgdp4(&orb);
if (imode==SGDP4_ERROR) {
@ -189,10 +189,12 @@ int main(int argc,char *argv[])
satpos_xyz(mjd+2400000.5,&satpos,&satvel);
compute_position(mjd,satpos,s,orb.satno,orb.desig,precess_flag);
} else {
file=fopen(fname,"r");
while (fgetline(file,line,LIM)>0) {
status=sscanf(line,"%lf",&mjd);
satpos_xyz(mjd+2400000.5,&satpos,&satvel);
strcpy(orb.desig,"14999A"); // FIX!
compute_position(mjd,satpos,s,orb.satno,orb.desig,precess_flag);
}
fclose(file);

View File

@ -49,6 +49,9 @@ propagate: propagate.o sgdp4.o satutl.o deep.o ferror.o
detect: detect.o
$(F77) -o detect detect.o -lm $(LFLAGS)
autodetect: autodetect.o
$(F77) -o autodetect autodetect.o -lm $(LFLAGS)
launchtle: launchtle.o sgdp4.o satutl.o deep.o ferror.o
$(CC) -o launchtle launchtle.o sgdp4.o satutl.o deep.o ferror.o -lm

View File

@ -816,6 +816,22 @@ int main(int argc,char *argv[])
obs.satno=trk.satno;
find_designation(obs.satno,obs.desig);
track(argv[1],obs,&img,frac);
x=frac*(trk.x1-trk.x0)+trk.x0;
y=frac*(trk.y1-trk.y0)+trk.y0;
width=100;
xmin=x-0.5*width;
xmax=x+0.5*width;
ymin=y-0.5*width*img.naxis2/img.naxis1;
ymax=y+0.5*width*img.naxis2/img.naxis1;
sigma=find_peak(img.zd,img.naxis1,img.naxis2,(int) xmin,(int) xmax,(int) ymin,(int) ymax,2.0,11,11,&x,&y);
printf("%f %f %f\n",x,y,sigma);
if (sigma>5.0) {
reduce_point(&obs,img,frac*img.exptime,x,y);
obs.x[0]=x;
obs.y[0]=y;
obs.state=2;
}
redraw=1;
layer=4;
}

View File

@ -113,7 +113,7 @@ int read_twoline(FILE *fp, long search_satno, orbit_t *orb)
// sscanf(st1+9,"%s",orb->desig);
strncpy(orb->desig,st1+9,8);
orb->desig[8]='\0';
return 0;
}

View File

@ -1627,7 +1627,7 @@ void skymap_plotsatellite(char *filename,int satno,double mjd0,double dt)
}
// Graves colouring
if (m.graves==1) {
if (m.graves!=0) {
if (s.illumg==1)
cpgsci(2);
else
@ -1841,7 +1841,7 @@ struct sat apparent_position(double mjd)
}
// Compute position at Graves
if (m.graves==1) {
if (m.graves!=0) {
graves_xyz(mjd,&grvpos,&grvvel);
dx=satpos.x-grvpos.x;
dy=satpos.y-grvpos.y;
@ -2375,7 +2375,7 @@ int plot_skymap(void)
plot_visibility(m.rvis);
// Graves visibility
if (m.graves==1)
if (m.graves==2)
plot_graves_visibility();
// Get time
@ -2573,9 +2573,8 @@ int plot_skymap(void)
// Toggle Graves illumination
if (c=='G') {
if (m.graves==0)
m.graves=1;
else if (m.graves==1)
m.graves++;
if (m.graves==3)
m.graves=0;
redraw=1;
}