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

Annotation of /yodaUtility/sgp4/cSite.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Sun Apr 30 11:08:15 2006 UTC (18 years, 7 months ago) by kusanagi
Branch point for: MAIN
Initial revision

1 kusanagi 1.1 //
2     // cSite.cpp
3     //
4     // Copyright (c) 2003 Michael F. Henry
5     //
6     //////////////////////////////////////////////////////////////////////////////
7     #include "stdafx.h"
8     #include "cSite.h"
9     #include "globals.h"
10    
11     //////////////////////////////////////////////////////////////////////////////
12     // Construction/Destruction
13     cSite::cSite(const cCoordGeo &geo) : m_geo(geo)
14     {}
15    
16     //////////////////////////////////////////////////////////////////////////////
17     // c'tor accepting:
18     // Latitude in degress (negative south)
19     // Longitude in degress (negative west)
20     // Altitude in km
21     cSite::cSite(double degLat, double degLon, double kmAlt) :
22     m_geo(deg2rad(degLat), deg2rad(degLon), kmAlt)
23     {}
24    
25     cSite::~cSite()
26     {}
27    
28     //////////////////////////////////////////////////////////////////////////////
29     // setGeo()
30     // Set a new geographic position
31     void cSite::setGeo(const cCoordGeo &geo)
32     {
33     m_geo = geo;
34     }
35    
36     //////////////////////////////////////////////////////////////////////////////
37     // getPosition()
38     // Return the ECI coordinate of the site at the given time.
39     cEci cSite::getPosition(const cJulian &date) const
40     {
41     return cEci(m_geo, date);
42     }
43    
44     //////////////////////////////////////////////////////////////////////////////
45     // getLookAngle()
46     // Return the topocentric (azimuth, elevation, etc.) coordinates for a target
47     // object described by the input ECI coordinates.
48     cCoordTopo cSite::getLookAngle(const cEci &eci) const
49     {
50     // Calculate the ECI coordinates for this cSite object at the time
51     // of interest.
52     cJulian date = eci.getDate();
53     cEci eciSite(m_geo, date);
54    
55     // The Site ECI units are km-based; ensure target ECI units are same
56     assert(eci.UnitsAreKm());
57    
58     cVector vecRgRate(eci.getVel().m_x - eciSite.getVel().m_x,
59     eci.getVel().m_y - eciSite.getVel().m_y,
60     eci.getVel().m_z - eciSite.getVel().m_z);
61    
62     double x = eci.getPos().m_x - eciSite.getPos().m_x;
63     double y = eci.getPos().m_y - eciSite.getPos().m_y;
64     double z = eci.getPos().m_z - eciSite.getPos().m_z;
65     double w = sqrt(sqr(x) + sqr(y) + sqr(z));
66    
67     cVector vecRange(x, y, z, w);
68    
69     // The site's Local Mean Sidereal Time at the time of interest.
70     double theta = date.toLMST(getLon());
71    
72     double sin_lat = sin(getLat());
73     double cos_lat = cos(getLat());
74     double sin_theta = sin(theta);
75     double cos_theta = cos(theta);
76    
77     double top_s = sin_lat * cos_theta * vecRange.m_x +
78     sin_lat * sin_theta * vecRange.m_y -
79     cos_lat * vecRange.m_z;
80     double top_e = -sin_theta * vecRange.m_x +
81     cos_theta * vecRange.m_y;
82     double top_z = cos_lat * cos_theta * vecRange.m_x +
83     cos_lat * sin_theta * vecRange.m_y +
84     sin_lat * vecRange.m_z;
85     double az = atan(-top_e / top_s);
86    
87     if (top_s > 0.0)
88     az += PI;
89    
90     if (az < 0.0)
91     az += 2.0*PI;
92    
93     double el = asin(top_z / vecRange.m_w);
94     double rate = (vecRange.m_x * vecRgRate.m_x +
95     vecRange.m_y * vecRgRate.m_y +
96     vecRange.m_z * vecRgRate.m_z) / vecRange.m_w;
97    
98     cCoordTopo topo(az, // azimuth, radians
99     el, // elevation, radians
100     vecRange.m_w, // range, km
101     rate); // rate, km / sec
102    
103     #ifdef WANT_ATMOSPHERIC_CORRECTION
104     // Elevation correction for atmospheric refraction.
105     // Reference: Astronomical Algorithms by Jean Meeus, pp. 101-104
106     // Note: Correction is meaningless when apparent elevation is below horizon
107     topo.m_El += deg2rad((1.02 /
108     tan(deg2rad(rad2deg(el) + 10.3 /
109     (rad2deg(el) + 5.11)))) / 60.0);
110     if (topo.m_El < 0.0)
111     topo.m_El = el; // Reset to true elevation
112    
113     if (topo.m_El > (PI / 2))
114     topo.m_El = (PI / 2);
115     #endif
116    
117     return topo;
118     }
119    
120     //////////////////////////////////////////////////////////////////////////////
121     // toString()
122     //
123     string cSite::toString() const
124     {
125     const int TEMP_SIZE = 128;
126     char sz[TEMP_SIZE];
127    
128     bool LatNorth = true;
129     bool LonEast = true;
130    
131     if (m_geo.m_Lat < 0.0)
132     {
133     LatNorth = false;
134     }
135    
136     if (m_geo.m_Lon < 0.0)
137     {
138     LonEast = false;
139     }
140    
141     snprintf(sz, TEMP_SIZE,
142     "%06.3f%c, ",
143     fabs(rad2deg(m_geo.m_Lat)),
144     (LatNorth ? 'N' : 'S'));
145    
146     string strLoc = sz;
147    
148     snprintf(sz, TEMP_SIZE,
149     "%07.3f%c, ",
150     fabs(rad2deg(m_geo.m_Lon)),
151     (LonEast ? 'E' : 'W'));
152     strLoc += sz;
153    
154     snprintf(sz, TEMP_SIZE,
155     "%.1fm\n",
156     (m_geo.m_Alt * 1000.0));
157     strLoc += sz;
158    
159     return strLoc;
160     }
161    

  ViewVC Help
Powered by ViewVC 1.1.23