/[PAMELA software]/yodaUtility/sgp4/cEci.cpp
ViewVC logotype

Contents of /yodaUtility/sgp4/cEci.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (show annotations) (download)
Sat Jun 10 05:35:03 2006 UTC (18 years, 5 months ago) by kusanagi
Branch: MAIN
CVS Tags: yodaUtility2_0/00, yodaUtility2_2/00, yodaUtility2_1/00, HEAD
Changes since 1.1: +1 -1 lines
Just added "newline" at end of file to avoid compilation warnings.

1 //
2 // cEci.cpp
3 //
4 // Copyright (c) 2002-2003 Michael F. Henry
5 //
6 #include "stdafx.h"
7
8 #include "cEci.h"
9 #include "globals.h"
10
11 //////////////////////////////////////////////////////////////////////
12 // cEci Class
13 //////////////////////////////////////////////////////////////////////
14 cEci::cEci(const cVector &pos,
15 const cVector &vel,
16 const cJulian &date,
17 bool IsAeUnits /* = true */)
18 {
19 m_pos = pos;
20 m_vel = vel;
21 m_date = date;
22 m_VecUnits = (IsAeUnits ? UNITS_AE : UNITS_NONE);
23 }
24
25 //////////////////////////////////////////////////////////////////////
26 // cEci(cCoordGeo&, cJulian&)
27 // Calculate the ECI coordinates of the location "geo" at time "date".
28 // Assumes geo coordinates are km-based.
29 // Assumes the earth is an oblate spheroid as defined in WGS '72.
30 // Reference: The 1992 Astronomical Almanac, page K11
31 // Reference: www.celestrak.com (Dr. TS Kelso)
32 cEci::cEci(const cCoordGeo &geo, const cJulian &date)
33 {
34 m_VecUnits = UNITS_KM;
35
36 double mfactor = TWOPI * (OMEGA_E / SEC_PER_DAY);
37 double lat = geo.m_Lat;
38 double lon = geo.m_Lon;
39 double alt = geo.m_Alt;
40
41 // Calculate Local Mean Sidereal Time (theta)
42 double theta = date.toLMST(lon);
43 double c = 1.0 / sqrt(1.0 + F * (F - 2.0) * sqr(sin(lat)));
44 double s = sqr(1.0 - F) * c;
45 double achcp = (XKMPER_WGS72 * c + alt) * cos(lat);
46
47 m_date = date;
48
49 m_pos.m_x = achcp * cos(theta); // km
50 m_pos.m_y = achcp * sin(theta); // km
51 m_pos.m_z = (XKMPER_WGS72 * s + alt) * sin(lat); // km
52 m_pos.m_w = sqrt(sqr(m_pos.m_x) +
53 sqr(m_pos.m_y) +
54 sqr(m_pos.m_z)); // range, km
55
56 m_vel.m_x = -mfactor * m_pos.m_y; // km / sec
57 m_vel.m_y = mfactor * m_pos.m_x;
58 m_vel.m_z = 0.0;
59 m_vel.m_w = sqrt(sqr(m_vel.m_x) + // range rate km/sec^2
60 sqr(m_vel.m_y));
61 }
62
63 //////////////////////////////////////////////////////////////////////////////
64 // toGeo()
65 // Return the corresponding geodetic position (based on the current ECI
66 // coordinates/Julian date).
67 // Assumes the earth is an oblate spheroid as defined in WGS '72.
68 // Side effects: Converts the position and velocity vectors to km-based units.
69 // Reference: The 1992 Astronomical Almanac, page K12.
70 // Reference: www.celestrak.com (Dr. TS Kelso)
71 cCoordGeo cEci::toGeo()
72 {
73 ae2km(); // Vectors must be in kilometer-based units
74
75 double theta = AcTan(m_pos.m_y, m_pos.m_x);
76 double lon = fmod(theta - m_date.toGMST(), TWOPI);
77
78 if (lon < 0.0)
79 lon += TWOPI; // "wrap" negative modulo
80
81 double r = sqrt(sqr(m_pos.m_x) + sqr(m_pos.m_y));
82 double e2 = F * (2.0 - F);
83 double lat = AcTan(m_pos.m_z, r);
84
85 const double delta = 1.0e-07;
86 double phi;
87 double c;
88
89 do
90 {
91 phi = lat;
92 c = 1.0 / sqrt(1.0 - e2 * sqr(sin(phi)));
93 lat = AcTan(m_pos.m_z + XKMPER_WGS72 * c * e2 * sin(phi), r);
94 }
95 while (fabs(lat - phi) > delta);
96
97 double alt = r / cos(lat) - XKMPER_WGS72 * c;
98
99 return cCoordGeo(lat, lon, alt); // radians, radians, kilometers
100 }
101
102 //////////////////////////////////////////////////////////////////////////////
103 // ae2km()
104 // Convert the position and velocity vector units from AE-based units
105 // to kilometer based units.
106 void cEci::ae2km()
107 {
108 if (UnitsAreAe())
109 {
110 MulPos(XKMPER_WGS72 / AE); // km
111 MulVel((XKMPER_WGS72 / AE) * (MIN_PER_DAY / 86400)); // km/sec
112 m_VecUnits = UNITS_KM;
113 }
114 }

  ViewVC Help
Powered by ViewVC 1.1.23