/* sat_id.cpp 8 March 2003 An example 'main' function illustrating how to find which satellite(s) are within a given radius of a given RA/dec, as seen from a given point. The code reads in a file of observations in MPC format (name provided as the first command-line argument). For example: sat_id mpc_obs.txt would hunt through the file 'mpc_obs.txt' for MPC-formatted observations. It would then read the file 'alldat.tle', looking for corresponding satellites within .2 degrees of said observations. It then spits out the original file, with satellite IDs added (when found) after each observation line. For each IDed satellite, the international and NORAD designations are given, along with its angular distance from the search point, position angle of motion, and apparent angular rate of motion in arcminutes/second (or, equivalently, degrees/minute). */ /* 2 July 2003: fixed the day/month/year to JD part of 'get_mpc_data()' so it will work for all years >= 0 (previously, it worked for years 2000 to 2099... plenty for the practical purpose of ID'ing recently-found satellites, but this is also an 'example' program.) */ #include #include #include #include #include "norad.h" #include "observe.h" #define PI 3.141592653589793238462643383279 #define TIME_EPSILON (1./86400.) static int get_mpc_data( const char *buff, double *jd, double *ra, double *dec) { int i1, i2, i, year, month; double tval, day; static const char month_len[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }; if( strlen( buff) < 80 || strlen( buff) > 82) return( -1); if( sscanf( buff + 32, "%d %d %lf", &i1, &i2, &tval) != 3) return( -2); *ra = ((double)i1 + (double)i2 / 60. + tval / 3600.) * (PI / 12.); if( sscanf( buff + 45, "%d %d %lf", &i1, &i2, &tval) != 3) return( -3); *dec = ((double)i1 + (double)i2 / 60. + tval / 3600.) * (PI / 180.); if( buff[44] == '-') *dec = -*dec; else if( buff[44] != '+') return( -4); /* Read in the day/month/year from the record... */ if( sscanf( buff + 15, "%d %d %lf", &year, &month, &day) != 3) return( -5); /* ...and convert to a JD: */ *jd = 1721059.5 + (double)( year * 365 + year / 4 - year / 100 + year / 400) + day; for( i = 0; i < month - 1; i++) *jd += (double)month_len[i]; if( month < 3 && !(year % 4)) /* leap years, January and February */ if( !(year % 400) || (year % 100)) (*jd)--; return( 0); } int main( int argc, char **argv) { const char *tle_file_name = "alldat.tle"; FILE *ifile = fopen( argv[1], "rb"), *tle_file; FILE *stations = fopen( "stations.txt", "rb"); char line1[100], line2[100], buff[90]; double search_radius = .2; /* default to .2-degree search */ double lon = 0., rho_sin_phi = 0., rho_cos_phi = 0.; char curr_mpc_code[4]; int i, debug_level = 0; if( !ifile) { printf( "Couldn't open input file %s\n", argv[1]); exit( -1); } curr_mpc_code[0] = curr_mpc_code[3] = '\0'; for( i = 1; i < argc; i++) if( argv[i][0] == '-') switch( argv[i][1]) { case 'r': search_radius = atof( argv[i] + 2); break; case 't': tle_file_name = argv[i] + 2; break; case 'd': debug_level = atoi( argv[i] + 2); break; default: printf( "Unrecognized command-line option '%s'\n", argv[i]); exit( -2); break; } tle_file = fopen( tle_file_name, "rb"); if( !tle_file) { printf( "Couldn't open TLE file %s\n", tle_file_name); exit( -1); } while( fgets( buff, sizeof( buff), ifile)) { double target_ra, target_dec, jd; printf( "%s", buff); if( !get_mpc_data( buff, &jd, &target_ra, &target_dec)) { double observer_loc[3], observer_loc2[3]; if( debug_level) printf( "JD = %.5lf; RA = %.5lf; dec = %.5lf\n", jd, target_ra * 180. / PI, target_dec * 180. / PI); if( memcmp( curr_mpc_code, buff + 77, 3)) { char tbuff[100]; int got_it = 0; memcpy( curr_mpc_code, buff + 77, 3); fseek( stations, 0L, SEEK_SET); while( !got_it && fgets( tbuff, sizeof( tbuff), stations)) got_it = !memcmp( tbuff, curr_mpc_code, 3); if( got_it) sscanf( tbuff + 3, "%lf %lf %lf", &lon, &rho_cos_phi, &rho_sin_phi); if( !got_it) printf( "FAILED to find MPC code %s\n", curr_mpc_code); } if( debug_level) printf( "lon = %.5lf rho cos phi = %.5lf rho sin phi = %.5lf\n", lon, rho_cos_phi, rho_sin_phi); observer_cartesian_coords( jd, lon * PI / 180., rho_cos_phi, rho_sin_phi, observer_loc); observer_cartesian_coords( jd + TIME_EPSILON, lon * PI / 180., rho_cos_phi, rho_sin_phi, observer_loc2); fseek( tle_file, 0L, SEEK_SET); fgets( line1, sizeof( line1), tle_file); while( fgets( line2, sizeof( line2), tle_file)) { tle_t tle; /* Structure for two-line elements set for satellite */ if( !parse_elements( line1, line2, &tle)) /* hey! we got a TLE! */ { int is_deep = select_ephemeris( &tle); double sat_params[N_SAT_PARAMS], radius, d_ra, d_dec; double ra, dec, dist_to_satellite, t_since; double pos[3]; /* Satellite position vector */ double unused_delta2; if( debug_level > 1) printf( "%s", line1); t_since = (jd - tle.epoch) * 1440.; if( is_deep) { SDP4_init( sat_params, &tle); SDP4( t_since, &tle, sat_params, pos, NULL); } else { SGP4_init( sat_params, &tle); SGP4( t_since, &tle, sat_params, pos, NULL); } if( debug_level > 1) printf( "%s", line2); if( debug_level > 2) printf( " %.5lf %.5lf %.5lf\n", pos[0], pos[1], pos[2]); get_satellite_ra_dec_delta( observer_loc, pos, &ra, &dec, &dist_to_satellite); if( debug_level > 3) printf( "RA: %.5lf dec: %.5lf\n", ra * 180. / PI, dec * 180. / PI); epoch_of_date_to_j2000( jd, &ra, &dec); d_ra = (ra - target_ra + PI * 4.); while( d_ra > PI) d_ra -= PI + PI; d_dec = dec - target_dec; radius = sqrt( d_ra * d_ra + d_dec * d_dec) * 180. / PI; if( radius < search_radius) /* good enough for us! */ { double speed, posn_ang_of_motion; line1[16] = '\0'; /* Compute position one second later, so we */ /* can show speed/PA of motion: */ t_since += TIME_EPSILON * 1440.; if( is_deep) SDP4( t_since, &tle, sat_params, pos, NULL); else SGP4( t_since, &tle, sat_params, pos, NULL); get_satellite_ra_dec_delta( observer_loc2, pos, &d_ra, &d_dec, &unused_delta2); epoch_of_date_to_j2000( jd, &d_ra, &d_dec); d_ra -= ra; d_dec -= dec; while( d_ra > PI) d_ra -= PI + PI; while( d_ra < -PI) d_ra += PI + PI; d_ra *= cos( dec); posn_ang_of_motion = atan2( d_ra, d_dec); if( posn_ang_of_motion < 0.) posn_ang_of_motion += PI + PI; speed = sqrt( d_ra * d_ra + d_dec * d_dec) * 180. / PI; /* Put RA into 0 to 2pi range: */ ra = fmod( ra + PI * 10., PI + PI); line1[8] = line1[16] = '\0'; memcpy( line1 + 30, line1 + 11, 6); line1[11] = '\0'; printf( " NORAD designation %s; international %s%s-%s\n", line1 + 2, (line1[9] >= '6' ? "19" : "20"), line1 + 9, line1 + 30); printf( " delta=%8.1lf km; offset=%5.2lf degrees; motion %5.2lf'/sec at PA=%3d\n", dist_to_satellite, radius, speed * 60., (int)(posn_ang_of_motion * 180 / PI)); /* "Speed" is displayed in arcminutes/second, or in degrees/minute */ } } strcpy( line1, line2); } } } fclose( tle_file); return( 0); } /* End of main() */