2013-05-18 11:54:11 -06:00
# include <stdio.h>
# include <string.h>
# include <stdlib.h>
# include <math.h>
2018-02-27 14:36:04 -07:00
# include <ctype.h>
# include <cpgplot.h>
# include <wcslib/cel.h>
2013-05-18 11:54:11 -06:00
# include "sgdp4h.h"
# include <getopt.h>
# define LIM 80
# define NMAX 256
# 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)
2014-09-24 01:27:45 -06:00
# define XKE 0.07436680 // Guassian Gravitational Constant
# define AE 1.0
# define XMNPDA 1440.0
2013-05-18 11:54:11 -06:00
long Isat = 0 ;
long Isatsel = 0 ;
extern double SGDP4_jd0 ;
struct point {
int flag , satno ;
2013-10-24 00:00:12 -06:00
char desig [ 10 ] ;
2013-05-18 11:54:11 -06:00
double mjd , ra , de , rac , dec ;
float st , sr ;
char iod_line [ LIM ] ;
double dx , dy , dr , dt ;
xyz_t obspos ;
} ;
2015-06-08 22:55:38 -06:00
struct doppler {
double mjd , freq , v ;
xyz_t obspos , obsvel ;
} ;
2013-05-18 11:54:11 -06:00
struct data {
2015-06-08 22:55:38 -06:00
int n , nsel , m ;
2013-05-18 11:54:11 -06:00
struct point * p ;
2015-06-08 22:55:38 -06:00
struct doppler * q ;
2013-05-18 11:54:11 -06:00
double chisq , rms ;
} d ;
struct site {
int id ;
double lng , lat ;
float alt ;
char observer [ 64 ] ;
} ;
orbit_t orb ;
struct site get_site ( int site_id ) ;
struct point decode_iod_observation ( char * iod_line ) ;
2015-06-08 22:55:38 -06:00
struct doppler decode_doppler_observation ( char * line ) ;
2013-05-18 11:54:11 -06:00
int fgetline ( FILE * file , char * s , int lim ) ;
double modulo ( double x , double y ) ;
double gmst ( double mjd ) ;
double dgmst ( double mjd ) ;
double date2mjd ( int year , int month , double day ) ;
double mjd2doy ( double mjd , int * yr ) ;
void mjd2date ( double mjd , int * year , int * month , double * day ) ;
void obspos_xyz ( double mjd , double lng , double lat , float alt , xyz_t * pos , xyz_t * vel ) ;
void precess ( double mjd0 , double ra0 , double de0 , double mjd , double * ra , double * de ) ;
void forward ( double ra0 , double de0 , double ra , double de , double * x , double * y ) ;
2015-06-08 22:55:38 -06:00
struct data read_data ( char * iodfname , char * dopfname ) ;
2013-05-18 11:54:11 -06:00
void versafit ( int m , int n , double * a , double * da , double ( * func ) ( double * ) , double dchisq , double tol , char * opt ) ;
double chisq ( double a [ ] ) ;
orbit_t read_tle ( char * filename , int satno ) ;
void format_tle ( orbit_t orb , char * line1 , char * line2 ) ;
void highlight ( float x0 , float y0 , float x , float y , int flag ) ;
void time_range ( double * mjdmin , double * mjdmax , int flag ) ;
void print_tle ( orbit_t orb , char * filename ) ;
void fit ( orbit_t orb , int * ia ) ;
void usage ( ) ;
2014-09-24 01:27:45 -06:00
double nfd2mjd ( char * date ) ;
double dot ( xyz_t a , xyz_t b ) ;
double magnitude ( xyz_t a ) ;
xyz_t cross ( xyz_t a , xyz_t b ) ;
orbit_t classel ( int ep_year , double ep_day , xyz_t r , xyz_t v ) ;
orbit_t rv2el ( int satno , double mjd , xyz_t r0 , xyz_t v0 ) ;
// Propagate elements
void propagate ( double mjd )
{
int imode ;
xyz_t r , v ;
char desig [ 20 ] , line1 [ 80 ] , line2 [ 80 ] ;
// Store designation
strcpy ( desig , orb . desig ) ;
// Propagate
imode = init_sgdp4 ( & orb ) ;
imode = satpos_xyz ( mjd + 2400000.5 , & r , & v ) ;
// Convert
orb = rv2el ( orb . satno , mjd , r , v ) ;
// Copy designation
strcpy ( orb . desig , desig ) ;
return ;
}
2013-05-18 11:54:11 -06:00
xyz_t get_position ( double r0 , int i0 )
{
int i ;
double rr , drr , r ;
xyz_t pos ;
double x , y , z ;
// Initial range
rr = 100.0 ;
do {
x = d . p [ i0 ] . obspos . x + rr * cos ( d . p [ i0 ] . ra * D2R ) * cos ( d . p [ i0 ] . de * D2R ) ;
y = d . p [ i0 ] . obspos . y + rr * sin ( d . p [ i0 ] . ra * D2R ) * cos ( d . p [ i0 ] . de * D2R ) ;
z = d . p [ i0 ] . obspos . z + rr * sin ( d . p [ i0 ] . de * D2R ) ;
r = sqrt ( x * x + y * y + z * z ) ;
drr = r - r0 ;
rr - = drr ;
} while ( fabs ( drr ) > 0.01 ) ;
pos . x = x ;
pos . y = y ;
pos . z = z ;
return pos ;
}
2018-02-27 14:36:04 -07:00
void period_search ( void )
2013-08-31 01:49:06 -06:00
{
int i , j , i1 , i2 ;
float dt ;
int nrev , nrevmin , nrevmax ;
char line1 [ 70 ] , line2 [ 70 ] ;
int ia [ 7 ] ;
// Set fitting parameters
for ( i = 0 ; i < 6 ; i + + )
ia [ i ] = 1 ;
ia [ 6 ] = 0 ;
// Select all points
for ( i = 0 ; i < d . n ; i + + )
d . p [ i ] . flag = 2 ;
// Print observations
printf ( " Present observations: \n " ) ;
for ( i = 0 ; i < d . n ; i + + )
printf ( " %3d: %s \n " , i + 1 , d . p [ i ] . iod_line ) ;
printf ( " \n Select center observations of both arcs: " ) ;
scanf ( " %d %d " , & i1 , & i2 ) ;
dt = d . p [ i2 ] . mjd - d . p [ i1 ] . mjd ;
printf ( " \n Time passed: %f days \n " , dt ) ;
printf ( " Provide revolution range to search over [min,max]: " ) ;
scanf ( " %d %d " , & nrevmin , & nrevmax ) ;
for ( nrev = nrevmin ; nrev < nrevmax + 1 ; nrev + + ) {
orb . satno = 79000 + nrev ;
orb . rev = nrev / dt ;
// Set parameters
for ( i = 0 ; i < 7 ; i + + )
ia [ i ] = 0 ;
// Loop over parameters
for ( i = 0 ; i < 6 ; i + + ) {
if ( i = = 0 ) ia [ 4 ] = 1 ;
if ( i = = 1 ) ia [ 1 ] = 1 ;
if ( i = = 2 ) ia [ 0 ] = 1 ;
if ( i = = 3 ) ia [ 3 ] = 1 ;
if ( i = = 4 ) ia [ 2 ] = 1 ;
if ( i = = 5 ) ia [ 5 ] = 1 ;
for ( j = 0 ; j < 5 ; j + + )
fit ( orb , ia ) ;
}
format_tle ( orb , line1 , line2 ) ;
printf ( " %s \n %s \n # %d revs, %f revs/day, %f \n " , line1 , line2 , nrev , nrev / dt , d . rms ) ;
}
return ;
}
2014-05-25 02:53:15 -06:00
int ewsearch ( void )
{
int i , satno = 88300 ;
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 , * tlefile ;
int year , month ;
double day ;
// Provide
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 2: get time range
time_range ( & mjdmin , & mjdmax , 2 ) ;
// Open files
file = fopen ( " search.dat " , " w " ) ;
tlefile = fopen ( " search.tle " , " w " ) ;
// Step 4: Loop over eccentricity
for ( ecc = eccmin ; ecc < eccmax ; ecc + = decc ) {
for ( argp = argpmin ; argp < argpmax ; argp + = dargp ) {
orb . satno = satno + + ;
orb . ecc = ecc ;
orb . argp = argp * D2R ;
// Set parameters
for ( i = 0 ; i < 7 ; i + + )
ia [ i ] = 0 ;
// Step 4: loop over parameters
for ( i = 0 ; i < 4 ; 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 ) ;
format_tle ( orb , line1 , line2 ) ;
fprintf ( tlefile , " OBJ \n %s \n %s \n " , line1 , line2 ) ;
mjd2date ( mjdmin , & year , & month , & day ) ;
fprintf ( tlefile , " # %4d%02d%05.2lf- " , year , month , day ) ;
mjd2date ( mjdmax , & year , & month , & day ) ;
fprintf ( tlefile , " %4d%02d%05.2lf, %d measurements, %.5lf deg rms \n " , year , month , day , d . n , d . rms ) ;
}
fprintf ( file , " \n " ) ;
}
fclose ( file ) ;
fclose ( tlefile ) ;
return orb . satno ;
}
2014-03-26 07:51:34 -06:00
2014-02-23 04:14:38 -07:00
int revsearch ( void )
{
2014-03-06 10:10:57 -07:00
int i , satno = 88300 ;
2014-02-23 04:14:38 -07:00
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 , * tlefile ;
int year , month ;
double day ;
// Provide
printf ( " Mean motion [min, max, stepsize]: \n " ) ;
scanf ( " %lf %lf %lf " , & revmin , & revmax , & drev ) ;
// 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 ) ;
// Open files
file = fopen ( " search.dat " , " w " ) ;
tlefile = fopen ( " search.tle " , " w " ) ;
// Step 4: Loop over eccentricity
for ( rev = revmin ; rev < revmax ; rev + = drev ) {
2014-03-04 09:10:05 -07:00
orb . satno = satno + + ;
2014-02-23 04:14:38 -07:00
orb . rev = rev ;
// 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 [ 2 ] = 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 ) ;
format_tle ( orb , line1 , line2 ) ;
fprintf ( tlefile , " OBJ \n %s \n %s \n " , line1 , line2 ) ;
mjd2date ( mjdmin , & year , & month , & day ) ;
fprintf ( tlefile , " # %4d%02d%05.2lf- " , year , month , day ) ;
mjd2date ( mjdmax , & year , & month , & day ) ;
fprintf ( tlefile , " %4d%02d%05.2lf, %d measurements, %.5lf deg rms \n " , year , month , day , d . n , d . rms ) ;
}
fclose ( file ) ;
fclose ( tlefile ) ;
return orb . satno ;
}
2013-07-10 03:54:14 -06:00
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 ;
}
2013-05-18 11:54:11 -06:00
int circular_fit ( void )
{
int i ;
double mjdmin , mjdmax ;
int ia [ 7 ] = { 0 , 0 , 0 , 0 , 0 , 0 , 0 } ;
// 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 ) ;
// Step 3: set initial orbit
orb . satno = d . p [ 0 ] . satno ;
orb . eqinc = 0.5 * M_PI ;
orb . ascn = 0.0 ;
2014-09-17 08:42:16 -06:00
orb . ecc = 0.0001 ;
2013-05-18 11:54:11 -06:00
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 ) ;
// Step 4: loop over parameters
for ( i = 0 ; i < 4 ; 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 ;
// Do fit
fit ( orb , ia ) ;
}
fit ( orb , ia ) ;
return orb . satno ;
}
2014-09-24 01:27:45 -06:00
int adjust_fit ( int depth )
2013-05-18 11:54:11 -06:00
{
int i ;
double mjdmin , mjdmax ;
int ia [ 6 ] = { 0 , 0 , 0 , 0 , 0 , 0 } ;
// Step 2: loop over parameters
2014-09-24 01:27:45 -06:00
for ( i = 0 ; i < depth ; i + + ) {
2013-05-18 11:54:11 -06:00
if ( i = = 0 ) ia [ 4 ] = 1 ;
if ( i = = 1 ) ia [ 1 ] = 1 ;
2014-09-24 01:27:45 -06:00
if ( i = = 2 ) ia [ 0 ] = 1 ;
if ( i = = 3 ) ia [ 5 ] = 1 ;
if ( i = = 4 ) ia [ 2 ] = 1 ;
if ( i = = 5 ) ia [ 3 ] = 1 ;
2013-05-18 11:54:11 -06:00
// Do fit
fit ( orb , ia ) ;
}
fit ( orb , ia ) ;
return orb . satno ;
}
void old_circular_fit ( void )
{
int i , j , i0 , i1 ;
float r0 = 6500 ;
xyz_t pos0 , pos1 ;
double ang , dt , w , w0 ;
// Get end points
for ( i = 0 , j = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . flag = = 2 ) {
if ( j = = 0 )
i0 = i ;
i1 = i ;
j + + ;
}
}
// Time difference
dt = 86400.0 * ( d . p [ i1 ] . mjd - d . p [ i0 ] . mjd ) ;
i = 0 ;
do {
w0 = 360.0 / ( 2.0 * M_PI * sqrt ( r0 * r0 * r0 / 398600 ) ) ;
// Get positions
pos0 = get_position ( r0 , i0 ) ;
pos1 = get_position ( r0 , i1 ) ;
// Compute angle
ang = acos ( ( pos0 . x * pos1 . x + pos0 . y * pos1 . y + pos0 . z * pos1 . z ) / ( r0 * r0 ) ) * R2D ;
// Angular motion (deg/sec);
w = ang / dt ;
r0 + = 1000.0 * ( w0 - w ) ;
i + + ;
} while ( fabs ( w0 - w ) > 1e-5 & & i < 1000 ) ;
printf ( " %f \n " , r0 ) ;
return ;
}
int main ( int argc , char * argv [ ] )
{
int i , j , nobs = 0 ;
2014-02-24 00:40:07 -07:00
int redraw = 1 , plot_residuals = 0 , adjust = 0 , quit = 0 , identify = 0 ;
2013-05-18 11:54:11 -06:00
int ia [ ] = { 0 , 0 , 0 , 0 , 0 , 0 , 0 } ;
float dx [ ] = { 0.1 , 0.1 , 0.35 , 0.35 , 0.6 , 0.6 , 0.85 } , dy [ ] = { 0.0 , - 0.25 , 0.0 , - 0.25 , 0.0 , - 0.25 , 0.0 } ;
char c ;
int mode = 0 , posn = 0 , click = 0 ;
float x0 , y0 , x , y ;
float xmin = 0.0 , xmax = 360.0 , ymin = - 90.0 , ymax = 90.0 ;
char string [ 64 ] , bstar [ 10 ] = " 50000-4 " , line0 [ 72 ] , line1 [ 72 ] , line2 [ 72 ] , text [ 10 ] ;
2014-09-24 01:27:45 -06:00
char filename [ 64 ] , nfd [ 32 ] ;
2013-05-18 11:54:11 -06:00
int satno = - 1 ;
2014-09-24 01:27:45 -06:00
double mjdmin , mjdmax , mjd = 0.0 ;
int length = 864000 ;
2013-07-10 03:54:14 -06:00
int arg = 0 , elset = 0 , circular = 0 , tleout = 0 , noplot = 0 ;
2015-06-08 22:55:38 -06:00
char * ioddatafile = NULL , * dopdatafile = NULL , * catalog , tlefile [ LIM ] ;
2013-05-18 11:54:11 -06:00
orbit_t orb0 ;
2014-02-24 00:40:07 -07:00
FILE * file ;
2013-05-18 11:54:11 -06:00
// Decode options
2015-06-08 22:55:38 -06:00
while ( ( arg = getopt ( argc , argv , " d:r:c:i:haCo:pIt:l:m: " ) ) ! = - 1 ) {
2013-05-18 11:54:11 -06:00
switch ( arg ) {
2015-06-08 22:55:38 -06:00
2013-05-18 11:54:11 -06:00
case ' d ' :
2015-06-08 22:55:38 -06:00
ioddatafile = optarg ;
break ;
case ' r ' :
dopdatafile = optarg ;
2013-05-18 11:54:11 -06:00
break ;
case ' c ' :
catalog = optarg ;
break ;
case ' i ' :
satno = atoi ( optarg ) ;
break ;
2014-09-24 01:27:45 -06:00
case ' t ' :
strcpy ( nfd , optarg ) ;
mjd = nfd2mjd ( nfd ) ;
break ;
case ' m ' :
mjd = ( double ) 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 ;
2013-07-10 03:54:14 -06:00
case ' C ' :
circular = 1 ;
break ;
case ' p ' :
noplot = 1 ;
break ;
case ' o ' :
tleout = 1 ;
strcpy ( tlefile , optarg ) ;
break ;
2013-05-18 11:54:11 -06:00
case ' h ' :
usage ( ) ;
return 0 ;
break ;
case ' a ' :
adjust = 1 ;
break ;
2014-02-24 00:40:07 -07:00
case ' I ' :
identify = 1 ;
break ;
2013-05-18 11:54:11 -06:00
default :
usage ( ) ;
return 0 ;
}
}
// Read data
2015-06-08 22:55:38 -06:00
d = read_data ( ioddatafile , dopdatafile ) ;
2014-09-24 01:27:45 -06:00
// Select observations based on time
if ( mjd > 0.0 ) {
mjdmin = mjd - length / 86400.0 ;
mjdmax = mjd ;
for ( i = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . mjd > mjdmin & & d . p [ i ] . mjd < = mjdmax & & d . p [ i ] . flag = = 1 )
d . p [ i ] . flag = 2 ;
else
d . p [ i ] . flag = 0 ;
}
}
2013-05-18 11:54:11 -06:00
time_range ( & mjdmin , & mjdmax , 1 ) ;
// Read TLE
2014-02-24 00:40:07 -07:00
if ( satno > = 0 )
2013-05-18 11:54:11 -06:00
orb = read_tle ( catalog , satno ) ;
2014-02-24 00:40:07 -07:00
// Open catalog
if ( identify = = 1 ) {
file = fopen ( catalog , " r " ) ;
if ( file = = NULL )
fatal_error ( " Failed to open %s \n " , catalog ) ;
2013-05-18 11:54:11 -06:00
}
2014-02-24 00:40:07 -07:00
// Reloop stderr
2013-05-18 11:54:11 -06:00
freopen ( " /tmp/stderr.txt " , " w " , stderr ) ;
2014-02-24 00:40:07 -07:00
2013-07-10 03:54:14 -06:00
// 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 ) ;
}
2014-02-24 00:40:07 -07:00
// Identify
if ( identify = = 1 ) {
2014-09-24 01:27:45 -06:00
// Select all points
for ( i = 0 ; i < d . n ; i + + )
d . p [ i ] . flag = 2 ;
2014-09-25 00:48:34 -06:00
// Reset satno
if ( satno = = - 1 )
satno = 0 ;
2014-02-24 00:40:07 -07:00
// Loop over file
2014-09-25 00:48:34 -06:00
while ( read_twoline ( file , satno , & orb ) = = 0 ) {
2014-02-24 00:40:07 -07:00
orb0 = orb ;
2014-09-24 01:27:45 -06:00
adjust_fit ( 2 ) ;
2014-02-24 00:40:07 -07:00
fit ( orb , ia ) ;
2015-06-08 22:55:38 -06:00
printf ( " %05d %8.3f %8.3f %8.3f %s %8.3f \n " , orb . satno , DEG ( orb . mnan - orb0 . mnan ) , DEG ( orb . ascn - orb0 . ascn ) , d . rms , ioddatafile , mjdmin - ( SGDP4_jd0 - 2400000.5 ) ) ;
2014-02-24 00:40:07 -07:00
}
fclose ( file ) ;
plot_residuals = 1 ;
redraw = 1 ;
quit = 1 ;
noplot = 1 ;
}
2013-05-18 11:54:11 -06:00
// Adjust
if ( adjust = = 1 ) {
2014-09-24 01:27:45 -06:00
// Count observations
2016-07-30 07:36:44 -06:00
for ( i = 0 , nobs = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . flag = = 1 ) {
d . p [ i ] . flag = 2 ;
nobs + + ;
}
}
2014-09-24 01:27:45 -06:00
// Not enough observations
if ( nobs < 10 )
return 0 ;
// Get last point
for ( i = 0 , j = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . flag = = 2 ) {
if ( j = = 0 ) mjdmax = d . p [ i ] . mjd ;
if ( d . p [ i ] . mjd > mjdmax ) mjdmax = d . p [ i ] . mjd ;
j + + ;
}
}
// Propagate orbit to time of last observation
propagate ( mjdmax ) ;
2013-05-18 11:54:11 -06:00
orb0 = orb ;
2014-09-24 01:27:45 -06:00
// Fit
adjust_fit ( 16 ) ;
// Print results
2016-07-30 07:36:44 -06:00
// 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);
2013-07-10 03:54:14 -06:00
// Dump tle
if ( tleout = = 1 )
print_tle ( orb , tlefile ) ;
2014-09-24 01:27:45 -06:00
// Set thingies
plot_residuals = 1 ;
redraw = 1 ;
quit = 1 ;
2013-07-10 03:54:14 -06:00
}
// Exit before plotting
if ( quit = = 1 & & noplot = = 1 ) {
2015-06-08 22:55:38 -06:00
if ( d . n ) free ( d . p ) ;
if ( d . m ) free ( d . q ) ;
2013-07-10 03:54:14 -06:00
fclose ( stderr ) ;
return 0 ;
2013-05-18 11:54:11 -06:00
}
cpgopen ( " /xs " ) ;
cpgask ( 0 ) ;
// For ever loop
for ( ; ; ) {
if ( redraw = = 1 ) {
cpgpage ( ) ;
cpgsvp ( 0.1 , 0.95 , 0.0 , 0.18 ) ;
cpgswin ( 0.0 , 1.0 , - 0.5 , 0.5 ) ;
// Buttons
cpgtext ( 0.12 , - 0.05 , " Inclination " ) ;
cpgtext ( 0.372 , - 0.05 , " Eccentricity " ) ;
cpgtext ( 0.62 , - 0.05 , " Mean Anomaly " ) ;
cpgtext ( 0.87 , - 0.05 , " B \\ u* \\ d " ) ;
cpgtext ( 0.12 , - 0.3 , " Ascending Node " ) ;
cpgtext ( 0.37 , - 0.3 , " Arg. of Perigee " ) ;
cpgtext ( 0.62 , - 0.3 , " Mean Motion " ) ;
// Toggles
for ( i = 0 ; i < 7 ; i + + ) {
cpgpt1 ( dx [ i ] , dy [ i ] , 19 ) ;
if ( ia [ i ] = = 1 ) {
cpgsci ( 2 ) ;
cpgpt1 ( dx [ i ] , dy [ i ] , 16 ) ;
cpgsci ( 1 ) ;
}
}
// Plot map
2014-09-25 00:48:34 -06:00
cpgsvp ( 0.1 , 0.9 , 0.2 , 0.9 ) ;
2013-05-18 11:54:11 -06:00
cpgswin ( xmax , xmin , ymin , ymax ) ;
cpgbox ( " BCTSN " , 0. , 0 , " BCTSN " , 0. , 0 ) ;
cpglab ( " Right Ascension " , " Declination " , " " ) ;
if ( satno > 0 ) {
// Plot tle
format_tle ( orb , line1 , line2 ) ;
cpgmtxt ( " T " , 2.0 , 0.0 , 0.0 , line1 ) ;
cpgmtxt ( " T " , 1.0 , 0.0 , 0.0 , line2 ) ;
}
// Plot points
for ( i = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . flag > = 1 ) {
cpgpt1 ( d . p [ i ] . ra , d . p [ i ] . de , 17 ) ;
sprintf ( text , " %d " , i + 1 ) ;
cpgtext ( d . p [ i ] . ra , d . p [ i ] . de , text ) ;
if ( plot_residuals = = 1 ) {
cpgmove ( d . p [ i ] . ra , d . p [ i ] . de ) ;
cpgdraw ( d . p [ i ] . rac , d . p [ i ] . dec ) ;
}
if ( d . p [ i ] . flag = = 2 ) {
cpgsci ( 2 ) ;
cpgpt1 ( d . p [ i ] . ra , d . p [ i ] . de , 4 ) ;
cpgsci ( 1 ) ;
}
}
}
}
// Quit
if ( quit = = 1 )
break ;
// Get cursor
cpgband ( mode , posn , x0 , y0 , & x , & y , & c ) ;
2013-10-24 00:00:12 -06:00
// Help
if ( c = = ' h ' | | c = = ' ? ' ) {
printf ( " q Quit \n w Write TLE file \n P Search mean motion space \n f Fit orbit \n s Select observations in current window \n u Unselect observations \n 1-7 Toggle parameter \n C Fit circular orbit \n S Search mean motion/eccentricity space \n c Change a parameter \n z Zoom \n r Unzoom \n t Toggle starting orbits (LEO, GTO, GEO, HEO) \n " ) ;
}
2013-05-18 11:54:11 -06:00
// Quit
if ( c = = ' q ' | | c = = ' Q ' )
break ;
2013-08-31 01:49:06 -06:00
// Period search
if ( c = = ' P ' ) {
period_search ( ) ;
}
2013-05-18 11:54:11 -06:00
// Fit
if ( c = = ' f ' ) {
// Count points
for ( i = 0 , nobs = 0 ; i < d . n ; i + + )
if ( d . p [ i ] . flag = = 2 )
nobs + + ;
if ( satno < 0 ) {
printf ( " No elements loaded! \n " ) ;
} else if ( nobs = = 0 ) {
printf ( " No points selected! \n " ) ;
} else {
fit ( orb , ia ) ;
2013-08-31 01:49:06 -06:00
printf ( " %d %.5f \n " , nobs , d . rms ) ;
2013-05-18 11:54:11 -06:00
plot_residuals = 1 ;
redraw = 1 ;
continue ;
}
}
// Write TLE
if ( c = = ' w ' ) {
printf ( " TLE filename to write: " ) ;
scanf ( " %s " , filename ) ;
print_tle ( orb , filename ) ;
printf ( " \n ================================================================================ \n " ) ;
continue ;
}
// Highlight
2013-10-24 00:00:12 -06:00
if ( c = = ' s ' ) {
2013-05-18 11:54:11 -06:00
highlight ( xmin , ymin , xmax , ymax , 2 ) ;
time_range ( & mjdmin , & mjdmax , 2 ) ;
for ( i = 0 , nobs = 0 ; i < d . n ; i + + )
if ( d . p [ i ] . flag = = 2 )
nobs + + ;
click = 0 ;
mode = 0 ;
redraw = 1 ;
continue ;
}
// Unselect
if ( c = = ' U ' ) {
for ( i = 0 ; i < d . n ; i + + )
d . p [ i ] . flag = 1 ;
time_range ( & mjdmin , & mjdmax , 1 ) ;
redraw = 1 ;
continue ;
}
// Unselect
if ( c = = ' u ' ) {
for ( i = 0 ; i < d . n ; i + + )
if ( d . p [ i ] . flag = = 2 )
d . p [ i ] . flag = 1 ;
redraw = 1 ;
continue ;
}
// Toggles
if ( isdigit ( c ) & & c - ' 0 ' > = 1 & & c - ' 0 ' < 8 ) {
if ( ia [ c - 49 ] = = 0 )
ia [ c - 49 ] = 1 ;
else if ( ia [ c - 49 ] = = 1 )
ia [ c - 49 ] = 0 ;
redraw = 1 ;
continue ;
}
// Circular fit
if ( c = = ' C ' ) {
satno = circular_fit ( ) ;
plot_residuals = 1 ;
printf ( " %.3f \n " , d . rms ) ;
ia [ 0 ] = ia [ 1 ] = ia [ 4 ] = ia [ 5 ] = 1 ;
redraw = 1 ;
}
2013-07-10 03:54:14 -06:00
// Search
if ( c = = ' S ' ) {
satno = psearch ( ) ;
plot_residuals = 1 ;
ia [ 0 ] = ia [ 1 ] = ia [ 4 ] = ia [ 5 ] = 1 ;
redraw = 1 ;
}
2014-02-23 04:14:38 -07:00
// Search
if ( c = = ' R ' ) {
satno = revsearch ( ) ;
plot_residuals = 1 ;
// ia[0]=ia[1]=ia[4]=ia[5]=ia[2]=1;
redraw = 1 ;
}
2014-05-25 02:53:15 -06:00
// EW search
if ( c = = ' e ' ) {
satno = ewsearch ( ) ;
plot_residuals = 1 ;
redraw = 1 ;
}
2013-05-18 11:54:11 -06:00
// Change
if ( c = = ' c ' ) {
2013-07-16 02:04:57 -06:00
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 (0) Sat ID \n Which parameter to change: " ) ;
2013-05-18 11:54:11 -06:00
scanf ( " %i " , & i ) ;
2013-07-16 02:04:57 -06:00
if ( i > = 0 & & i < = 9 ) {
2013-05-18 11:54:11 -06:00
printf ( " \n New value: " ) ;
fgets ( string , 64 , stdin ) ;
scanf ( " %s " , string ) ;
2013-10-24 00:00:12 -06:00
if ( i = = 0 ) strcpy ( orb . desig , string ) ;
2013-05-18 11:54:11 -06:00
if ( i = = 1 ) orb . eqinc = RAD ( atof ( string ) ) ;
if ( i = = 2 ) orb . ascn = RAD ( atof ( string ) ) ;
if ( i = = 3 ) orb . ecc = atof ( string ) ;
if ( i = = 4 ) orb . argp = RAD ( atof ( string ) ) ;
if ( i = = 5 ) orb . mnan = RAD ( atof ( string ) ) ;
if ( i = = 6 ) orb . rev = atof ( string ) ;
if ( i = = 7 ) orb . bstar = atof ( string ) ;
if ( i = = 8 ) {
orb . ep_year = 2000 + ( int ) floor ( atof ( string ) / 1000.0 ) ;
orb . ep_day = atof ( string ) - 1000 * floor ( atof ( string ) / 1000.0 ) ;
}
if ( i = = 9 ) orb . satno = atoi ( string ) ;
redraw = 1 ;
continue ;
}
printf ( " \n ================================================================================ \n " ) ;
}
// Zoom
if ( c = = ' z ' ) {
click = 1 ;
mode = 2 ;
}
// Execute zoom, or box delete
if ( c = = ' A ' ) {
if ( click = = 0 ) {
click = 1 ;
} else if ( click = = 1 & & mode = = 2 ) {
xmin = ( x0 < x ) ? x0 : x ;
xmax = ( x0 > x ) ? x0 : x ;
ymin = ( y0 < y ) ? y0 : y ;
ymax = ( y0 > y ) ? y0 : y ;
click = 0 ;
mode = 0 ;
redraw = 1 ;
continue ;
} else {
click = 0 ;
mode = 0 ;
redraw = 1 ;
continue ;
}
}
// Unzoom
if ( c = = ' r ' ) {
xmin = 0.0 ;
xmax = 360.0 ;
ymin = - 90.0 ;
ymax = 90.0 ;
mode = 0 ;
click = 0 ;
redraw = 1 ;
continue ;
}
// Default tle
if ( c = = ' t ' ) {
2013-10-24 00:00:12 -06:00
orb . satno = d . p [ 0 ] . satno ;
strcpy ( orb . desig , d . p [ 0 ] . desig ) ;
2013-05-18 11:54:11 -06:00
orb . ep_day = mjd2doy ( 0.5 * ( mjdmin + mjdmax ) , & orb . ep_year ) ;
2013-10-24 00:00:12 -06:00
satno = orb . satno ;
2013-07-10 03:54:14 -06:00
if ( elset = = 0 ) {
orb . eqinc = 0.5 * M_PI ;
orb . ascn = 0.0 ;
2014-09-17 08:42:16 -06:00
orb . ecc = 0.0001 ;
2013-07-10 03:54:14 -06:00
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 ;
2014-09-17 08:42:16 -06:00
orb . ecc = 0.0001 ;
2013-07-10 03:54:14 -06:00
orb . argp = 0.0 ;
orb . mnan = 0.0 ;
orb . rev = 1.0027 ;
orb . bstar = 0.0 ;
printf ( " GSO orbit \n " ) ;
} else if ( elset = = 3 ) {
2013-11-11 14:20:17 -07:00
orb . eqinc = 63.434 * D2R ;
2013-07-10 03:54:14 -06:00
orb . ascn = 0.0 ;
2013-11-11 14:20:17 -07:00
orb . ecc = 0.71 ;
2014-02-23 04:14:38 -07:00
orb . argp = 270 * D2R ;
orb . mnan = 0.0 ;
2013-11-11 14:20:17 -07:00
orb . rev = 2.006 ;
2013-07-10 03:54:14 -06:00
orb . bstar = 0.0 ;
printf ( " HEO orbit \n " ) ;
}
elset + + ;
if ( elset > 3 )
elset = 0 ;
2013-05-18 11:54:11 -06:00
print_orb ( & orb ) ;
printf ( " \n ================================================================================ \n " ) ;
click = 0 ;
redraw = 1 ;
continue ;
}
// Save
x0 = x ;
y0 = y ;
}
cpgend ( ) ;
2015-06-08 22:55:38 -06:00
if ( d . n ) free ( d . p ) ;
if ( d . m ) free ( d . q ) ;
2013-05-18 11:54:11 -06:00
fclose ( stderr ) ;
return 0 ;
}
// Get observing site
struct site 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 ] ;
struct site s ;
char * env , filename [ LIM ] ;
env = getenv ( " ST_DATADIR " ) ;
sprintf ( filename , " %s/data/sites.txt " , env ) ;
file = fopen ( filename , " r " ) ;
if ( file = = NULL ) {
printf ( " File with site information not found! \n " ) ;
2018-02-27 14:36:04 -07:00
return s ;
2013-05-18 11:54:11 -06:00
}
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 ;
// Copy site
if ( id = = site_id ) {
s . lat = lat ;
s . lng = lng ;
s . alt = alt ;
s . id = id ;
strcpy ( s . observer , observer ) ;
}
}
fclose ( file ) ;
2016-07-30 07:36:44 -06:00
if ( id ! = site_id )
s . id = = - 1 ;
2013-05-18 11:54:11 -06:00
return s ;
}
// 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 ;
}
// 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 ;
}
// Observer position
void obspos_xyz ( double mjd , double lng , double lat , float alt , xyz_t * pos , xyz_t * vel )
{
double ff , gc , gs , theta , s , dtheta ;
s = sin ( lat * D2R ) ;
ff = sqrt ( 1.0 - FLAT * ( 2.0 - FLAT ) * s * s ) ;
gc = 1.0 / ff + alt / XKMPER ;
gs = ( 1.0 - FLAT ) * ( 1.0 - FLAT ) / ff + alt / XKMPER ;
theta = gmst ( mjd ) + lng ;
dtheta = dgmst ( mjd ) * D2R / 86400 ;
pos - > x = gc * cos ( lat * D2R ) * cos ( theta * D2R ) * XKMPER ;
pos - > y = gc * cos ( lat * D2R ) * sin ( theta * D2R ) * XKMPER ;
pos - > z = gs * sin ( lat * D2R ) * XKMPER ;
vel - > x = - gc * cos ( lat * D2R ) * sin ( theta * D2R ) * XKMPER * dtheta ;
vel - > y = gc * cos ( lat * D2R ) * cos ( theta * D2R ) * XKMPER * dtheta ;
vel - > z = 0.0 ;
return ;
}
// Precess a celestial position
void precess ( double mjd0 , double ra0 , double de0 , double mjd , double * ra , double * de )
{
double t0 , t ;
double zeta , z , theta ;
double a , b , c ;
// Angles in radians
ra0 * = D2R ;
de0 * = D2R ;
// Time in centuries
t0 = ( mjd0 - 51544.5 ) / 36525.0 ;
t = ( mjd - mjd0 ) / 36525.0 ;
// Precession angles
zeta = ( 2306.2181 + 1.39656 * t0 - 0.000139 * t0 * t0 ) * t ;
zeta + = ( 0.30188 - 0.000344 * t0 ) * t * t + 0.017998 * t * t * t ;
zeta * = D2R / 3600.0 ;
z = ( 2306.2181 + 1.39656 * t0 - 0.000139 * t0 * t0 ) * t ;
z + = ( 1.09468 + 0.000066 * t0 ) * t * t + 0.018203 * t * t * t ;
z * = D2R / 3600.0 ;
theta = ( 2004.3109 - 0.85330 * t0 - 0.000217 * t0 * t0 ) * t ;
theta + = - ( 0.42665 + 0.000217 * t0 ) * t * t - 0.041833 * t * t * t ;
theta * = D2R / 3600.0 ;
a = cos ( de0 ) * sin ( ra0 + zeta ) ;
b = cos ( theta ) * cos ( de0 ) * cos ( ra0 + zeta ) - sin ( theta ) * sin ( de0 ) ;
c = sin ( theta ) * cos ( de0 ) * cos ( ra0 + zeta ) + cos ( theta ) * sin ( de0 ) ;
* ra = ( atan2 ( a , b ) + z ) * R2D ;
* de = asin ( c ) * R2D ;
if ( * ra < 360.0 )
* ra + = 360.0 ;
if ( * ra > 360.0 )
* ra - = 360.0 ;
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 ;
2015-04-16 01:30:49 -06:00
if ( year = = 1582 & & month = = 10 & & day < = 4 ) b = 0 ;
2013-05-18 11:54:11 -06:00
jd = floor ( 365.25 * ( year + 4716 ) ) + floor ( 30.6001 * ( month + 1 ) ) + day + b - 1524.5 ;
return jd - 2400000.5 ;
}
// Decode IOD Observations
struct point decode_iod_observation ( char * iod_line )
{
int year , month , iday , hour , min ;
int format , epoch , me , xe , sign ;
int site_id ;
double sec , ra , mm , ss , de , dd , ds , day , mjd0 ;
2013-10-24 00:00:12 -06:00
char secbuf [ 6 ] , sn [ 2 ] , degbuf [ 3 ] , buf1 [ 3 ] , buf2 [ 6 ] ;
2013-05-18 11:54:11 -06:00
struct point p ;
struct site s ;
xyz_t vel ;
// Strip newline
iod_line [ strlen ( iod_line ) - 1 ] = ' \0 ' ;
// Copy full line
strcpy ( p . iod_line , iod_line ) ;
// Set flag
p . flag = 1 ;
// Get SSN
sscanf ( iod_line , " %5d " , & p . satno ) ;
2013-10-24 00:00:12 -06:00
// Get desig
sscanf ( iod_line + 6 , " %s %s " , buf1 , buf2 ) ;
sprintf ( p . desig , " %s%s " , buf1 , buf2 ) ;
2013-05-18 11:54:11 -06:00
// Get site
sscanf ( iod_line + 16 , " %4d " , & site_id ) ;
s = get_site ( site_id ) ;
2016-07-30 07:36:44 -06:00
// Skip if site not found
if ( s . id < 0 ) {
fprintf ( stderr , " Site %d not found! \n " , site_id ) ;
p . flag = 0 ;
}
2013-05-18 11:54:11 -06:00
// Decode date/time
sscanf ( iod_line + 23 , " %4d%2d%2d%2d%2d%5s " , & year , & month , & iday , & hour , & min , secbuf ) ;
sec = atof ( secbuf ) ;
sec / = pow ( 10 , strlen ( secbuf ) - 2 ) ;
day = ( double ) iday + ( double ) hour / 24.0 + ( double ) min / 1440.0 + ( double ) sec / 86400.0 ;
p . mjd = date2mjd ( year , month , day ) ;
// Get uncertainty in time
sscanf ( iod_line + 41 , " %1d%1d " , & me , & xe ) ;
p . st = ( float ) me * pow ( 10 , xe - 8 ) ;
// Get observer position
obspos_xyz ( p . mjd , s . lng , s . lat , s . alt , & p . obspos , & vel ) ;
// Skip empty observations
if ( strlen ( iod_line ) < 64 | | ( iod_line [ 54 ] ! = ' + ' & & iod_line [ 54 ] ! = ' - ' ) )
p . flag = 0 ;
// Get format, epoch
sscanf ( iod_line + 44 , " %1d%1d " , & format , & epoch ) ;
// Read position
sscanf ( iod_line + 47 , " %2lf%2lf%3lf%1s " , & ra , & mm , & ss , sn ) ;
sscanf ( iod_line + 55 , " %2lf%2lf%2s " , & de , & dd , degbuf ) ;
ds = atof ( degbuf ) ;
if ( strlen ( degbuf ) = = 1 )
ds * = 10 ;
sign = ( sn [ 0 ] = = ' - ' ) ? - 1 : 1 ;
sscanf ( iod_line + 62 , " %1d%1d " , & me , & xe ) ;
p . sr = ( float ) me * pow ( 10 , xe - 8 ) ;
// Decode position
switch ( format )
{
// Format 1: RA/DEC = HHMMSSs+DDMMSS MX (MX in seconds of arc)
case 1 :
ra + = mm / 60 + ss / 36000 ;
de = sign * ( de + dd / 60 + ds / 3600 ) ;
p . sr / = 3600.0 ;
break ;
// Format 2: RA/DEC = HHMMmmm+DDMMmm MX (MX in minutes of arc)
case 2 :
ra + = mm / 60 + ss / 60000 ;
de = sign * ( de + dd / 60 + ds / 6000 ) ;
p . sr / = 60.0 ;
break ;
// Format 3: RA/DEC = HHMMmmm+DDdddd MX (MX in degrees of arc)
case 3 :
ra + = mm / 60 + ss / 60000 ;
de = sign * ( de + dd / 100 + ds / 10000 ) ;
break ;
// Format 7: RA/DEC = HHMMSSs+DDdddd MX (MX in degrees of arc)
case 7 :
ra + = mm / 60 + ss / 36000 ;
de = sign * ( de + dd / 100 + ds / 10000 ) ;
break ;
default :
printf ( " IOD Format not implemented \n " ) ;
p . flag = 0 ;
break ;
}
// Convert to degrees
ra * = 15.0 ;
// Get precession epoch
if ( epoch = = 0 ) {
p . ra = ra ;
p . de = de ;
return p ;
} else if ( epoch = = 4 ) {
mjd0 = 33281.9235 ;
} else if ( epoch = = 5 ) {
mjd0 = 51544.5 ;
} else {
printf ( " Observing epoch not implemented \n " ) ;
p . flag = 0 ;
}
// Precess position
precess ( mjd0 , ra , de , p . mjd , & p . ra , & p . de ) ;
return p ;
}
2015-06-08 22:55:38 -06:00
// Decode doppler Observations
struct doppler decode_doppler_observation ( char * line )
{
struct doppler q ;
struct site s ;
float flux ;
int site_id ;
// Read line
sscanf ( line , " %lf %lf %f %d " , & q . mjd , & q . freq , & flux , & site_id ) ;
// Get site
s = get_site ( site_id ) ;
// Get observer position
obspos_xyz ( q . mjd , s . lng , s . lat , s . alt , & q . obspos , & q . obsvel ) ;
return q ;
}
2013-05-18 11:54:11 -06:00
// 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 = = ' \t ' )
c = ' ' ;
if ( c = = ' \n ' )
s [ i + + ] = c ;
s [ i ] = ' \0 ' ;
return i ;
}
2015-06-08 22:55:38 -06:00
// Read IOD data
struct data read_data ( char * iodfname , char * dopfname )
2013-05-18 11:54:11 -06:00
{
int i = 0 ;
char line [ LIM ] ;
FILE * file ;
struct data d ;
2015-06-08 22:55:38 -06:00
d . n = 0 ;
d . m = 0 ;
2013-05-18 11:54:11 -06:00
2015-06-08 22:55:38 -06:00
// Read IOD data
if ( iodfname ! = NULL ) {
// Open file
file = fopen ( iodfname , " r " ) ;
if ( file = = NULL ) {
fprintf ( stderr , " Failed to open %s \n " , iodfname ) ;
exit ( 1 ) ;
}
// Count lines
while ( fgetline ( file , line , LIM ) > 0 )
i + + ;
d . n = i ;
2013-05-18 11:54:11 -06:00
2015-06-08 22:55:38 -06:00
// Allocate
d . p = ( struct point * ) malloc ( sizeof ( struct point ) * d . n ) ;
// Rewind file
rewind ( file ) ;
// Read data
i = 0 ;
while ( fgetline ( file , line , LIM ) > 0 )
d . p [ i + + ] = decode_iod_observation ( line ) ;
// Close file
fclose ( file ) ;
}
2013-05-18 11:54:11 -06:00
2015-06-08 22:55:38 -06:00
// Read doppler data
if ( dopfname ! = NULL ) {
// Open file
file = fopen ( iodfname , " r " ) ;
if ( file = = NULL ) {
fprintf ( stderr , " Failed to open %s \n " , dopfname ) ;
exit ( 1 ) ;
}
// Count lines
i = 0 ;
while ( fgetline ( file , line , LIM ) > 0 )
i + + ;
d . m = i ;
// Allocate
d . q = ( struct doppler * ) malloc ( sizeof ( struct doppler ) * d . m ) ;
2013-05-18 11:54:11 -06:00
2015-06-08 22:55:38 -06:00
// Rewind file
rewind ( file ) ;
// Read data
i = 0 ;
while ( fgetline ( file , line , LIM ) > 0 )
d . q [ i + + ] = decode_doppler_observation ( line ) ;
// Close file
fclose ( file ) ;
}
2013-05-18 11:54:11 -06:00
return d ;
}
// Chi-squared
double chisq ( double a [ ] )
{
int i , imode , nsel ;
double chisq , rms ;
xyz_t satpos , satvel ;
double dx , dy , dz ;
double r ;
// Construct struct
// a[0]: inclination
// a[1]: RA of ascending node
// a[2]: eccentricity
// a[3]: argument of periastron
// a[4]: mean anomaly
// a[5]: revs per day
if ( a [ 2 ] < 0.0 )
a [ 2 ] = 0.0 ;
if ( a [ 0 ] < 0.0 ) {
a [ 0 ] * = - 1 ;
a [ 1 ] + = 180.0 ;
} else if ( a [ 0 ] > 180.0 ) {
a [ 0 ] = 180.0 ;
}
2013-10-01 10:53:57 -06:00
if ( a [ 5 ] > 20.0 )
a [ 5 ] = 20.0 ;
2014-11-09 07:41:30 -07:00
if ( a [ 5 ] < 0.05 )
a [ 5 ] = 0.05 ;
2013-05-18 11:54:11 -06:00
// Set parameters
orb . eqinc = RAD ( a [ 0 ] ) ;
orb . ascn = RAD ( modulo ( a [ 1 ] , 360.0 ) ) ;
orb . ecc = a [ 2 ] ;
orb . argp = RAD ( modulo ( a [ 3 ] , 360.0 ) ) ;
orb . mnan = RAD ( modulo ( a [ 4 ] , 360.0 ) ) ;
orb . rev = a [ 5 ] ;
orb . bstar = a [ 6 ] ;
// Initialize
imode = init_sgdp4 ( & orb ) ;
if ( imode = = SGDP4_ERROR )
printf ( " Error \n " ) ;
// Loop over points
for ( i = 0 , nsel = 0 , chisq = 0.0 , rms = 0.0 ; i < d . n ; i + + ) {
// Get satellite position
satpos_xyz ( d . p [ i ] . mjd + 2400000.5 , & satpos , & satvel ) ;
// compute difference vector
dx = satpos . x - d . p [ i ] . obspos . x ;
dy = satpos . y - d . p [ i ] . obspos . y ;
dz = satpos . z - d . p [ i ] . obspos . z ;
// Celestial position
r = sqrt ( dx * dx + dy * dy + dz * dz ) ;
d . p [ i ] . rac = modulo ( atan2 ( dy , dx ) * R2D , 360.0 ) ;
d . p [ i ] . dec = asin ( dz / r ) * R2D ;
2018-03-01 09:15:04 -07:00
2013-05-18 11:54:11 -06:00
// Compute offset
forward ( d . p [ i ] . ra , d . p [ i ] . de , d . p [ i ] . rac , d . p [ i ] . dec , & d . p [ i ] . dx , & d . p [ i ] . dy ) ;
d . p [ i ] . dr = sqrt ( d . p [ i ] . dx * d . p [ i ] . dx + d . p [ i ] . dy * d . p [ i ] . dy ) ;
if ( d . p [ i ] . flag = = 2 ) {
// Compute chi-squared
chisq + = pow ( d . p [ i ] . dr / d . p [ i ] . sr , 2 ) ;
// Compute rms
rms + = pow ( d . p [ i ] . dr , 2 ) ;
// Count selected points
nsel + + ;
}
}
if ( nsel > 0 )
rms = sqrt ( rms / ( float ) nsel ) ;
2015-04-16 01:30:49 -06:00
2013-05-18 11:54:11 -06:00
d . chisq = chisq ;
d . rms = rms ;
d . nsel = nsel ;
return chisq ;
}
// Read tle
orbit_t read_tle ( char * filename , int satno )
{
int i ;
FILE * file ;
orbit_t orb ;
file = fopen ( filename , " r " ) ;
if ( file = = NULL )
fatal_error ( " Failed to open %s \n " , filename ) ;
// Read TLE
read_twoline ( file , satno , & orb ) ;
fclose ( file ) ;
return orb ;
}
// MJD to DOY
double mjd2doy ( double mjd , int * yr )
{
int year , month , k = 2 ;
double day , doy ;
mjd2date ( mjd , & year , & month , & day ) ;
if ( year % 4 = = 0 & & year % 400 ! = 0 )
k = 1 ;
doy = floor ( 275.0 * month / 9.0 ) - k * floor ( ( month + 9.0 ) / 12.0 ) + day - 30 ;
* yr = year ;
return doy ;
}
// Compute Date from Julian Day
void mjd2date ( double mjd , int * year , int * month , double * day )
{
double f , jd ;
int z , alpha , a , b , c , d , e ;
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 ) ;
* day = 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 ;
return ;
}
// Format TLE
void format_tle ( orbit_t orb , char * line1 , char * line2 )
{
int i , csum ;
char sbstar [ ] = " 00000-0 " , bstar [ 13 ] ;
// Format Bstar term
if ( fabs ( orb . bstar ) > 1e-9 ) {
sprintf ( bstar , " %11.4e " , 10 * orb . bstar ) ;
sbstar [ 0 ] = bstar [ 0 ] ; sbstar [ 1 ] = bstar [ 1 ] ; sbstar [ 2 ] = bstar [ 3 ] ; sbstar [ 3 ] = bstar [ 4 ] ;
sbstar [ 4 ] = bstar [ 5 ] ; sbstar [ 5 ] = bstar [ 6 ] ; sbstar [ 6 ] = bstar [ 8 ] ; sbstar [ 7 ] = bstar [ 10 ] ; sbstar [ 8 ] = ' \0 ' ;
}
// Print lines
2013-10-24 00:00:12 -06:00
sprintf ( line1 , " 1 %05dU %-8s %2d%012.8f .00000000 00000-0 %8s 0 0 " , orb . satno , orb . desig , orb . ep_year - 2000 , orb . ep_day , sbstar ) ;
2013-05-18 11:54:11 -06:00
sprintf ( line2 , " 2 %05d %8.4f %8.4f %07.0f %8.4f %8.4f %11.8f 0 " , orb . satno , DEG ( orb . eqinc ) , DEG ( orb . ascn ) , 1E7 * orb . ecc , DEG ( orb . argp ) , DEG ( orb . mnan ) , orb . rev ) ;
// Compute checksums
for ( i = 0 , csum = 0 ; i < strlen ( line1 ) ; i + + ) {
if ( isdigit ( line1 [ i ] ) )
csum + = line1 [ i ] - ' 0 ' ;
else if ( line1 [ i ] = = ' - ' )
csum + + ;
}
sprintf ( line1 , " %s%d " , line1 , csum % 10 ) ;
for ( i = 0 , csum = 0 ; i < strlen ( line2 ) ; i + + ) {
if ( isdigit ( line2 [ i ] ) )
csum + = line2 [ i ] - ' 0 ' ;
else if ( line2 [ i ] = = ' - ' )
csum + + ;
}
sprintf ( line2 , " %s%d " , line2 , csum % 10 ) ;
return ;
}
// Highlight
void highlight ( float x0 , float y0 , float x , float y , int flag )
{
int i ;
float xmin , xmax , ymin , ymax ;
xmin = ( x0 < x ) ? x0 : x ;
xmax = ( x0 > x ) ? x0 : x ;
ymin = ( y0 < y ) ? y0 : y ;
ymax = ( y0 > y ) ? y0 : y ;
for ( i = 0 ; i < d . n ; i + + )
if ( d . p [ i ] . ra > xmin & & d . p [ i ] . ra < xmax & & d . p [ i ] . de > ymin & & d . p [ i ] . de < ymax & & d . p [ i ] . flag ! = 0 )
d . p [ i ] . flag = flag ;
return ;
}
// Select time range
void time_range ( double * mjdmin , double * mjdmax , int flag )
{
int i , n ;
float c ;
for ( i = 0 , n = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . flag = = flag ) {
if ( n = = 0 ) {
* mjdmin = d . p [ i ] . mjd ;
* mjdmax = d . p [ i ] . mjd ;
}
if ( d . p [ i ] . mjd < * mjdmin ) * mjdmin = d . p [ i ] . mjd ;
if ( d . p [ i ] . mjd > * mjdmax ) * mjdmax = d . p [ i ] . mjd ;
n + + ;
}
}
c = 0.1 * ( * mjdmax - * mjdmin ) ;
* mjdmin - = c ;
* mjdmax + = c ;
return ;
}
// Print TLE
void print_tle ( orbit_t orb , char * filename )
{
int i , n ;
FILE * file ;
double mjdmin , mjdmax ;
int year , month ;
double day ;
char line1 [ 70 ] , line2 [ 70 ] ;
// Count number of points
for ( i = 0 , n = 0 ; i < d . n ; i + + ) {
if ( d . p [ i ] . flag = = 2 ) {
if ( n = = 0 ) {
mjdmin = d . p [ i ] . mjd ;
mjdmax = d . p [ i ] . mjd ;
}
if ( d . p [ i ] . mjd < mjdmin ) mjdmin = d . p [ i ] . mjd ;
if ( d . p [ i ] . mjd > mjdmax ) mjdmax = d . p [ i ] . mjd ;
n + + ;
}
}
// Write TLE
2014-09-24 01:27:45 -06:00
file = fopen ( filename , " a " ) ;
2013-05-18 11:54:11 -06:00
format_tle ( orb , line1 , line2 ) ;
fprintf ( file , " OBJ \n %s \n %s \n " , line1 , line2 ) ;
mjd2date ( mjdmin , & year , & month , & day ) ;
fprintf ( file , " # %4d%02d%05.2lf- " , year , month , day ) ;
mjd2date ( mjdmax , & year , & month , & day ) ;
2013-10-24 00:00:12 -06:00
fprintf ( file , " %4d%02d%05.2lf, %d measurements, %.3lf deg rms \n " , year , month , day , n , d . rms ) ;
2013-05-18 11:54:11 -06:00
fclose ( file ) ;
return ;
}
// Fit
void fit ( orbit_t orb , int * ia )
{
int i , n ;
double a [ 7 ] , da [ 7 ] ;
2013-08-31 01:49:06 -06:00
// double db[7]={5.0,5.0,0.1,5.0,5.0,0.5,0.0001};
2014-03-26 07:51:34 -06:00
// double db[7]={1.0,1.0,0.02,1.0,1.0,0.1,0.0001};
double db [ 7 ] = { 0.1 , 0.1 , 0.002 , 0.1 , 0.1 , 0.01 , 0.00001 } ;
2013-05-18 11:54:11 -06:00
a [ 0 ] = orb . eqinc * R2D ;
2013-07-10 03:54:14 -06:00
da [ 0 ] = da [ 0 ] * R2D ;
2013-05-18 11:54:11 -06:00
a [ 1 ] = orb . ascn * R2D ;
2013-07-10 03:54:14 -06:00
da [ 1 ] = da [ 1 ] * R2D ;
2013-05-18 11:54:11 -06:00
a [ 2 ] = orb . ecc ;
a [ 3 ] = orb . argp * R2D ;
2013-07-10 03:54:14 -06:00
da [ 3 ] = da [ 3 ] * R2D ;
2013-05-18 11:54:11 -06:00
a [ 4 ] = orb . mnan * R2D ;
2013-07-10 03:54:14 -06:00
da [ 4 ] = da [ 4 ] * R2D ;
2013-05-18 11:54:11 -06:00
a [ 5 ] = orb . rev ;
a [ 6 ] = orb . bstar ;
for ( i = 0 ; i < 7 ; i + + ) {
if ( ia [ i ] = = 1 )
da [ i ] = db [ i ] ;
else
da [ i ] = 0.0 ;
}
// Construct struct
// a[0]: inclination
// a[1]: RA of ascending node
// a[2]: eccentricity
// a[3]: argument of periastron
// a[4]: mean anomaly
// a[5]: revs per day
// a[6]: bstar
// Count highlighted points
for ( i = 0 , n = 0 ; i < d . n ; i + + )
if ( d . p [ i ] . flag = = 2 )
n + + ;
if ( n > 0 )
versafit ( n , 7 , a , da , chisq , 0.0 , 1e-6 , " n " ) ;
// Return parameters
orb . eqinc = RAD ( a [ 0 ] ) ;
orb . ascn = RAD ( modulo ( a [ 1 ] , 360.0 ) ) ;
orb . ecc = a [ 2 ] ;
orb . argp = RAD ( modulo ( a [ 3 ] , 360.0 ) ) ;
orb . mnan = RAD ( modulo ( a [ 4 ] , 360.0 ) ) ;
orb . rev = a [ 5 ] ;
orb . bstar = a [ 6 ] ;
return ;
}
void usage ( )
{
2013-10-24 00:00:12 -06:00
printf ( " satfit d:c:i:haCo:p \n \n " ) ;
printf ( " d IOD observations \n " ) ;
printf ( " c TLE catalog \n " ) ;
printf ( " i Satellite ID (NORAD) \n " ) ;
printf ( " C Fit circular orbit \n " ) ;
printf ( " p No plotting \n " ) ;
printf ( " o Output TLE file \n " ) ;
printf ( " a Adjust MA and RAAN \n " ) ;
printf ( " h This help \n " ) ;
2013-05-18 11:54:11 -06:00
return ;
}
2014-09-24 01:27:45 -06:00
// 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 ;
}
// Dot product
double dot ( xyz_t a , xyz_t b )
{
return a . x * b . x + a . y * b . y + a . z * b . z ;
}
// Magnitude
double magnitude ( xyz_t a )
{
return sqrt ( dot ( a , a ) ) ;
}
// Cross product
xyz_t cross ( xyz_t a , xyz_t b )
{
xyz_t c ;
c . x = a . y * b . z - a . z * b . y ;
c . y = a . z * b . x - a . x * b . z ;
c . z = a . x * b . y - a . y * b . x ;
return c ;
}
// Clasical elements
orbit_t classel ( int ep_year , double ep_day , xyz_t r , xyz_t v )
{
int i ;
double rm , vm , vm2 , rvm , mu = 1.0 ; ;
double chi , xp , yp , sx , cx , b , ee ;
double a , ecc , incl , node , peri , mm , n ;
xyz_t h , e , kk , nn ;
orbit_t orb ;
r . x / = XKMPER ;
r . y / = XKMPER ;
r . z / = XKMPER ;
v . x / = ( XKE * XKMPER / AE * XMNPDA / 86400.0 ) ;
v . y / = ( XKE * XKMPER / AE * XMNPDA / 86400.0 ) ;
v . z / = ( XKE * XKMPER / AE * XMNPDA / 86400.0 ) ;
rm = magnitude ( r ) ;
vm2 = dot ( v , v ) ;
rvm = dot ( r , v ) ;
h = cross ( r , v ) ;
chi = dot ( h , h ) / mu ;
e . x = ( vm2 / mu - 1.0 / rm ) * r . x - rvm / mu * v . x ;
e . y = ( vm2 / mu - 1.0 / rm ) * r . y - rvm / mu * v . y ;
e . z = ( vm2 / mu - 1.0 / rm ) * r . z - rvm / mu * v . z ;
a = pow ( 2.0 / rm - vm2 / mu , - 1 ) ;
ecc = magnitude ( e ) ;
incl = acos ( h . z / magnitude ( h ) ) * R2D ;
kk . x = 0.0 ;
kk . y = 0.0 ;
kk . z = 1.0 ;
nn = cross ( kk , h ) ;
if ( nn . x = = 0.0 & & nn . y = = 0.0 )
nn . x = 1.0 ;
node = atan2 ( nn . y , nn . x ) * R2D ;
if ( node < 0.0 )
node + = 360.0 ;
peri = acos ( dot ( nn , e ) / ( magnitude ( nn ) * ecc ) ) * R2D ;
if ( e . z < 0.0 )
peri = 360.0 - peri ;
if ( peri < 0.0 )
peri + = 360.0 ;
// Elliptic motion
if ( ecc < 1.0 ) {
xp = ( chi - rm ) / ecc ;
yp = rvm / ecc * sqrt ( chi / mu ) ;
b = a * sqrt ( 1.0 - ecc * ecc ) ;
cx = xp / a + ecc ;
sx = yp / b ;
ee = atan2 ( sx , cx ) ;
n = XKE * sqrt ( mu / ( a * a * a ) ) ;
mm = ( ee - ecc * sx ) * R2D ;
}
if ( mm < 0.0 )
mm + = 360.0 ;
// Fill
orb . satno = 0 ;
orb . eqinc = incl * D2R ;
orb . ascn = node * D2R ;
orb . argp = peri * D2R ;
orb . mnan = mm * D2R ;
orb . ecc = ecc ;
orb . rev = XKE * pow ( a , - 3.0 / 2.0 ) * XMNPDA / ( 2.0 * M_PI ) ;
orb . bstar = 0.0 ;
orb . ep_year = ep_year ;
orb . ep_day = ep_day ;
orb . norb = 0 ;
return orb ;
}
// Position and velocity to elements
orbit_t rv2el ( int satno , double mjd , xyz_t r0 , xyz_t v0 )
{
int i , imode ;
orbit_t orb [ 5 ] , orb1 [ 5 ] ;
xyz_t r , v ;
kep_t kep ;
char line1 [ 70 ] , line2 [ 70 ] ;
int ep_year ;
double ep_day ;
// Epoch
ep_day = mjd2doy ( mjd , & ep_year ) ;
// Initial guess
orb [ 0 ] = classel ( ep_year , ep_day , r0 , v0 ) ;
orb [ 0 ] . satno = satno ;
for ( i = 0 ; i < 4 ; i + + ) {
// Propagate
imode = init_sgdp4 ( & orb [ i ] ) ;
imode = satpos_xyz ( mjd + 2400000.5 , & r , & v ) ;
// Compute initial orbital elements
orb1 [ i ] = classel ( ep_year , ep_day , r , v ) ;
// Adjust
orb [ i + 1 ] . rev = orb [ i ] . rev + orb [ 0 ] . rev - orb1 [ i ] . rev ;
orb [ i + 1 ] . ascn = orb [ i ] . ascn + orb [ 0 ] . ascn - orb1 [ i ] . ascn ;
orb [ i + 1 ] . argp = orb [ i ] . argp + orb [ 0 ] . argp - orb1 [ i ] . argp ;
orb [ i + 1 ] . mnan = orb [ i ] . mnan + orb [ 0 ] . mnan - orb1 [ i ] . mnan ;
orb [ i + 1 ] . eqinc = orb [ i ] . eqinc + orb [ 0 ] . eqinc - orb1 [ i ] . eqinc ;
orb [ i + 1 ] . ecc = orb [ i ] . ecc + orb [ 0 ] . ecc - orb1 [ i ] . ecc ;
orb [ i + 1 ] . ep_year = orb [ i ] . ep_year ;
orb [ i + 1 ] . ep_day = orb [ i ] . ep_day ;
orb [ i + 1 ] . satno = orb [ i ] . satno ;
orb [ i + 1 ] . norb = orb [ i ] . norb ;
orb [ i + 1 ] . bstar = orb [ i ] . bstar ;
// Keep in range
if ( orb [ i + 1 ] . ecc < 0.0 )
orb [ i + 1 ] . ecc = 0.0 ;
if ( orb [ i + 1 ] . eqinc < 0.0 )
orb [ i + 1 ] . eqinc = 0.0 ;
}
2015-04-16 01:30:49 -06:00
orb [ i ] . mnan = modulo ( orb [ i ] . mnan , 2.0 * M_PI ) ;
orb [ i ] . ascn = modulo ( orb [ i ] . ascn , 2.0 * M_PI ) ;
orb [ i ] . argp = modulo ( orb [ i ] . argp , 2.0 * M_PI ) ;
2014-09-24 01:27:45 -06:00
return orb [ i ] ;
}
2018-03-01 09:15:04 -07:00
// 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 , status ;
double phi , theta ;
struct celprm cel ;
// Initialize Reference Angles
celini ( & cel ) ;
cel . ref [ 0 ] = ra0 ;
cel . ref [ 1 ] = de0 ;
cel . ref [ 2 ] = 999. ;
cel . ref [ 3 ] = 999. ;
cel . flag = 0. ;
strcpy ( cel . prj . code , " STG " ) ;
if ( celset ( & cel ) ) {
printf ( " Error in Projection (celset) \n " ) ;
return ;
}
cels2x ( & cel , 1 , 0 , 1 , 1 , & ra , & de , & phi , & theta , x , y , & status ) ;
return ;
}