1
0
Fork 0
pull/3/head
Cees Bassa 2016-07-30 17:20:46 +02:00
commit b08c14b2e6
16 changed files with 1998 additions and 701 deletions

View File

@ -1,12 +1,28 @@
Satellite Tracking Toolkit
=========
Required libraries:
qfits-5.2.0: ftp://ftp.eso.org/pub/qfits/qfits-5.2.0.tar.gz
pgplot-5.2.2: http://www.astro.caltech.edu/~tjp/pgplot/
gsl-1.15: ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz
wcslib-2.9
Sattools is a collection of tools to facilitate Photographic and Video satellite tracking.
Required environment variables:
ST_COSPAR: COSPAR number
ST_DATADIR: path to sattools directory
ST_TLEDIR: path to TLE directory
Build
------
* Clone locally the code repository
* Install common dependencies
* gfortran
* gcc
* libpng-dev
* libx11-dev
* libjpeg-dev
* libexif-dev
* Build & install required libraries
* qfits-5.2.0: ftp://ftp.eso.org/pub/qfits/qfits-5.2.0.tar.gz
* pgplot-5.2.2: http://www.astro.caltech.edu/~tjp/pgplot/
* gsl-1.15: ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz
* wcslib-2.9: http://www.epta.eu.org/~bassa/wcslib-2.9.tar
* Run `make` on the sattools folder
Run
-----
You will need to set the following environment variables to run sattools:
* `ST_COSPAR` COSPAR number
* `ST_DATADIR` path to sattools directory
* `ST_TLEDIR` path to TLE directory

View File

@ -86,63 +86,68 @@ int main(int argc,char *argv[])
sprintf(starfile,"%s/data/tycho2.dat",env);
// Decode options
while ((arg=getopt(argc,argv,"f:R:D:s:q:m:r:M:"))!=-1) {
switch(arg) {
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;
if (argc>1) {
while ((arg=getopt(argc,argv,"f:R:D:s:q:m:r:M:"))!=-1) {
switch(arg) {
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;
}
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;
}
}
} else {
usage();
return 0;
}
// Read image
@ -224,6 +229,22 @@ int main(int argc,char *argv[])
if (c=='q')
break;
// Help
if (c=='h') {
printf("Calibrates astrometry. Initially requires manual matching of at least three stars. Use 'a' to select star on the image, then 'b' to select star from the catalog (blue circles). If at least three stars are matched, use 'f' and 'm' to converge the calibration.\n\n");
printf("q Quit\n");
printf("a Select star on image\n");
printf("b Select star from catalog\n");
printf("c Center image on pixel\n");
printf("f Fit calibration\n");
printf("m Match stars using current calibration\n");
printf("z/+ Zoom in on cursor\n");
printf("x/- Zoom out on cursor\n");
printf("p Plot sextractor catalog\n");
printf("r Reset zoom\n");
}
// Select pixel catalog
if (c=='a' && click==0) {
i=select_nearest(cat,x,y);
@ -259,7 +280,7 @@ int main(int argc,char *argv[])
}
// Zoom
if (c=='z' || c=='+') {
if (c=='z' || c=='+' || c=='=') {
width/=1.25;
xmin=x-0.5*width;
xmax=x+0.5*width;
@ -276,7 +297,7 @@ int main(int argc,char *argv[])
}
// Unzoom
if (c=='x' || c=='+' || c=='=') {
if (c=='x' || c=='-') {
width*=1.25;
xmin=x-0.5*width;
xmax=x+0.5*width;

File diff suppressed because it is too large Load Diff

View File

@ -15,6 +15,7 @@
0434 IR -26.1030 27.9288 1646 Ian Roberts
0710 LS 52.3261 10.6756 85 Lutz Schindler
1056 MK 57.0122 23.9833 4 Martins Keruss
1086 NK 46.4778 30.7556 56 Nikolay Koshkin
1244 AM 44.3932 33.9701 69 Andriy Makeyev
1747 DD 45.7275 -72.3526 191 Daniel Deak
1775 KF 44.6062 -75.6910 200 Kevin Fetter
@ -33,6 +34,7 @@
4541 AR 41.9639 12.4531 80 Alberto Rango
4542 AR 41.9683 12.4545 80 Alberto Rango
4641 AR 41.1060 16.9010 70 Alberto Rango
5555 SG 56.1019 94.5533 154 Sergey Guryanov
5918 BG 59.2985 18.1045 40 Bjorn Gimle
5919 BG 59.2615 18.6206 30 Bjorn Gimle
6226 SC 28.4861 -97.8194 110 Scott Campbell
@ -45,14 +47,17 @@
8536 TL 36.8479 -76.5010 4 Tim Luton
8539 SN 39.4707 -79.3388 839 Steve Newcomb
8597 TB -34.9638 138.6333 100 Tony Beresford
8600 PC -32.9770 151.6477 18 Paul Camilleri
8730 EC 30.3086 -97.7279 150 Ed Cannon
8739 DB 37.1133 -121.7028 282 Derek Breit
9461 KO 47.9175 19.89365 938 Konkoly Obs
9633 JN 33.0206 -96.9982 153 Jim Nix
9730 MM 30.3150 -97.8660 280 Mike McCants
6242 JM 42.9453 -2.8284 623 Jon Mikel
6241 JM 42.9565 -2.8203 619 Jon Mikel
4160 BD 51.2793 5.4768 35 Bram Dorreman
9001 MS -32.3793 20.810682 1760 Master K95
9000 PP 40.0848 21.1581 1433 Pierros Papadeas
9001 AF 44.5700 6.68 1850 Alain Figer
9002 AO 62.2544 26.5962 150 Arto Oksanen
9999 GR 47.348 5.5151 100 Graves
9000 NJ 51.4198 5.4086 25 Nico Janssen
9905 DR -26.8724 26.6481 1200 Deon van Rooyen

View File

@ -684,13 +684,14 @@ int main(int argc,char *argv[])
printf("%d points above %g sigma\n",n,sigma);
// Exit if too many
/*
if (n>2500) {
file=fopen("skipped.dat","a");
fprintf(file,"%s : %d points above %g sigma\n",fitsfile,n,sigma);
fclose(file);
return 0;
}
*/
// Fill points
p=(struct point *) malloc(sizeof(struct point)*n);
for (i=0,l=0;i<ff.naxis1;i++) {

View File

@ -13,6 +13,12 @@ F77 = gfortran
all:
make addwcs angular calibrate dec2sex faketle fitsheader fitskey imgstat jpg2fits jpgstack measure pgm2fits plotfits pstrack rde2iod reduce residuals runsched satfit satid satmap satorbit sex2dec skymap tle2ole tleinfo uk2iod viewer wcsfit deproject slewto waitfor pass detect launchtle propagate fakeiod csv2tle normal posmatch posvel xyz2tle mvtle
selectiod: selectiod.o
$(CC) -o selectiod selectiod.o -lm
planscan: planscan.o sgdp4.o satutl.o deep.o ferror.o
$(CC) -o planscan planscan.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS)
tle2rv: tle2rv.o sgdp4.o satutl.o deep.o ferror.o
$(CC) -o tle2rv tle2rv.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS)

View File

@ -562,7 +562,7 @@ int main(int argc,char *argv[])
cpgmtxt("T",6.0,0.0,0.0,text);
sprintf(text,"R.A.: %10.5f (%4.1f'') Decl.: %10.5f (%4.1f'')",img[iimg].ra0,img[iimg].xrms,img[iimg].de0,img[iimg].yrms);
cpgmtxt("T",4.8,0.0,0.0,text);
sprintf(text,"FoV: %.2f\\(2218)x%.2f\\(2218) Scale: %.2f''x%.2f'' pix\\u-1\\d",img[iimg].naxis1*sqrt(img[iimg].a[1]*img[iimg].a[1]+img[iimg].b[1]*img[iimg].b[1])/3600.0,img[iimg].naxis2*sqrt(img[iimg].a[2]*img[iimg].a[2]+img[iimg].b[2]*img[iimg].b[2])/3600.0,sqrt(img[iimg].a[1]*img[iimg].a[1]+img[iimg].b[1]*img[iimg].b[1]),sqrt(img[iimg].a[2]*img[iimg].a[2]+img[iimg].b[2]*img[iimg].b[2]));
sprintf(text,"FoV: %.2f\\(2218)x%.2f\\(2218) Scale: %.2f''x%.2f'' pix\\u-1\\d Fraction: %.1f",img[iimg].naxis1*sqrt(img[iimg].a[1]*img[iimg].a[1]+img[iimg].b[1]*img[iimg].b[1])/3600.0,img[iimg].naxis2*sqrt(img[iimg].a[2]*img[iimg].a[2]+img[iimg].b[2]*img[iimg].b[2])/3600.0,sqrt(img[iimg].a[1]*img[iimg].a[1]+img[iimg].b[1]*img[iimg].b[1]),sqrt(img[iimg].a[2]*img[iimg].a[2]+img[iimg].b[2]*img[iimg].b[2]),frac);
cpgmtxt("T",3.6,0.0,0.0,text);
zmin=img[iimg].zmin;
@ -632,6 +632,7 @@ int main(int argc,char *argv[])
else if (frac<0.49)
frac=0.5;
printf("Fraction: %.1f\n",frac);
redraw=1;
}
if (c=='p' || c=='X') {
@ -822,6 +823,27 @@ int main(int argc,char *argv[])
redraw=1;
continue;
}
// Help
if (c=='h') {
printf("q Quit\n");
printf("e Change of exposure (0.0 = start, 0.5 = middle, 1.0 = end)\n");
printf("p (right) Plot objects\n");
printf("l Plot star catalog\n");
printf("d Set object NORAD designation\n");
printf("c Center on cursor\n");
printf("0-9 Set zoom level\n");
printf("s/] Cycle through images forward\n");
printf("a/[ Cycle through images backward\n");
printf("o Show maximum pixels of all images\n");
printf("TAB Cycle through frame at current zoom level\n");
printf("z/+ Zoom in on cursor\n");
printf("x/- Zoom out on cursor\n");
printf("r Reset zoom\n");
printf("R Reset setup\n");
printf("w Write IOD observation\n");
printf("m (mid) measure position\n");
}
}
cpgend();

579
planscan.c 100644
View File

@ -0,0 +1,579 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <getopt.h>
#include <time.h>
#include "sgdp4h.h"
#define LIM 128
#define D2R M_PI/180.0
#define R2D 180.0/M_PI
#define XKMPER 6378.135 // Earth radius in km
#define XKMPAU 149597879.691 // AU in km
#define FLAT (1.0/298.257)
#define STDMAG 6.0
#define MMAX 1024
long Isat=0;
long Isatsel=0;
extern double SGDP4_jd0;
struct map {
double lat,lng;
float alt;
char observer[32];
int site_id;
} m;
void get_site(int site_id);
double nfd2mjd(char *date);
double date2mjd(int year,int month,double day);
void nfd_now(char *s);
void obspos_xyz(double mjd,xyz_t *pos,xyz_t *vel);
void sunpos_xyz(double mjd,xyz_t *pos,double *ra,double *de);
double gmst(double mjd);
double dgmst(double mjd);
double modulo(double x,double y);
void equatorial2horizontal(double mjd,double ra,double de,double *azi,double *alt);
void mjd2date(double mjd,char *date);
int properties(kep_t K,xyz_t obspos,xyz_t sunpos,float radius,float t,double *ra,double *de,double *r,float *mag);
void dec2sex(double x,char *s,int f,int len);
void usage(void)
{
return;
}
int main(int argc,char *argv[])
{
int arg=0,imode,i;
int satno=-1;
orbit_t orb;
FILE *fp;
char tlefile[256]="classfd.tle";
char nfd[32];
double mjd0,mjd,jd,tsince;
int rv,withvel;
kep_t K;
double radius=1080.0;
xyz_t obspos,obsvel,sunpos;
double p,pmin,p0,p1,r,ra,de,azi,alt,sazi,salt,altmin=10.0,saltmin=-6.0;
float mag,mmin;
int state,pstate,nstate;
float t,length=86400.0,dt=60.0;
char sra[16],sde[16],type[32];
// Initialize
nfd_now(nfd);
mjd0=nfd2mjd(nfd);
get_site(4171);
// Decode options
while ((arg=getopt(argc,argv,"t:c:i:s:l:hS:A:r:"))!=-1) {
switch (arg) {
case 't':
strcpy(nfd,optarg);
mjd0=nfd2mjd(nfd);
break;
case 'c':
strcpy(tlefile,optarg);
break;
case 's':
get_site(atoi(optarg));
break;
case 'i':
satno=atoi(optarg);
break;
case 'r':
radius=atof(optarg);
break;
case 'l':
length=atoi(optarg);
if (strchr(optarg,'h')!=NULL)
length*=3600;
else if (strchr(optarg,'m')!=NULL)
length*=60;
else if (strchr(optarg,'d')!=NULL)
length*=86400;
break;
case 'S':
saltmin=atof(optarg);
break;
case 'A':
altmin=atof(optarg);
break;
case 'h':
usage();
return 0;
break;
default:
usage();
return 0;
}
}
// Error checking
if (satno<=0) {
fprintf(stderr,"ERROR: NORAD satellite number not provided!\n");
return -1;
}
// Get TLE
fp=fopen(tlefile,"rb");
if (fp==NULL) {
fprintf(stderr,"ERROR: Failed to open file with TLEs (%s)!\n",tlefile);
return -1;
}
// Read TLE
while (read_twoline(fp,satno,&orb)==0) {
Isat=orb.satno;
imode=init_sgdp4(&orb);
if (imode==SGDP4_ERROR)
continue;
}
fclose(fp);
// Object found?
if (orb.satno!=satno) {
fprintf(stderr,"ERROR: Object %d not found in %s!\n",satno,tlefile);
return -1;
}
// Object found?
if (orb.rev<10.0) {
fprintf(stderr,"ERROR: Object %d is not in a LEO orbit.\n",satno);
return -1;
}
// Print header
printf("Observer: %s (%04d) [%+.4f, %+.4f, %.0fm]\n",m.observer,m.site_id,m.lat,m.lng,m.alt*1000.0);
printf("Elements: %s\n",tlefile);
printf("Object: %d\n",satno);
printf("Start UT Date/Time: %s for %g h \n\n",nfd,length/3600.0);
printf("UT Date/Time R.A. Decl. Azi. Alt. Range Mag Sun Alt. Type\n");
printf(" (deg) (deg) (km) (deg)\n");
printf("=====================================================================================\n");
for (t=0.0;t<length;t+=dt) {
// (Modified) Julian Date
mjd=mjd0+t/86400.0;
jd=mjd+2400000.5;
// Get kepler
tsince=1440.0*(jd-SGDP4_jd0);
rv=sgdp4(tsince,1,&K);
// Get positions
obspos_xyz(mjd,&obspos,&obsvel);
sunpos_xyz(mjd,&sunpos,&ra,&de);
equatorial2horizontal(mjd,ra,de,&sazi,&salt);
p0=0.0;
p1=2.0*M_PI;
for (i=0,pmin=0.0,mmin=15.0;i<MMAX;i++) {
p=p0+(p1-p0)*(float) i/(float) (MMAX-1);
state=properties(K,obspos,sunpos,radius,p,&ra,&de,&r,&mag);
if (mag<mmin) {
pmin=p;
mmin=mag;
}
}
p0=pmin-4.0*M_PI/(float) MMAX;
p1=pmin+4.0*M_PI/(float) MMAX;
for (i=0,pmin=0.0,mmin=15.0;i<MMAX;i++) {
p=p0+(p1-p0)*(float) i/(float) (MMAX-1);
state=properties(K,obspos,sunpos,radius,p,&ra,&de,&r,&mag);
if (mag<mmin) {
pmin=p;
mmin=mag;
}
}
// Get properties before and after maximum
pstate=properties(K,obspos,sunpos,radius,pmin-0.01,&ra,&de,&r,&mag);
nstate=properties(K,obspos,sunpos,radius,pmin+0.01,&ra,&de,&r,&mag);
state=properties(K,obspos,sunpos,radius,pmin,&ra,&de,&r,&mag);
if (pstate<state && state==nstate)
strcpy(type,"Egress");
else if (pstate==state && state>nstate)
strcpy(type,"Ingress");
else
strcpy(type,"Maximum");
ra=modulo(ra,360.0);
equatorial2horizontal(mjd,ra,de,&azi,&alt);
azi=modulo(azi+180.0,360.0);
mjd2date(mjd,nfd);
dec2sex(ra/15.0,sra,0,2);
dec2sex(de,sde,0,2);
if (alt>altmin && salt<saltmin)
printf("%s %s %s %6.2f %6.2f %7.1f %5.2f %6.2f %s\n",nfd,sra,sde,azi,alt,r,mag,salt,type);
}
return 0;
}
// Get observing site
void get_site(int site_id)
{
int i=0;
char line[LIM];
FILE *file;
int id;
double lat,lng;
float alt;
char abbrev[3],observer[64],filename[LIM],*env;
env=getenv("ST_DATADIR");
sprintf(filename,"%s/data/sites.txt",env);
file=fopen(filename,"r");
if (file==NULL) {
fprintf(stderr,"File with site information not found (%s)!\n",filename);
return;
}
while (fgets(line,LIM,file)!=NULL) {
// Skip
if (strstr(line,"#")!=NULL)
continue;
// Strip newline
line[strlen(line)-1]='\0';
// Read data
sscanf(line,"%4d %2s %lf %lf %f",
&id,abbrev,&lat,&lng,&alt);
strcpy(observer,line+38);
// Change to km
alt/=1000.0;
if (id==site_id) {
m.lat=lat;
m.lng=lng;
m.alt=alt;
m.site_id=id;
strcpy(m.observer,observer);
}
}
fclose(file);
return;
}
// 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;
}
// 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==1582 && 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;
}
// Present nfd
void nfd_now(char *s)
{
time_t rawtime;
struct tm *ptm;
// Get UTC time
time(&rawtime);
ptm=gmtime(&rawtime);
sprintf(s,"%04d-%02d-%02dT%02d:%02d:%02d",ptm->tm_year+1900,ptm->tm_mon+1,ptm->tm_mday,ptm->tm_hour,ptm->tm_min,ptm->tm_sec);
return;
}
// Observer position
void obspos_xyz(double mjd,xyz_t *pos,xyz_t *vel)
{
double ff,gc,gs,theta,s,dtheta;
s=sin(m.lat*D2R);
ff=sqrt(1.0-FLAT*(2.0-FLAT)*s*s);
gc=1.0/ff+m.alt/XKMPER;
gs=(1.0-FLAT)*(1.0-FLAT)/ff+m.alt/XKMPER;
theta=gmst(mjd)+m.lng;
dtheta=dgmst(mjd)*D2R/86400;
pos->x=gc*cos(m.lat*D2R)*cos(theta*D2R)*XKMPER;
pos->y=gc*cos(m.lat*D2R)*sin(theta*D2R)*XKMPER;
pos->z=gs*sin(m.lat*D2R)*XKMPER;
vel->x=-gc*cos(m.lat*D2R)*sin(theta*D2R)*XKMPER*dtheta;
vel->y=gc*cos(m.lat*D2R)*cos(theta*D2R)*XKMPER*dtheta;
vel->z=0.0;
return;
}
// Solar position
void sunpos_xyz(double mjd,xyz_t *pos,double *ra,double *de)
{
double jd,t,l0,m,e,c,r;
double n,s,ecl;
jd=mjd+2400000.5;
t=(jd-2451545.0)/36525.0;
l0=modulo(280.46646+t*(36000.76983+t*0.0003032),360.0)*D2R;
m=modulo(357.52911+t*(35999.05029-t*0.0001537),360.0)*D2R;
e=0.016708634+t*(-0.000042037-t*0.0000001267);
c=(1.914602+t*(-0.004817-t*0.000014))*sin(m)*D2R;
c+=(0.019993-0.000101*t)*sin(2.0*m)*D2R;
c+=0.000289*sin(3.0*m)*D2R;
r=1.000001018*(1.0-e*e)/(1.0+e*cos(m+c));
n=modulo(125.04-1934.136*t,360.0)*D2R;
s=l0+c+(-0.00569-0.00478*sin(n))*D2R;
ecl=(23.43929111+(-46.8150*t-0.00059*t*t+0.001813*t*t*t)/3600.0+0.00256*cos(n))*D2R;
*ra=atan2(cos(ecl)*sin(s),cos(s))*R2D;
*de=asin(sin(ecl)*sin(s))*R2D;
pos->x=r*cos(*de*D2R)*cos(*ra*D2R)*XKMPAU;
pos->y=r*cos(*de*D2R)*sin(*ra*D2R)*XKMPAU;
pos->z=r*sin(*de*D2R)*XKMPAU;
return;
}
// Greenwich Mean Sidereal Time
double dgmst(double mjd)
{
double t,dgmst;
t=(mjd-51544.5)/36525.0;
dgmst=360.98564736629+t*(0.000387933-t/38710000);
return dgmst;
}
// Greenwich Mean Sidereal Time
double gmst(double mjd)
{
double t,gmst;
t=(mjd-51544.5)/36525.0;
gmst=modulo(280.46061837+360.98564736629*(mjd-51544.5)+t*t*(0.000387933-t/38710000),360.0);
return gmst;
}
// Return x modulo y [0,y)
double modulo(double x,double y)
{
x=fmod(x,y);
if (x<0.0) x+=y;
return x;
}
// Convert equatorial into horizontal coordinates
void equatorial2horizontal(double mjd,double ra,double de,double *azi,double *alt)
{
double h;
h=gmst(mjd)+m.lng-ra;
*azi=modulo(atan2(sin(h*D2R),cos(h*D2R)*sin(m.lat*D2R)-tan(de*D2R)*cos(m.lat*D2R))*R2D,360.0);
*alt=asin(sin(m.lat*D2R)*sin(de*D2R)+cos(m.lat*D2R)*cos(de*D2R)*cos(h*D2R))*R2D;
return;
}
// Compute Date from Julian Day
void mjd2date(double mjd,char *date)
{
double f,jd,dday;
int z,alpha,a,b,c,d,e;
int year,month,day,hour,min;
float sec,x;
jd=mjd+2400000.5;
jd+=0.5;
z=floor(jd);
f=fmod(jd,1.);
if (z<2299161)
a=z;
else {
alpha=floor((z-1867216.25)/36524.25);
a=z+1+alpha-floor(alpha/4.);
}
b=a+1524;
c=floor((b-122.1)/365.25);
d=floor(365.25*c);
e=floor((b-d)/30.6001);
dday=b-d-floor(30.6001*e)+f;
if (e<14)
month=e-1;
else
month=e-13;
if (month>2)
year=c-4716;
else
year=c-4715;
day=(int) floor(dday);
x=24.0*(dday-day);
x=3600.*fabs(x)+0.0001;
sec=fmod(x,60.);
x=(x-sec)/60.;
min=fmod(x,60.);
x=(x-min)/60.;
hour=x;
sec=floor(1000.0*sec)/1000.0;
sprintf(date,"%04d-%02d-%02dT%02d:%02d:%02.0f",year,month,day,hour,min,sec);
return;
}
int properties(kep_t K,xyz_t obspos,xyz_t sunpos,float radius,float t,double *ra,double *de,double *r,float *mag)
{
double st,ct,sn,cn,si,ci;
xyz_t satpos;
double dx,dy,dz,rsun,rearth,psun,pearth,p;
float phase;
int state;
// Angles
sn=sin(K.ascn);
cn=cos(K.ascn);
si=sin(K.eqinc);
ci=cos(K.eqinc);
st=sin(K.theta+t);
ct=cos(K.theta+t);
satpos.x=-sn*ci*st+cn*ct;
satpos.y=cn*ci*st+sn*ct;
satpos.z=si*st;
satpos.x*=(radius+XKMPER);
satpos.y*=(radius+XKMPER);
satpos.z*=(radius+XKMPER);
// 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);
// Angles
psun=asin(696.0e3/rsun)*R2D;
pearth=asin(6378.135/rearth)*R2D;
pearth=asin(6378.135/rearth)*R2D;
p=acos((-dx*satpos.x-dy*satpos.y-dz*satpos.z)/(rsun*rearth))*R2D;
p-=pearth;
// Position differences
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=atan2(dy,dx)*R2D;
*de=asin(dz/ *r)*R2D;
// Visibility
if (p<-psun) {
// strcpy(state,"eclipsed");
state=0;
} else if (p>-psun && p<psun) {
// strcpy(state,"umbra");
state=1;
} else if (p>psun) {
// strcpy(state,"sunlit");
state=2;
}
// Phase
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* *r))*R2D;
// Magnitude
if (state==2)
*mag=STDMAG-15.0+5*log10( *r)-2.5*log10(sin(phase*D2R)+(M_PI-phase*D2R)*cos(phase*D2R));
else
*mag=15;
return state;
}
// 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;
}

View File

@ -241,6 +241,12 @@ struct point decode_iod_observation(char *iod_line)
sscanf(iod_line+16,"%4d",&site_id);
s=get_site(site_id);
// Skip if site not found
if (s.id<0) {
fprintf(stderr,"Site %d not found!\n",site_id);
p.flag=0;
}
// Decode date/time
sscanf(iod_line+23,"%4d%2d%2d%2d%2d%5s",&year,&month,&iday,&hour,&min,secbuf);
sec=atof(secbuf);
@ -298,7 +304,7 @@ struct point decode_iod_observation(char *iod_line)
de=sign*(de+dd/100+ds/10000);
break;
default :
printf("IOD Format not implemented\n");
fprintf(stderr,"IOD Format not implemented\n");
p.flag=0;
break;
}
@ -315,7 +321,7 @@ struct point decode_iod_observation(char *iod_line)
} else if (epoch==5) {
mjd0=51544.5;
} else {
printf("Observing epoch not implemented\n");
fprintf(stderr,"Observing epoch not implemented\n");
p.flag=0;
}
@ -373,6 +379,9 @@ struct site get_site(int site_id)
}
fclose(file);
if (id!=site_id)
s.id==-1;
return s;
}
// Return x modulo y [0,y)

View File

@ -7,7 +7,7 @@
#include "satutl.h"
#include <getopt.h>
#define LIM 128
#define LIM 1024
#define XKE 0.07436680 // Guassian Gravitational Constant
#define XKMPER 6378.135
#define AE 1.0
@ -342,7 +342,7 @@ int fgetline(FILE *file,char *s,int lim)
int main(int argc,char *argv[])
{
int imode,satno=99000,arg;
int imode,satno=99000,arg,gmat=0;
FILE *file;
orbit_t orb;
xyz_t r,v;
@ -352,7 +352,7 @@ int main(int argc,char *argv[])
char *env;
// Decode options
while ((arg=getopt(argc,argv,"p:i:d:"))!=-1) {
while ((arg=getopt(argc,argv,"p:i:d:g"))!=-1) {
switch (arg) {
case 'p':
@ -372,6 +372,10 @@ int main(int argc,char *argv[])
return 0;
break;
case 'g':
gmat=1;
break;
default:
usage();
return 0;
@ -384,6 +388,11 @@ int main(int argc,char *argv[])
while (fgetline(file,line,LIM)>0) {
sscanf(line,"%lf %lf %lf %lf %lf %lf %lf",&mjd,&r.x,&r.y,&r.z,&v.x,&v.y,&v.z);
// Convert to MJD
if (gmat==1)
mjd+=29999.5;
// Convert
orb=rv2el(orb.satno,mjd,r,v);
orb.satno=satno;

View File

@ -654,9 +654,12 @@ int main(int argc,char *argv[])
// Adjust
if (adjust==1) {
// Count observations
for (i=0,nobs=0;i<d.n;i++)
if (d.p[i].flag==2)
nobs++;
for (i=0,nobs=0;i<d.n;i++) {
if (d.p[i].flag==1) {
d.p[i].flag=2;
nobs++;
}
}
// Not enough observations
if (nobs<10)
@ -679,7 +682,7 @@ int main(int argc,char *argv[])
adjust_fit(16);
// Print results
printf("%05d %8.3f %8.3f %8.3f %s %8.3f %d\n",satno,DEG(orb.mnan-orb0.mnan),DEG(orb.ascn-orb0.ascn),d.rms,ioddatafile,mjdmin-(SGDP4_jd0-2400000.5),d.nsel);
// printf("%05d %8.3f %8.3f %8.3f %s %8.3f %d\n",satno,DEG(orb.mnan-orb0.mnan),DEG(orb.ascn-orb0.ascn),d.rms,ioddatafile,mjdmin-(SGDP4_jd0-2400000.5),d.nsel);
// Dump tle
if (tleout==1)
@ -1065,6 +1068,9 @@ struct site get_site(int site_id)
}
fclose(file);
if (id!=site_id)
s.id==-1;
return s;
}
@ -1220,6 +1226,12 @@ struct point decode_iod_observation(char *iod_line)
sscanf(iod_line+16,"%4d",&site_id);
s=get_site(site_id);
// Skip if site not found
if (s.id<0) {
fprintf(stderr,"Site %d not found!\n",site_id);
p.flag=0;
}
// Decode date/time
sscanf(iod_line+23,"%4d%2d%2d%2d%2d%5s",&year,&month,&iday,&hour,&min,secbuf);
sec=atof(secbuf);

View File

@ -127,7 +127,7 @@ void plot_satellites(char *tlefile,struct image img,long satno,double mjd0,float
// Open TLE file
fp=fopen(tlefile,"rb");
if (fp==NULL)
fatal_error("File open failed for reading %s\n",tlefile);
return;
cpgsci(color);
@ -339,7 +339,7 @@ int main(int argc,char *argv[])
plot_satellites(filename,img,0,img.mjd,img.exptime,3);
sprintf(filename,"%s/catalog.tle",env);
plot_satellites(filename,img,0,img.mjd,img.exptime,0);
sprintf(filename,"/home/bassa/jsc.txt");
sprintf(filename,"%s/jsc.txt",env);
plot_satellites(filename,img,0,img.mjd,img.exptime,5);
cpgend();

View File

@ -967,7 +967,7 @@ void plot_notam(char *filename)
return;
}
void plot_map(void)
void plot_map(int plotflag)
{
int redraw=1,status;
char text[256];
@ -1060,6 +1060,11 @@ void plot_map(void)
// Reset redraw
redraw=0;
if (plotflag==1) {
cpgend();
exit(0);
}
// Get cursor
cpgcurs(&x,&y,&c);
@ -1186,14 +1191,20 @@ void plot_map(void)
int main(int argc,char *argv[])
{
int arg=0;
int arg=0,plotflag=0;
char plottype[128]="/xs";
// Initialize setup
initialize_setup();
// Decode options
while ((arg=getopt(argc,argv,"t:c:i:s:l:hN:p:mL:B:R:Sq"))!=-1) {
while ((arg=getopt(argc,argv,"t:c:i:s:l:hN:p:mL:B:R:Sqg:"))!=-1) {
switch (arg) {
case 'g':
strcpy(plottype,optarg);
plotflag=1;
break;
case 't':
strcpy(m.nfd,optarg);
@ -1263,9 +1274,9 @@ int main(int argc,char *argv[])
read_globe();
cpgopen("/xs");
cpgopen(plottype);
plot_map();
plot_map(plotflag);
cpgend();

View File

@ -15,8 +15,8 @@ PARAMETERS_NAME $ST_DATADIR/sextractor/default.param # name of the file containi
DETECT_TYPE CCD # "CCD" or "PHOTO" (*)
FLAG_IMAGE flag.fits # filename for an input FLAG-image
DETECT_MINAREA 5 # minimum number of pixels above threshold
DETECT_THRESH 1.5 # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
ANALYSIS_THRESH 1.5 # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
DETECT_THRESH 2.5 # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
ANALYSIS_THRESH 2.5 # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
FILTER Y # apply filter for detection ("Y" or "N")?
FILTER_NAME $ST_DATADIR/sextractor/default.conv # name of the file containing the filter

View File

@ -1078,6 +1078,7 @@ double modulo(double x,double y)
return x;
}
// Read a line of maximum length int lim from file FILE into string s
int fgetline(FILE *file,char *s,int lim)
{

View File

@ -12,12 +12,16 @@
#define AE 1.0 /* Earth radius in "chosen units". */
#define XKE 0.743669161e-1
#define CK2 5.413080e-4 /* (0.5 * XJ2 * AE * AE) */
#define D2R M_PI/180.0
#define R2D 180.0/M_PI
extern double SGDP4_jd0;
void usage(void)
{
return;
}
double modulo(double x,double y);
// Compute Julian Day from Date
double date2mjd(int year,int month,double day)
@ -151,6 +155,23 @@ void mjd2nfd(double mjd,char *nfd)
return;
}
float orbital_longitude_at_midnight(orbit_t orb,double mjd0)
{
int rv,imode;
double jd,tsince,mjd;
kep_t K;
imode=init_sgdp4(&orb);
mjd=floor(mjd0);
jd=mjd+2400000.5;
tsince=1440.0*(jd-SGDP4_jd0);
rv=sgdp4(tsince,1,&K);
return modulo(K.theta*R2D,360.0);
}
int main(int argc,char *argv[])
{
int arg=0,satno=0,header=0,oneline=0,no,name=0,desig=0;
@ -158,7 +179,7 @@ int main(int argc,char *argv[])
char line0[LIM],line1[LIM],line2[LIM],nfd[32];
FILE *file;
orbit_t orb;
float aodp,perigee,apogee,period;
float aodp,perigee,apogee,period,lng;
int info=0;
double mjd;
char *env;
@ -167,7 +188,7 @@ int main(int argc,char *argv[])
sprintf(tlefile,"%s/bulk.tle",env);
// Decode options
while ((arg=getopt(argc,argv,"c:i:aH1ftnd"))!=-1) {
while ((arg=getopt(argc,argv,"c:i:aH1ftndb"))!=-1) {
switch (arg) {
case 'c':
@ -198,6 +219,10 @@ int main(int argc,char *argv[])
info=1;
break;
case 'b':
info=2;
break;
case 'H':
header=1;
break;
@ -281,6 +306,10 @@ int main(int argc,char *argv[])
mjd2nfd(mjd,nfd);
if (info==0) printf("%05d %10.4lf %8.4f %8.4f %8.4f %8.4f %8.6f %8.5f\n",orb.satno,mjd,DEG(orb.eqinc),DEG(orb.ascn),DEG(orb.argp),DEG(orb.mnan),orb.ecc,orb.rev);
if (info==1) printf("%05d %6.0f x %6.0f x %6.2f %8.2f %8.6f %14.8lf\n",orb.satno,perigee,apogee,DEG(orb.eqinc),period,orb.ecc,mjd);
if (info==2) {
lng=orbital_longitude_at_midnight(orb,mjd);
printf("%05d %10.4lf %8.4f %8.4f %8.4f %8.4f %8.6f %8.5f %10.4lf %8.4f\n",orb.satno,mjd,DEG(orb.eqinc),DEG(orb.ascn),DEG(orb.argp),DEG(orb.mnan),orb.ecc,orb.rev,floor(mjd),lng);
}
}
fclose(file);
} else if (oneline==2) {
@ -302,3 +331,12 @@ int main(int argc,char *argv[])
return 0;
}
// Return x modulo y [0,y)
double modulo(double x,double y)
{
x=fmod(x,y);
if (x<0.0) x+=y;
return x;
}