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