sgp_obs.c

00001 /*
00002  * Unit SGP_Obs
00003  *           Author:  Dr TS Kelso 
00004  * Original Version:  1992 Jun 02 
00005  * Current Revision:  1992 Sep 28 
00006  *          Version:  1.40 
00007  *        Copyright:  1992, All Rights Reserved 
00008  *
00009  *   Ported to C by:  Neoklis Kyriazis  April 9 2001
00010  */
00011 
00012 #include "sgp4sdp4.h"
00013 
00014 /* Procedure Calculate_User_PosVel passes the user's geodetic position */
00015 /* and the time of interest and returns the ECI position and velocity  */
00016 /* of the observer. The velocity calculation assumes the geodetic      */
00017 /* position is stationary relative to the earth's surface.             */
00018 void
00019 Calculate_User_PosVel(double _time,
00020                       geodetic_t *geodetic,
00021                       vector_t *obs_pos,
00022                       vector_t *obs_vel)
00023 {
00024 /* Reference:  The 1992 Astronomical Almanac, page K11. */
00025 
00026   double c,sq,achcp;
00027 
00028   geodetic->theta = FMod2p(ThetaG_JD(_time) + geodetic->lon);/*LMST*/
00029   c = 1/sqrt(1 + __f*(__f - 2)*Sqr(sin(geodetic->lat)));
00030   sq = Sqr(1 - __f)*c;
00031   achcp = (xkmper*c + geodetic->alt)*cos(geodetic->lat);
00032   obs_pos->x = achcp*cos(geodetic->theta);/*kilometers*/
00033   obs_pos->y = achcp*sin(geodetic->theta);
00034   obs_pos->z = (xkmper*sq + geodetic->alt)*sin(geodetic->lat);
00035   obs_vel->x = -mfactor*obs_pos->y;/*kilometers/second*/
00036   obs_vel->y =  mfactor*obs_pos->x;
00037   obs_vel->z =  0;
00038   SgpMagnitude(obs_pos);
00039   SgpMagnitude(obs_vel);
00040 } /*Procedure Calculate_User_PosVel*/
00041 
00042 /*------------------------------------------------------------------*/
00043 
00044 /* Procedure Calculate_LatLonAlt will calculate the geodetic  */
00045 /* position of an object given its ECI position pos and time. */
00046 /* It is intended to be used to determine the ground track of */
00047 /* a satellite.  The calculations  assume the earth to be an  */
00048 /* oblate spheroid as defined in WGS '72.                     */
00049 void
00050 Calculate_LatLonAlt(double _time, vector_t *pos,  geodetic_t *geodetic)
00051 {
00052   /* Reference:  The 1992 Astronomical Almanac, page K12. */
00053 
00054   double r,e2,phi,c;
00055 
00056   geodetic->theta = AcTan(pos->y,pos->x);/*radians*/
00057   geodetic->lon = FMod2p(geodetic->theta - ThetaG_JD(_time));/*radians*/
00058   r = sqrt(Sqr(pos->x) + Sqr(pos->y));
00059   e2 = __f*(2 - __f);
00060   geodetic->lat = AcTan(pos->z,r);/*radians*/
00061 
00062   do
00063     {
00064       phi = geodetic->lat;
00065       c = 1/sqrt(1 - e2*Sqr(sin(phi)));
00066       geodetic->lat = AcTan(pos->z + xkmper*c*e2*sin(phi),r);
00067     }
00068   while(fabs(geodetic->lat - phi) >= 1E-10);
00069 
00070   geodetic->alt = r/cos(geodetic->lat) - xkmper*c;/*kilometers*/
00071 
00072   if( geodetic->lat > pio2 ) geodetic->lat -= twopi;
00073   
00074 } /*Procedure Calculate_LatLonAlt*/
00075 
00076 /*------------------------------------------------------------------*/
00077 
00078 /* The procedures Calculate_Obs and Calculate_RADec calculate         */
00079 /* the *topocentric* coordinates of the object with ECI position,     */
00080 /* {pos}, and velocity, {vel}, from location {geodetic} at {time}.    */
00081 /* The {obs_set} returned for Calculate_Obs consists of azimuth,      */
00082 /* elevation, range, and range rate (in that order) with units of     */
00083 /* radians, radians, kilometers, and kilometers/second, respectively. */
00084 /* The WGS '72 geoid is used and the effect of atmospheric refraction */
00085 /* (under standard temperature and pressure) is incorporated into the */
00086 /* elevation calculation; the effect of atmospheric refraction on     */
00087 /* range and range rate has not yet been quantified.                  */
00088 
00089 /* The {obs_set} for Calculate_RADec consists of right ascension and  */
00090 /* declination (in that order) in radians.  Again, calculations are   */
00091 /* based on *topocentric* position using the WGS '72 geoid and        */
00092 /* incorporating atmospheric refraction.                              */
00093 
00094 void
00095 Calculate_Obs(double _time,
00096                 vector_t *pos,
00097                 vector_t *vel,
00098                 geodetic_t *geodetic,
00099                 vector_t *obs_set)
00100   {
00101    double
00102      sin_lat,cos_lat,
00103      sin_theta,cos_theta,
00104      el,azim,
00105      top_s,top_e,top_z;
00106 
00107    vector_t
00108      obs_pos,obs_vel,range,rgvel;
00109 
00110   Calculate_User_PosVel(_time, geodetic, &obs_pos, &obs_vel);
00111 
00112     range.x = pos->x - obs_pos.x;
00113     range.y = pos->y - obs_pos.y;
00114     range.z = pos->z - obs_pos.z;
00115 
00116     rgvel.x = vel->x - obs_vel.x;
00117     rgvel.y = vel->y - obs_vel.y;
00118     rgvel.z = vel->z - obs_vel.z;
00119 
00120   SgpMagnitude(&range);
00121 
00122   sin_lat = sin(geodetic->lat);
00123   cos_lat = cos(geodetic->lat);
00124   sin_theta = sin(geodetic->theta);
00125   cos_theta = cos(geodetic->theta);
00126   top_s = sin_lat*cos_theta*range.x
00127          + sin_lat*sin_theta*range.y
00128          - cos_lat*range.z;
00129   top_e = -sin_theta*range.x
00130          + cos_theta*range.y;
00131   top_z = cos_lat*cos_theta*range.x
00132          + cos_lat*sin_theta*range.y
00133          + sin_lat*range.z;
00134   azim = atan(-top_e/top_s); /*Azimuth*/
00135   if( top_s > 0 ) 
00136     azim = azim + pi;
00137   if( azim < 0 )
00138     azim = azim + twopi;
00139   el = ArcSin(top_z/range.w);
00140   obs_set->x = azim;      /* Azimuth (radians)  */
00141   obs_set->y = el;        /* Elevation (radians)*/
00142   obs_set->z = range.w; /* Range (kilometers) */
00143 
00144  /*Range Rate (kilometers/second)*/
00145   obs_set->w = Dot(&range, &rgvel)/range.w;
00146 
00147 /* Corrections for atmospheric refraction */
00148 /* Reference:  Astronomical Algorithms by Jean Meeus, pp. 101-104    */
00149 /* Correction is meaningless when apparent elevation is below horizon */
00150   obs_set->y = obs_set->y + Radians((1.02/tan(Radians(Degrees(el)+
00151                10.3/(Degrees(el)+5.11))))/60);
00152   if( obs_set->y >= 0 )
00153     SetFlag(VISIBLE_FLAG);
00154   else
00155     {
00156     obs_set->y = el;  /*Reset to true elevation*/
00157     ClearFlag(VISIBLE_FLAG);
00158     } /*else*/
00159   } /*Procedure Calculate_Obs*/
00160 
00161 /*------------------------------------------------------------------*/
00162 
00163 void
00164 Calculate_RADec( double _time,
00165                  vector_t *pos,
00166                  vector_t *vel,
00167                  geodetic_t *geodetic,
00168                  vector_t *obs_set)
00169 {
00170 /* Reference:  Methods of Orbit Determination by  */
00171 /*                Pedro Ramon Escobal, pp. 401-402 */
00172 
00173 double
00174     phi,theta,sin_theta,cos_theta,sin_phi,cos_phi,
00175     az,el,Lxh,Lyh,Lzh,Sx,Ex,Zx,Sy,Ey,Zy,Sz,Ez,Zz,
00176     Lx,Ly,Lz,cos_delta,sin_alpha,cos_alpha;
00177 
00178   Calculate_Obs(_time,pos,vel,geodetic,obs_set);
00179 
00180 /*  if( isFlagSet(VISIBLE_FLAG) )
00181     {*/
00182     az = obs_set->x;
00183     el = obs_set->y;
00184     phi   = geodetic->lat;
00185     theta = FMod2p(ThetaG_JD(_time) + geodetic->lon);
00186     sin_theta = sin(theta);
00187     cos_theta = cos(theta);
00188     sin_phi = sin(phi);
00189     cos_phi = cos(phi);
00190     Lxh = -cos(az)*cos(el);
00191     Lyh =  sin(az)*cos(el);
00192     Lzh =  sin(el);
00193     Sx = sin_phi*cos_theta;
00194     Ex = -sin_theta;
00195     Zx = cos_theta*cos_phi;
00196     Sy = sin_phi*sin_theta;
00197     Ey = cos_theta;
00198     Zy = sin_theta*cos_phi;
00199     Sz = -cos_phi;
00200     Ez = 0;
00201     Zz = sin_phi;
00202     Lx = Sx*Lxh + Ex*Lyh + Zx*Lzh;
00203     Ly = Sy*Lxh + Ey*Lyh + Zy*Lzh;
00204     Lz = Sz*Lxh + Ez*Lyh + Zz*Lzh;
00205     obs_set->y = ArcSin(Lz);  /*Declination (radians)*/
00206     cos_delta = sqrt(1 - Sqr(Lz));
00207     sin_alpha = Ly/cos_delta;
00208     cos_alpha = Lx/cos_delta;
00209     obs_set->x = AcTan(sin_alpha,cos_alpha); /*Right Ascension (radians)*/
00210     obs_set->x = FMod2p(obs_set->x);
00211     /*}*/  /*if*/
00212   } /* Procedure Calculate_RADec */
00213 
00214 /*------------------------------------------------------------------*/