spacecruft
/
sattools
1
0
Fork 0
Code Releases Activity

Various changes

pull/2/head
Cees Bassa 2013-07-10 10:54:14 +01:00
parent c812b5bb56
commit 59014a18d5
4 changed files with 327 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
reduce.c
View File

@ -301,8 +301,8 @@ void fit(struct observation *obs,struct image img)
// Fit x-pixel position
chi2y=linear_fit(dt,y,n,ay,say);
printf("x: %6.2f +- %.2f %7.3f +- %.3f %8.2f\n",ax[0],sax[0],ax[1],sax[1],chi2x);
printf("y: %6.2f +- %.2f %7.3f +- %.3f %8.2f\n",ay[0],say[0],ay[1],say[1],chi2y);
printf("x: %6.2f +- %.2f %7.3f +- %.3f %8.2f; %f pix/s\n",ax[0],sax[0],ax[1],sax[1],chi2x,ax[1]/img.nframes*img.exptime);
printf("y: %6.2f +- %.2f %7.3f +- %.3f %8.2f; %f pix/s\n",ay[0],say[0],ay[1],say[1],chi2y,ay[1]/img.nframes*img.exptime);
obs->x[0]=ax[0];
obs->y[0]=ay[0];
@ -721,6 +721,15 @@ int main(int argc,char *argv[])
iobject=0;
}
// Change fraction
if (c=='E') {
frac+=0.1;
if (frac>1.0)
frac=0.0;
printf("Fraction: %.1f\n",frac);
iobject=0;
}
// Reduce
if (c=='M' || c=='D') {
reduce_point(&obs,img,frac*img.exptime,x,y);

53
residuals.c
View File

@ -108,22 +108,54 @@ void compute_residual(char *filename,struct point p)
dr=sqrt(pow(dry*rx[0]-drx*ry[0],2)/(drx*drx+dry*dry));
if ((-rx[0]*drx-ry[0]*dry)<0.0)
dr*=-1;
printf("%s | %8.3f deg %8.3f sec %5.1f day, %.1f km\n",p.iod_line,dr,dt,age,r[0]);
printf("%s | %8.5f deg %8.3f sec %5.1f day, %.1f km\n",p.iod_line,dr,dt,age,r[0]);
return;
}
void split_file(struct data d,float dtmax)
{
int i,j,flag=0;
FILE *file;
char filename[LIM];
double mjd0,dt;
for (i=0,j=0;i<d.n;i++) {
if (flag==1) {
dt=86400*(d.p[i].mjd-mjd0);
if (dt>dtmax) {
if (file!=NULL)
fclose(file);
flag=0;
j++;
}
}
if (flag==0) {
mjd0=d.p[i].mjd;
flag=1;
sprintf(filename,"split%04d.dat",j+1);
file=fopen(filename,"w");
}
fprintf(file,"%s\n",d.p[i].iod_line);
}
if (file!=NULL)
fclose(file);
return;
}
int main(int argc,char *argv[])
{
int i,arg=0;
int i,arg=0,split=0;
struct data d;
char *datafile,catalog[LIM];
char *env;
float dt;
env=getenv("ST_TLEDIR");
sprintf(catalog,"%s/classfd.tle",env);
// Decode options
while ((arg=getopt(argc,argv,"d:c:h"))!=-1) {
while ((arg=getopt(argc,argv,"d:c:hs:"))!=-1) {
switch(arg) {
case 'd':
datafile=optarg;
@ -133,6 +165,11 @@ int main(int argc,char *argv[])
strcpy(catalog,optarg);
break;
case 's':
dt=atof(optarg);
split=1;
break;
case 'h':
usage();
return 0;
@ -146,8 +183,14 @@ int main(int argc,char *argv[])
// Read data
d=read_data(datafile);
for (i=0;i<d.n;i++)
compute_residual(catalog,d.p[i]);
if (split==1) {
split_file(d,dt);
} else {
for (i=0;i<d.n;i++)
compute_residual(catalog,d.p[i]);
}
return 0;
}

174
satfit.c
View File

@ -87,6 +87,74 @@ xyz_t get_position(double r0,int i0)
return pos;
}
int psearch(void)
{
int i,satno=99300;
double mjdmin,mjdmax;
int ia[7]={0,0,0,0,0,0,0};
double ecc,eccmin,eccmax,decc;
double rev,revmin,revmax,drev;
double argp,argpmin,argpmax,dargp;
char line1[70],line2[70];
FILE *file;
// Provide
printf("Mean motion [min, max, stepsize]: \n");
scanf("%lf %lf %lf",&revmin,&revmax,&drev);
printf("Eccentricity [min, max, stepsize]: \n");
scanf("%lf %lf %lf",&eccmin,&eccmax,&decc);
// printf("Argument of perigee [min, max, stepsize]: \n");
// scanf("%lf %lf %lf",&argpmin,&argpmax,&dargp);
// Step 1: select all points
// for (i=0;i<d.n;i++)
// d.p[i].flag=2;
// Step 2: get time range
time_range(&mjdmin,&mjdmax,2);
file=fopen("search.dat","w");
// Step 4: Loop over eccentricity
for (rev=revmin;rev<revmax;rev+=drev) {
for (ecc=eccmin;ecc<eccmax;ecc+=decc) {
// for (argp=argpmin;argp<argpmax;argp+=dargp) {
orb.satno=satno;
orb.ecc=ecc;
orb.rev=rev;
//orb.argp=argp*D2R;
// Set parameters
for (i=0;i<7;i++)
ia[i]=0;
// Step 4: loop over parameters
for (i=0;i<5;i++) {
if (i==0) ia[4]=1;
if (i==1) ia[1]=1;
if (i==2) ia[0]=1;
// if (i==3) ia[5]=1;
if (i==4) ia[3]=1;
// Do fit
fit(orb,ia);
}
fit(orb,ia);
fit(orb,ia);
fit(orb,ia);
printf("%8.5lf %8.6lf %8.3lf %8.3lf %8.3lf %8.3lf %8.5lf\n",orb.rev,orb.ecc,orb.argp*R2D,orb.ascn*R2D,orb.mnan*R2D,orb.eqinc*R2D,d.rms);
fprintf(file,"%8.5lf %8.6lf %8.3lf %8.3lf %8.3lf %8.3lf %8.5lf\n",orb.rev,orb.ecc,orb.argp*R2D,orb.ascn*R2D,orb.mnan*R2D,orb.eqinc*R2D,d.rms);
}
fprintf(file,"\n");
// }
}
fclose(file);
return orb.satno;
}
int circular_fit(void)
{
int i;
@ -202,12 +270,12 @@ int main(int argc,char *argv[])
char filename[64];
int satno=-1;
double mjdmin,mjdmax;
int arg=0;
char *datafile,*catalog;
int arg=0,elset=0,circular=0,tleout=0,noplot=0;
char *datafile,*catalog,tlefile[LIM];
orbit_t orb0;
// Decode options
while ((arg=getopt(argc,argv,"d:c:i:ha"))!=-1) {
while ((arg=getopt(argc,argv,"d:c:i:haCo:p"))!=-1) {
switch(arg) {
case 'd':
datafile=optarg;
@ -221,6 +289,19 @@ int main(int argc,char *argv[])
satno=atoi(optarg);
break;
case 'C':
circular=1;
break;
case 'p':
noplot=1;
break;
case 'o':
tleout=1;
strcpy(tlefile,optarg);
break;
case 'h':
usage();
return 0;
@ -247,6 +328,19 @@ int main(int argc,char *argv[])
freopen("/tmp/stderr.txt","w",stderr);
// Fit circular orbit
if (circular==1) {
for (i=0;i<d.n;i++)
d.p[i].flag=2;
satno=circular_fit();
plot_residuals=1;
quit=1;
// Dump tle
if (tleout==1)
print_tle(orb,tlefile);
}
// Adjust
if (adjust==1) {
orb0=orb;
@ -256,6 +350,17 @@ int main(int argc,char *argv[])
plot_residuals=1;
redraw=1;
quit=1;
// Dump tle
if (tleout==1)
print_tle(orb,tlefile);
}
// Exit before plotting
if (quit==1 && noplot==1) {
free(d.p);
fclose(stderr);
return 0;
}
cpgopen("/xs");
@ -408,6 +513,14 @@ int main(int argc,char *argv[])
redraw=1;
}
// Search
if (c=='S') {
satno=psearch();
plot_residuals=1;
ia[0]=ia[1]=ia[4]=ia[5]=1;
redraw=1;
}
// Change
if (c=='c') {
printf("(1) Inclination, (2) Ascending Node, (3) Eccentricity,\n(4) Arg. of Perigee, (5) Mean Anomaly, (6) Mean Motion,\n(7) B* drag, (8) Epoch, (9) Satellite ID\n\nWhich parameter to change: ");
@ -478,15 +591,48 @@ int main(int argc,char *argv[])
// Default tle
if (c=='t') {
orb.satno=99999;
orb.eqinc=0.5*M_PI;
orb.ascn=0.0;
orb.ecc=0.0;
orb.argp=0.0;
orb.mnan=0.0;
orb.rev=14.0;
orb.bstar=0.5e-4;
orb.ep_day=mjd2doy(0.5*(mjdmin+mjdmax),&orb.ep_year);
satno=99999;
if (elset==0) {
orb.eqinc=0.5*M_PI;
orb.ascn=0.0;
orb.ecc=0.0;
orb.argp=0.0;
orb.mnan=0.0;
orb.rev=14.0;
orb.bstar=0.5e-4;
printf("LEO orbit\n");
} else if (elset==1) {
orb.eqinc=20.0*D2R;
orb.ascn=0.0;
orb.ecc=0.7;
orb.argp=0.0;
orb.mnan=0.0;
orb.rev=2.25;
orb.bstar=0.0;
printf("GTO orbit\n");
} else if (elset==2) {
orb.eqinc=10.0*D2R;
orb.ascn=0.0;
orb.ecc=0.0;
orb.argp=0.0;
orb.mnan=0.0;
orb.rev=1.0027;
orb.bstar=0.0;
printf("GSO orbit\n");
} else if (elset==3) {
orb.eqinc=63.4*D2R;
orb.ascn=0.0;
orb.ecc=0.7;
orb.argp=0.0;
orb.mnan=0.0;
orb.rev=2.0;
orb.bstar=0.0;
printf("HEO orbit\n");
}
elset++;
if (elset>3)
elset=0;
print_orb(&orb);
printf("\n================================================================================\n");
click=0;
@ -901,6 +1047,10 @@ double chisq(double a[])
} else if (a[0]>180.0) {
a[0]=180.0;
}
if (a[5]>17.0)
a[5]=17.0;
if (a[5]<0.1)
a[5]=0.1;
// Set parameters
@ -1151,10 +1301,14 @@ void fit(orbit_t orb,int *ia)
double db[7]={5.0,5.0,0.1,5.0,5.0,0.5,0.0001};
a[0]=orb.eqinc*R2D;
da[0]=da[0]*R2D;
a[1]=orb.ascn*R2D;
da[1]=da[1]*R2D;
a[2]=orb.ecc;
a[3]=orb.argp*R2D;
da[3]=da[3]*R2D;
a[4]=orb.mnan*R2D;
da[4]=da[4]*R2D;
a[5]=orb.rev;
a[6]=orb.bstar;

108
skymap.c
View File

@ -28,14 +28,14 @@ struct map {
float length;
float minmag,maxmag,minrad,maxrad;
char orientation[LIM],projection[4],observer[32];
char nfd[LIM],starfile[LIM],tlefile[LIM],iodfile[LIM];
char nfd[LIM],starfile[LIM],tlefile[LIM],iodfile[LIM],xyzstring[LIM];
char datadir[LIM],tledir[LIM];
double lat,lng;
double h,sra,sde,sazi,salt;
float alt,timezone;
float fw,fh;
int level,grid,site_id;
int leoflag,iodflag,iodpoint,visflag,planar,pssatno,psnr;
int leoflag,iodflag,iodpoint,visflag,planar,pssatno,psnr,xyzflag;
float psrmin,psrmax,rvis;
} m;
struct sat {
@ -102,7 +102,7 @@ void get_site(int site_id);
void usage()
{
printf("skymap t:c:i:R:D:hs:d:l:P:r:V:\n\n");
printf("skymap t:c:i:R:D:hs:d:l:P:r:V:p:\n\n");
printf("t date/time (yyyy-mm-ddThh:mm:ss.sss) [default: now]\n");
printf("c TLE catalog file [default: classfd.tle]\n");
printf("i satellite ID (NORAD) [default: all]\n");
@ -115,6 +115,7 @@ void usage()
printf("P planar search satellite ID\n");
printf("r planar search altitude\n");
printf("V altitude for visibility contours\n");
printf("p xyz position of object\n");
return;
}
@ -148,6 +149,7 @@ void init_skymap(void)
m.mjd=-1.0;
m.leoflag=1;
m.iodflag=0;
m.xyzflag=0;
m.visflag=0;
m.planar=0;
@ -260,6 +262,95 @@ void read_iod(char *filename,int iobs)
return;
}
void plot_xyz(double mjd,char *string)
{
struct sat s;
double jd,rsun,rearth,rsat;
double dx,dy,dz,dvx,dvy,dvz;
xyz_t satpos,obspos,obsvel,satvel,sunpos;
double sra,sde;
sscanf(string,"%lf,%lf,%lf",&satpos.x,&satpos.y,&satpos.z);
// Get positions
obspos_xyz(mjd,&obspos,&obsvel);
sunpos_xyz(mjd,&sunpos,&sra,&sde);
// Age
s.age=0.0;
// Sat positions
s.x=satpos.x;
s.y=satpos.y;
s.z=satpos.z;
s.vx=satvel.x;
s.vy=satvel.y;
s.vz=satvel.z;
// Sun position from satellite
dx=-satpos.x+sunpos.x;
dy=-satpos.y+sunpos.y;
dz=-satpos.z+sunpos.z;
// Distances
rsun=sqrt(dx*dx+dy*dy+dz*dz);
rearth=sqrt(satpos.x*satpos.x+satpos.y*satpos.y+satpos.z*satpos.z);
s.h=rearth-XKMPER;
// Angles
s.psun=asin(696.0e3/rsun)*R2D;
s.pearth=asin(6378.135/rearth)*R2D;
s.p=acos((-dx*satpos.x-dy*satpos.y-dz*satpos.z)/(rsun*rearth))*R2D;
// Visibility state
if (s.p-s.pearth<-s.psun)
strcpy(s.state,"eclipsed");
else if (s.p-s.pearth>-s.psun && s.p-s.pearth<s.psun)
strcpy(s.state,"umbra");
else if (s.p-s.pearth>s.psun)
strcpy(s.state,"sunlit");
// Position differences
dx=satpos.x-obspos.x;
dy=satpos.y-obspos.y;
dz=satpos.z-obspos.z;
dvx=satvel.x-obsvel.x;
dvy=satvel.y-obsvel.y;
dvz=satvel.z-obsvel.z;
// Celestial position
s.r=sqrt(dx*dx+dy*dy+dz*dz);
s.v=(dvx*dx+dvy*dy+dvz*dz)/s.r;
s.ra=modulo(atan2(dy,dx)*R2D,360.0);
s.de=asin(dz/s.r)*R2D;
// Phase
s.phase=acos(((obspos.x-satpos.x)*(sunpos.x-satpos.x)+(obspos.y-satpos.y)*(sunpos.y-satpos.y)+(obspos.z-satpos.z)*(sunpos.z-satpos.z))/(rsun*s.r))*R2D;
// Magnitude
if (strcmp(s.state,"sunlit")==0)
s.mag=STDMAG-15.0+5*log10(s.r)-2.5*log10(sin(s.phase*D2R)+(M_PI-s.phase*D2R)*cos(s.phase*D2R));
else
s.mag=15;
// Convert and project
if (strcmp(m.orientation,"horizontal")==0) {
equatorial2horizontal(mjd,s.ra,s.de,&s.azi,&s.alt);
forward(s.azi,s.alt,&s.rx,&s.ry);
} else if (strcmp(m.orientation,"equatorial")==0) {
forward(s.ra,s.de,&s.rx,&s.ry);
}
cpgsci(3);
cpgpt1(s.rx,s.ry,17);
cpgsch(0.6);
cpgtext(s.rx,s.ry," xyz");
cpgsch(1.0);
cpgsci(1);
printf("%s %6.3f %6.3f %.1f km\n",m.nfd,s.ra,s.de,s.r);
return;
}
int main(int argc,char *argv[])
{
@ -271,7 +362,7 @@ int main(int argc,char *argv[])
init_skymap();
// Decode options
while ((arg=getopt(argc,argv,"t:c:i:R:D:hs:d:l:P:r:V:"))!=-1) {
while ((arg=getopt(argc,argv,"t:c:i:R:D:hs:d:l:P:r:V:p:"))!=-1) {
switch(arg) {
case 't':
@ -312,6 +403,11 @@ int main(int argc,char *argv[])
m.psnr=8;
break;
case 'p':
strcpy(m.xyzstring,optarg);
m.xyzflag=1;
break;
case 'r':
m.psrmin=atof(optarg);
m.psrmax=atof(optarg);
@ -1812,6 +1908,10 @@ int plot_skymap(void)
if (m.iodflag==1)
plot_iod(m.iodfile);
// Plot XYZ position
if (m.xyzflag==1)
plot_xyz(m.mjd,m.xyzstring);
// Plot visibility
if (m.visflag==1 && strcmp(m.orientation,"horizontal")==0)
plot_visibility(m.rvis);