sgp_obs.c
1 /* 2 * Unit SGP_Obs 3 * Author: Dr TS Kelso 4 * Original Version: 1992 Jun 02 5 * Current Revision: 1992 Sep 28 6 * Version: 1.40 7 * Copyright: 1992, All Rights Reserved 8 * 9 * Ported to C by: Neoklis Kyriazis April 9 2001 10 */ 11 12 #include "sgp4sdp4.h" 13 14 /* Procedure Calculate_User_PosVel passes the user's geodetic position */ 15 /* and the time of interest and returns the ECI position and velocity */ 16 /* of the observer. The velocity calculation assumes the geodetic */ 17 /* position is stationary relative to the earth's surface. */ 18 void Calculate_User_PosVel(double _time, geodetic_t * geodetic, 19 vector_t * obs_pos, vector_t * obs_vel) 20 { 21 /* Reference: The 1992 Astronomical Almanac, page K11. */ 22 23 double c, sq, achcp; 24 25 geodetic->theta = FMod2p(ThetaG_JD(_time) + geodetic->lon); 26 /* LMST */ 27 c = 1 / sqrt(1 + __f * (__f - 2) * Sqr(sin(geodetic->lat))); 28 sq = Sqr(1 - __f) * c; 29 achcp = (xkmper * c + geodetic->alt) * cos(geodetic->lat); 30 obs_pos->x = achcp * cos(geodetic->theta); /* km */ 31 obs_pos->y = achcp * sin(geodetic->theta); 32 obs_pos->z = (xkmper * sq + geodetic->alt) * sin(geodetic->lat); 33 obs_vel->x = -mfactor * obs_pos->y; /* km/sec */ 34 obs_vel->y = mfactor * obs_pos->x; 35 obs_vel->z = 0; 36 Magnitude(obs_pos); 37 Magnitude(obs_vel); 38 } 39 40 /* Procedure Calculate_LatLonAlt will calculate the geodetic */ 41 /* position of an object given its ECI position pos and time. */ 42 /* It is intended to be used to determine the ground track of */ 43 /* a satellite. The calculations assume the earth to be an */ 44 /* oblate spheroid as defined in WGS '72. */ 45 void Calculate_LatLonAlt(double _time, vector_t * pos, geodetic_t * geodetic) 46 { 47 /* Reference: The 1992 Astronomical Almanac, page K12. */ 48 49 double r, e2, phi, c; 50 51 geodetic->theta = AcTan(pos->y, pos->x); /* rad */ 52 geodetic->lon = FMod2p(geodetic->theta - ThetaG_JD(_time)); /* rad */ 53 r = sqrt(Sqr(pos->x) + Sqr(pos->y)); 54 e2 = __f * (2 - __f); 55 geodetic->lat = AcTan(pos->z, r); /* rad */ 56 57 do 58 { 59 phi = geodetic->lat; 60 c = 1 / sqrt(1 - e2 * Sqr(sin(phi))); 61 geodetic->lat = AcTan(pos->z + xkmper * c * e2 * sin(phi), r); 62 } 63 while (fabs(geodetic->lat - phi) >= 1E-10); 64 65 geodetic->alt = r / cos(geodetic->lat) - xkmper * c; /* km */ 66 67 if (geodetic->lat > pio2) 68 geodetic->lat -= twopi; 69 } 70 71 /* The procedures Calculate_Obs and Calculate_RADec calculate */ 72 /* the *topocentric* coordinates of the object with ECI position, */ 73 /* {pos}, and velocity, {vel}, from location {geodetic} at {time}. */ 74 /* The {obs_set} returned for Calculate_Obs consists of azimuth, */ 75 /* elevation, range, and range rate (in that order) with units of */ 76 /* radians, radians, kilometers, and kilometers/second, respectively. */ 77 /* The WGS '72 geoid is used and the effect of atmospheric refraction */ 78 /* (under standard temperature and pressure) is incorporated into the */ 79 /* elevation calculation; the effect of atmospheric refraction on */ 80 /* range and range rate has not yet been quantified. */ 81 82 /* The {obs_set} for Calculate_RADec consists of right ascension and */ 83 /* declination (in that order) in radians. Again, calculations are */ 84 /* based on *topocentric* position using the WGS '72 geoid and */ 85 /* incorporating atmospheric refraction. */ 86 void Calculate_Obs(double _time, vector_t * pos, 87 vector_t * vel, geodetic_t * geodetic, obs_set_t * obs_set) 88 { 89 double sin_lat, cos_lat, sin_theta, cos_theta; 90 double el, azim, top_s, top_e, top_z; 91 92 vector_t obs_pos, obs_vel, range, rgvel; 93 94 Calculate_User_PosVel(_time, geodetic, &obs_pos, &obs_vel); 95 96 range.x = pos->x - obs_pos.x; 97 range.y = pos->y - obs_pos.y; 98 range.z = pos->z - obs_pos.z; 99 100 rgvel.x = vel->x - obs_vel.x; 101 rgvel.y = vel->y - obs_vel.y; 102 rgvel.z = vel->z - obs_vel.z; 103 104 Magnitude(&range); 105 106 sin_lat = sin(geodetic->lat); 107 cos_lat = cos(geodetic->lat); 108 sin_theta = sin(geodetic->theta); 109 cos_theta = cos(geodetic->theta); 110 top_s = sin_lat * cos_theta * range.x 111 + sin_lat * sin_theta * range.y - cos_lat * range.z; 112 top_e = -sin_theta * range.x + cos_theta * range.y; 113 top_z = cos_lat * cos_theta * range.x 114 + cos_lat * sin_theta * range.y + sin_lat * range.z; 115 azim = atan(-top_e / top_s); /*Azimuth */ 116 if (top_s > 0) 117 azim = azim + pi; 118 if (azim < 0) 119 azim = azim + twopi; 120 el = ArcSin(top_z / range.w); 121 obs_set->az = azim; /* Azimuth (radians) */ 122 obs_set->el = el; /* Elevation (radians) */ 123 obs_set->range = range.w; /* Range (kilometers) */ 124 125 /* Range Rate (kilometers/second) */ 126 obs_set->range_rate = Dot(&range, &rgvel) / range.w; 127 128 /* Corrections for atmospheric refraction */ 129 /* Reference: Astronomical Algorithms by Jean Meeus, pp. 101-104 */ 130 /* Correction is meaningless when apparent elevation is below horizon */ 131 // obs_set->el = obs_set->el + Radians((1.02/tan(Radians(Degrees(el)+ 132 // 10.3/(Degrees(el)+5.11))))/60); 133 if (obs_set->el >= 0) 134 SetFlag(VISIBLE_FLAG); 135 else 136 { 137 obs_set->el = el; /*Reset to true elevation */ 138 ClearFlag(VISIBLE_FLAG); 139 } 140 } 141 142 void Calculate_RADec_and_Obs(double _time, vector_t * pos, vector_t * vel, 143 geodetic_t * geodetic, obs_astro_t * obs_set) 144 { 145 /* Reference: Methods of Orbit Determination by */ 146 /* Pedro Ramon Escobal, pp. 401-402 */ 147 148 double phi, theta, sin_theta, cos_theta, sin_phi, cos_phi, 149 az, el, Lxh, Lyh, Lzh, Sx, Ex, Zx, Sy, Ey, Zy, Sz, Ez, Zz, 150 Lx, Ly, Lz, cos_delta, sin_alpha, cos_alpha; 151 152 obs_set_t obs; 153 154 Calculate_Obs(_time, pos, vel, geodetic, &obs); 155 156 az = obs.az; 157 el = obs.el; 158 phi = geodetic->lat; 159 theta = FMod2p(ThetaG_JD(_time) + geodetic->lon); 160 sin_theta = sin(theta); 161 cos_theta = cos(theta); 162 sin_phi = sin(phi); 163 cos_phi = cos(phi); 164 Lxh = -cos(az) * cos(el); 165 Lyh = sin(az) * cos(el); 166 Lzh = sin(el); 167 Sx = sin_phi * cos_theta; 168 Ex = -sin_theta; 169 Zx = cos_theta * cos_phi; 170 Sy = sin_phi * sin_theta; 171 Ey = cos_theta; 172 Zy = sin_theta * cos_phi; 173 Sz = -cos_phi; 174 Ez = 0; 175 Zz = sin_phi; 176 Lx = Sx * Lxh + Ex * Lyh + Zx * Lzh; 177 Ly = Sy * Lxh + Ey * Lyh + Zy * Lzh; 178 Lz = Sz * Lxh + Ez * Lyh + Zz * Lzh; 179 obs_set->dec = ArcSin(Lz); /* Declination (radians) */ 180 cos_delta = sqrt(1 - Sqr(Lz)); 181 sin_alpha = Ly / cos_delta; 182 cos_alpha = Lx / cos_delta; 183 obs_set->ra = AcTan(sin_alpha, cos_alpha); /* Right Ascension (radians) */ 184 obs_set->ra = FMod2p(obs_set->ra); 185 }