Satellite Observers Workbench. NOT yet complete, just published for forum posters to \"cherry pick\" pieces of code as requiered as an example.
Diff: sgp4sdp4/sgp_obs.c
- Revision:
- 0:0a841b89d614
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sgp4sdp4/sgp_obs.c Mon Oct 11 10:34:55 2010 +0000 @@ -0,0 +1,214 @@ +/* + * Unit SGP_Obs + * Author: Dr TS Kelso + * Original Version: 1992 Jun 02 + * Current Revision: 1992 Sep 28 + * Version: 1.40 + * Copyright: 1992, All Rights Reserved + * + * Ported to C by: Neoklis Kyriazis April 9 2001 + */ + +#include "sgp4sdp4.h" + +/* Procedure Calculate_User_PosVel passes the user's geodetic position */ +/* and the time of interest and returns the ECI position and velocity */ +/* of the observer. The velocity calculation assumes the geodetic */ +/* position is stationary relative to the earth's surface. */ +void +Calculate_User_PosVel(double _time, + geodetic_t *geodetic, + vector_t *obs_pos, + vector_t *obs_vel) +{ +/* Reference: The 1992 Astronomical Almanac, page K11. */ + + double c,sq,achcp; + + geodetic->theta = FMod2p(ThetaG_JD(_time) + geodetic->lon);/*LMST*/ + c = 1/sqrt(1 + __f*(__f - 2)*Sqr(sin(geodetic->lat))); + sq = Sqr(1 - __f)*c; + achcp = (xkmper*c + geodetic->alt)*cos(geodetic->lat); + obs_pos->x = achcp*cos(geodetic->theta);/*kilometers*/ + obs_pos->y = achcp*sin(geodetic->theta); + obs_pos->z = (xkmper*sq + geodetic->alt)*sin(geodetic->lat); + obs_vel->x = -mfactor*obs_pos->y;/*kilometers/second*/ + obs_vel->y = mfactor*obs_pos->x; + obs_vel->z = 0; + SgpMagnitude(obs_pos); + SgpMagnitude(obs_vel); +} /*Procedure Calculate_User_PosVel*/ + +/*------------------------------------------------------------------*/ + +/* Procedure Calculate_LatLonAlt will calculate the geodetic */ +/* position of an object given its ECI position pos and time. */ +/* It is intended to be used to determine the ground track of */ +/* a satellite. The calculations assume the earth to be an */ +/* oblate spheroid as defined in WGS '72. */ +void +Calculate_LatLonAlt(double _time, vector_t *pos, geodetic_t *geodetic) +{ + /* Reference: The 1992 Astronomical Almanac, page K12. */ + + double r,e2,phi,c; + + geodetic->theta = AcTan(pos->y,pos->x);/*radians*/ + geodetic->lon = FMod2p(geodetic->theta - ThetaG_JD(_time));/*radians*/ + r = sqrt(Sqr(pos->x) + Sqr(pos->y)); + e2 = __f*(2 - __f); + geodetic->lat = AcTan(pos->z,r);/*radians*/ + + do + { + phi = geodetic->lat; + c = 1/sqrt(1 - e2*Sqr(sin(phi))); + geodetic->lat = AcTan(pos->z + xkmper*c*e2*sin(phi),r); + } + while(fabs(geodetic->lat - phi) >= 1E-10); + + geodetic->alt = r/cos(geodetic->lat) - xkmper*c;/*kilometers*/ + + if( geodetic->lat > pio2 ) geodetic->lat -= twopi; + +} /*Procedure Calculate_LatLonAlt*/ + +/*------------------------------------------------------------------*/ + +/* The procedures Calculate_Obs and Calculate_RADec calculate */ +/* the *topocentric* coordinates of the object with ECI position, */ +/* {pos}, and velocity, {vel}, from location {geodetic} at {time}. */ +/* The {obs_set} returned for Calculate_Obs consists of azimuth, */ +/* elevation, range, and range rate (in that order) with units of */ +/* radians, radians, kilometers, and kilometers/second, respectively. */ +/* The WGS '72 geoid is used and the effect of atmospheric refraction */ +/* (under standard temperature and pressure) is incorporated into the */ +/* elevation calculation; the effect of atmospheric refraction on */ +/* range and range rate has not yet been quantified. */ + +/* The {obs_set} for Calculate_RADec consists of right ascension and */ +/* declination (in that order) in radians. Again, calculations are */ +/* based on *topocentric* position using the WGS '72 geoid and */ +/* incorporating atmospheric refraction. */ + +void +Calculate_Obs(double _time, + vector_t *pos, + vector_t *vel, + geodetic_t *geodetic, + vector_t *obs_set) + { + double + sin_lat,cos_lat, + sin_theta,cos_theta, + el,azim, + top_s,top_e,top_z; + + vector_t + obs_pos,obs_vel,range,rgvel; + + Calculate_User_PosVel(_time, geodetic, &obs_pos, &obs_vel); + + range.x = pos->x - obs_pos.x; + range.y = pos->y - obs_pos.y; + range.z = pos->z - obs_pos.z; + + rgvel.x = vel->x - obs_vel.x; + rgvel.y = vel->y - obs_vel.y; + rgvel.z = vel->z - obs_vel.z; + + SgpMagnitude(&range); + + sin_lat = sin(geodetic->lat); + cos_lat = cos(geodetic->lat); + sin_theta = sin(geodetic->theta); + cos_theta = cos(geodetic->theta); + top_s = sin_lat*cos_theta*range.x + + sin_lat*sin_theta*range.y + - cos_lat*range.z; + top_e = -sin_theta*range.x + + cos_theta*range.y; + top_z = cos_lat*cos_theta*range.x + + cos_lat*sin_theta*range.y + + sin_lat*range.z; + azim = atan(-top_e/top_s); /*Azimuth*/ + if( top_s > 0 ) + azim = azim + pi; + if( azim < 0 ) + azim = azim + twopi; + el = ArcSin(top_z/range.w); + obs_set->x = azim; /* Azimuth (radians) */ + obs_set->y = el; /* Elevation (radians)*/ + obs_set->z = range.w; /* Range (kilometers) */ + + /*Range Rate (kilometers/second)*/ + obs_set->w = Dot(&range, &rgvel)/range.w; + +/* Corrections for atmospheric refraction */ +/* Reference: Astronomical Algorithms by Jean Meeus, pp. 101-104 */ +/* Correction is meaningless when apparent elevation is below horizon */ + obs_set->y = obs_set->y + Radians((1.02/tan(Radians(Degrees(el)+ + 10.3/(Degrees(el)+5.11))))/60); + if( obs_set->y >= 0 ) + SetFlag(VISIBLE_FLAG); + else + { + obs_set->y = el; /*Reset to true elevation*/ + ClearFlag(VISIBLE_FLAG); + } /*else*/ + } /*Procedure Calculate_Obs*/ + +/*------------------------------------------------------------------*/ + +void +Calculate_RADec( double _time, + vector_t *pos, + vector_t *vel, + geodetic_t *geodetic, + vector_t *obs_set) +{ +/* Reference: Methods of Orbit Determination by */ +/* Pedro Ramon Escobal, pp. 401-402 */ + +double + phi,theta,sin_theta,cos_theta,sin_phi,cos_phi, + az,el,Lxh,Lyh,Lzh,Sx,Ex,Zx,Sy,Ey,Zy,Sz,Ez,Zz, + Lx,Ly,Lz,cos_delta,sin_alpha,cos_alpha; + + Calculate_Obs(_time,pos,vel,geodetic,obs_set); + +/* if( isFlagSet(VISIBLE_FLAG) ) + {*/ + az = obs_set->x; + el = obs_set->y; + phi = geodetic->lat; + theta = FMod2p(ThetaG_JD(_time) + geodetic->lon); + sin_theta = sin(theta); + cos_theta = cos(theta); + sin_phi = sin(phi); + cos_phi = cos(phi); + Lxh = -cos(az)*cos(el); + Lyh = sin(az)*cos(el); + Lzh = sin(el); + Sx = sin_phi*cos_theta; + Ex = -sin_theta; + Zx = cos_theta*cos_phi; + Sy = sin_phi*sin_theta; + Ey = cos_theta; + Zy = sin_theta*cos_phi; + Sz = -cos_phi; + Ez = 0; + Zz = sin_phi; + Lx = Sx*Lxh + Ex*Lyh + Zx*Lzh; + Ly = Sy*Lxh + Ey*Lyh + Zy*Lzh; + Lz = Sz*Lxh + Ez*Lyh + Zz*Lzh; + obs_set->y = ArcSin(Lz); /*Declination (radians)*/ + cos_delta = sqrt(1 - Sqr(Lz)); + sin_alpha = Ly/cos_delta; + cos_alpha = Lx/cos_delta; + obs_set->x = AcTan(sin_alpha,cos_alpha); /*Right Ascension (radians)*/ + obs_set->x = FMod2p(obs_set->x); + /*}*/ /*if*/ + } /* Procedure Calculate_RADec */ + +/*------------------------------------------------------------------*/