Compute geostationary longitudes
parent
ff5456a053
commit
cba5394ebf
|
@ -0,0 +1,204 @@
|
|||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <time.h>
|
||||
#include <getopt.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
|
||||
|
||||
long Isat=0;
|
||||
long Isatsel=0;
|
||||
extern double SGDP4_jd0;
|
||||
struct map {
|
||||
long satno;
|
||||
double mjd;
|
||||
char nfd[LIM],tlefile[LIM],observer[32];
|
||||
char datadir[LIM],tledir[LIM];
|
||||
} m;
|
||||
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
void usage()
|
||||
{
|
||||
printf("Usage goes here.\n");
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
// Return x modulo y [0,y)
|
||||
double modulo(double x,double y)
|
||||
{
|
||||
x=fmod(x,y);
|
||||
if (x<0.0) x+=y;
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
// Compute longitude
|
||||
void compute_longitude(char *tlefile,long satno,double mjd)
|
||||
{
|
||||
FILE *fp=NULL;
|
||||
orbit_t orb;
|
||||
xyz_t satpos,satvel;
|
||||
double jd,h,l,b,r;
|
||||
long imode;
|
||||
|
||||
// Open TLE file
|
||||
fp = fopen(tlefile, "rb");
|
||||
if(fp == NULL)
|
||||
fatal_error("File open failed for reading \"%s\"", tlefile);
|
||||
|
||||
// Loop over elements
|
||||
while(read_twoline(fp, satno, &orb) == 0) {
|
||||
Isat = orb.satno;
|
||||
imode = init_sgdp4(&orb);
|
||||
if(imode == SGDP4_ERROR) continue;
|
||||
|
||||
// Skip objects with mean motions outside of 0.8 to 1.2 revs/day
|
||||
if (orb.rev<0.8 || orb.rev>1.2)
|
||||
continue;
|
||||
|
||||
// Get Julian Date
|
||||
jd=mjd+2400000.5;
|
||||
|
||||
// Get positions
|
||||
satpos_xyz(jd,&satpos,&satvel);
|
||||
|
||||
// Greenwich Mean Sidereal time
|
||||
h=gmst(mjd);
|
||||
|
||||
// Celestial position
|
||||
r=sqrt(satpos.x*satpos.x+satpos.y*satpos.y+satpos.z*satpos.z);
|
||||
l=atan2(satpos.y,satpos.x)*R2D;
|
||||
l=modulo(l-h,360.0);
|
||||
b=asin(satpos.z/r)*R2D;
|
||||
if (l>180.0)
|
||||
l-=360.0;
|
||||
if (l<-180.0)
|
||||
l+=360.0;
|
||||
|
||||
printf("%05d %10s %8.3lf %8.3lf %6.0lf\n",orb.satno,orb.desig,l,b,r-XKMPER);
|
||||
}
|
||||
fclose(fp);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
int main(int argc,char *argv[])
|
||||
{
|
||||
int arg=0;
|
||||
long satno=0;
|
||||
double mjd;
|
||||
char nfd[LIM],tlefile[LIM];
|
||||
|
||||
nfd_now(nfd);
|
||||
mjd=nfd2mjd(nfd);
|
||||
|
||||
|
||||
// Decode options
|
||||
if (argc>1) {
|
||||
while ((arg=getopt(argc,argv,"t:c:i:h"))!=-1) {
|
||||
switch (arg) {
|
||||
|
||||
case 't':
|
||||
strcpy(nfd,optarg);
|
||||
mjd=nfd2mjd(nfd);
|
||||
break;
|
||||
|
||||
case 'c':
|
||||
strcpy(tlefile,optarg);
|
||||
break;
|
||||
|
||||
case 'i':
|
||||
satno=atoi(optarg);
|
||||
break;
|
||||
|
||||
case 'h':
|
||||
usage();
|
||||
return 0;
|
||||
break;
|
||||
|
||||
default:
|
||||
usage();
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
usage();
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Compute longitudes of satellites
|
||||
compute_longitude(tlefile,satno,mjd);
|
||||
|
||||
return 0;
|
||||
}
|
Loading…
Reference in New Issue