Satellite Observers Workbench. NOT yet complete, just published for forum posters to \"cherry pick\" pieces of code as requiered as an example.
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 /*------------------------------------------------------------------*/
Generated on Tue Jul 12 2022 18:05:35 by 1.7.2