Andy K / Mbed 2 deprecated SOWB

Satellite Observers Workbench. NOT yet complete, just published for forum posters to \"cherry pick\" pieces of code as requiered as an example.

Dependencies:   mbed

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?

UserRevisionLine numberNew 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 /*------------------------------------------------------------------*/