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 |
}
|