|
1
|
|
|
<?php |
|
2
|
|
|
|
|
3
|
|
|
/** |
|
4
|
|
|
* Predict_SGPObs |
|
5
|
|
|
* |
|
6
|
|
|
* Ported to PHP by Bill Shupp. Original comments below |
|
7
|
|
|
*/ |
|
8
|
|
|
|
|
9
|
|
|
require_once 'Math.php'; |
|
10
|
|
|
require_once 'Time.php'; |
|
11
|
|
|
//require_once 'Predict.php'; |
|
12
|
|
|
require_once 'Vector.php'; |
|
13
|
|
|
require_once 'ObsSet.php'; |
|
14
|
|
|
|
|
15
|
|
|
/* |
|
16
|
|
|
* Unit SGP_Obs |
|
17
|
|
|
* Author: Dr TS Kelso |
|
18
|
|
|
* Original Version: 1992 Jun 02 |
|
19
|
|
|
* Current Revision: 1992 Sep 28 |
|
20
|
|
|
* Version: 1.40 |
|
21
|
|
|
* Copyright: 1992, All Rights Reserved |
|
22
|
|
|
* |
|
23
|
|
|
* Ported to C by: Neoklis Kyriazis April 9 2001 |
|
24
|
|
|
* Ported to PHP by Bill Shupp August, 2011 |
|
25
|
|
|
*/ |
|
26
|
|
|
class Predict_SGPObs |
|
27
|
|
|
{ |
|
28
|
|
|
/* Procedure Calculate_User_PosVel passes the user's geodetic position */ |
|
29
|
|
|
/* and the time of interest and returns the ECI position and velocity */ |
|
30
|
|
|
/* of the observer. The velocity calculation assumes the geodetic */ |
|
31
|
|
|
/* position is stationary relative to the earth's surface. */ |
|
32
|
|
|
public static function Calculate_User_PosVel( |
|
33
|
|
|
$_time, Predict_Geodetic $geodetic, Predict_Vector $obs_pos, Predict_Vector $obs_vel |
|
34
|
|
|
) |
|
35
|
|
|
{ |
|
36
|
|
|
/* Reference: The 1992 Astronomical Almanac, page K11. */ |
|
37
|
|
|
|
|
38
|
|
|
$sinGeodeticLat = sin($geodetic->lat); /* Only run sin($geodetic->lat) once */ |
|
39
|
|
|
|
|
40
|
|
|
$geodetic->theta = Predict_Math::FMod2p(Predict_Time::ThetaG_JD($_time) + $geodetic->lon);/*LMST*/ |
|
41
|
|
|
$c = 1 / sqrt(1 + Predict::__f * (Predict::__f - 2) * $sinGeodeticLat * $sinGeodeticLat); |
|
42
|
|
|
$sq = (1 - Predict::__f) * (1 - Predict::__f) * $c; |
|
43
|
|
|
$achcp = (Predict::xkmper * $c + $geodetic->alt) * cos($geodetic->lat); |
|
44
|
|
|
$obs_pos->x = $achcp * cos($geodetic->theta); /*kilometers*/ |
|
|
|
|
|
|
45
|
|
|
$obs_pos->y = $achcp * sin($geodetic->theta); |
|
|
|
|
|
|
46
|
|
|
$obs_pos->z = (Predict::xkmper * $sq + $geodetic->alt) * $sinGeodeticLat; |
|
|
|
|
|
|
47
|
|
|
$obs_vel->x = -Predict::mfactor * $obs_pos->y; /*kilometers/second*/ |
|
|
|
|
|
|
48
|
|
|
$obs_vel->y = Predict::mfactor * $obs_pos->x; |
|
|
|
|
|
|
49
|
|
|
$obs_vel->z = 0; |
|
50
|
|
|
$obs_pos->w = sqrt($obs_pos->x * $obs_pos->x + $obs_pos->y * $obs_pos->y + $obs_pos->z * $obs_pos->z); |
|
|
|
|
|
|
51
|
|
|
$obs_vel->w = sqrt($obs_vel->x * $obs_vel->x + $obs_vel->y * $obs_vel->y + $obs_vel->z * $obs_vel->z); |
|
|
|
|
|
|
52
|
|
|
} |
|
53
|
|
|
|
|
54
|
|
|
/* Procedure Calculate_LatLonAlt will calculate the geodetic */ |
|
55
|
|
|
/* position of an object given its ECI position pos and time. */ |
|
56
|
|
|
/* It is intended to be used to determine the ground track of */ |
|
57
|
|
|
/* a satellite. The calculations assume the earth to be an */ |
|
58
|
|
|
/* oblate spheroid as defined in WGS '72. */ |
|
59
|
|
|
public static function Calculate_LatLonAlt($_time, Predict_Vector $pos, Predict_Geodetic $geodetic) |
|
60
|
|
|
{ |
|
61
|
|
|
/* Reference: The 1992 Astronomical Almanac, page K12. */ |
|
62
|
|
|
|
|
63
|
|
|
/* double r,e2,phi,c; */ |
|
64
|
|
|
|
|
65
|
|
|
$geodetic->theta = Predict_Math::AcTan($pos->y, $pos->x); /*radians*/ |
|
66
|
|
|
$geodetic->lon = Predict_Math::FMod2p($geodetic->theta - Predict_Time::ThetaG_JD($_time)); /*radians*/ |
|
67
|
|
|
$r = sqrt(($pos->x * $pos->x) + ($pos->y * $pos->y)); |
|
68
|
|
|
$e2 = Predict::__f * (2 - Predict::__f); |
|
69
|
|
|
$geodetic->lat = Predict_Math::AcTan($pos->z, $r); /*radians*/ |
|
70
|
|
|
|
|
71
|
|
|
do { |
|
72
|
|
|
$phi = $geodetic->lat; |
|
73
|
|
|
$sinPhi = sin($phi); |
|
74
|
|
|
$c = 1 / sqrt(1 - $e2 * ($sinPhi * $sinPhi)); |
|
75
|
|
|
$geodetic->lat = Predict_Math::AcTan($pos->z + Predict::xkmper * $c * $e2 * $sinPhi, $r); |
|
76
|
|
|
} while (abs($geodetic->lat - $phi) >= 1E-10); |
|
77
|
|
|
|
|
78
|
|
|
$geodetic->alt = $r / cos($geodetic->lat) - Predict::xkmper * $c;/*kilometers*/ |
|
79
|
|
|
|
|
80
|
|
|
if ($geodetic->lat > Predict::pio2) { |
|
81
|
|
|
$geodetic->lat -= Predict::twopi; |
|
82
|
|
|
} |
|
83
|
|
|
} |
|
84
|
|
|
|
|
85
|
|
|
/* The procedures Calculate_Obs and Calculate_RADec calculate */ |
|
86
|
|
|
/* the *topocentric* coordinates of the object with ECI position, */ |
|
87
|
|
|
/* {pos}, and velocity, {vel}, from location {geodetic} at {time}. */ |
|
88
|
|
|
/* The {obs_set} returned for Calculate_Obs consists of azimuth, */ |
|
89
|
|
|
/* elevation, range, and range rate (in that order) with units of */ |
|
90
|
|
|
/* radians, radians, kilometers, and kilometers/second, respectively. */ |
|
91
|
|
|
/* The WGS '72 geoid is used and the effect of atmospheric refraction */ |
|
92
|
|
|
/* (under standard temperature and pressure) is incorporated into the */ |
|
93
|
|
|
/* elevation calculation; the effect of atmospheric refraction on */ |
|
94
|
|
|
/* range and range rate has not yet been quantified. */ |
|
95
|
|
|
|
|
96
|
|
|
/* The {obs_set} for Calculate_RADec consists of right ascension and */ |
|
97
|
|
|
/* declination (in that order) in radians. Again, calculations are */ |
|
98
|
|
|
/* based on *topocentric* position using the WGS '72 geoid and */ |
|
99
|
|
|
/* incorporating atmospheric refraction. */ |
|
100
|
|
|
public static function Calculate_Obs($_time, Predict_Vector $pos, Predict_Vector $vel, Predict_Geodetic $geodetic, Predict_ObsSet $obs_set) |
|
101
|
|
|
{ |
|
102
|
|
|
$obs_pos = new Predict_Vector(); |
|
103
|
|
|
$obs_vel = new Predict_Vector(); |
|
104
|
|
|
$range = new Predict_Vector(); |
|
105
|
|
|
$rgvel = new Predict_Vector(); |
|
106
|
|
|
|
|
107
|
|
|
self::Calculate_User_PosVel($_time, $geodetic, $obs_pos, $obs_vel); |
|
108
|
|
|
|
|
109
|
|
|
$range->x = $pos->x - $obs_pos->x; |
|
110
|
|
|
$range->y = $pos->y - $obs_pos->y; |
|
111
|
|
|
$range->z = $pos->z - $obs_pos->z; |
|
112
|
|
|
|
|
113
|
|
|
$rgvel->x = $vel->x - $obs_vel->x; |
|
114
|
|
|
$rgvel->y = $vel->y - $obs_vel->y; |
|
115
|
|
|
$rgvel->z = $vel->z - $obs_vel->z; |
|
116
|
|
|
|
|
117
|
|
|
$range->w = sqrt($range->x * $range->x + $range->y * $range->y + $range->z * $range->z); |
|
|
|
|
|
|
118
|
|
|
|
|
119
|
|
|
$sin_lat = sin($geodetic->lat); |
|
120
|
|
|
$cos_lat = cos($geodetic->lat); |
|
121
|
|
|
$sin_theta = sin($geodetic->theta); |
|
122
|
|
|
$cos_theta = cos($geodetic->theta); |
|
123
|
|
|
$top_s = $sin_lat * $cos_theta * $range->x |
|
124
|
|
|
+ $sin_lat * $sin_theta * $range->y |
|
125
|
|
|
- $cos_lat * $range->z; |
|
126
|
|
|
$top_e = -$sin_theta * $range->x |
|
127
|
|
|
+ $cos_theta * $range->y; |
|
128
|
|
|
$top_z = $cos_lat * $cos_theta * $range->x |
|
129
|
|
|
+ $cos_lat * $sin_theta * $range->y |
|
130
|
|
|
+ $sin_lat * $range->z; |
|
131
|
|
|
$azim = atan(-$top_e / $top_s); /*Azimuth*/ |
|
132
|
|
|
if ($top_s > 0) { |
|
133
|
|
|
$azim = $azim + Predict::pi; |
|
134
|
|
|
} |
|
135
|
|
|
if ($azim < 0 ) { |
|
136
|
|
|
$azim = $azim + Predict::twopi; |
|
137
|
|
|
} |
|
138
|
|
|
$el = Predict_Math::ArcSin($top_z / $range->w); |
|
139
|
|
|
$obs_set->az = $azim; /* Azimuth (radians) */ |
|
140
|
|
|
$obs_set->el = $el; /* Elevation (radians)*/ |
|
141
|
|
|
$obs_set->range = $range->w; /* Range (kilometers) */ |
|
142
|
|
|
|
|
143
|
|
|
/* Range Rate (kilometers/second)*/ |
|
144
|
|
|
$obs_set->range_rate = Predict_Math::Dot($range, $rgvel) / $range->w; |
|
145
|
|
|
|
|
146
|
|
|
/* Corrections for atmospheric refraction */ |
|
147
|
|
|
/* Reference: Astronomical Algorithms by Jean Meeus, pp. 101-104 */ |
|
148
|
|
|
/* Correction is meaningless when apparent elevation is below horizon */ |
|
149
|
|
|
// obs_set->el = obs_set->el + Radians((1.02/tan(Radians(Degrees(el)+ |
|
150
|
|
|
// 10.3/(Degrees(el)+5.11))))/60); |
|
151
|
|
|
if ($obs_set->el < 0) { |
|
152
|
|
|
$obs_set->el = $el; /*Reset to true elevation*/ |
|
153
|
|
|
} |
|
154
|
|
|
} |
|
155
|
|
|
} |
|
156
|
|
|
|
Ysurac /
FlightAirMap
Loading history...
This check looks for assignments to scalar types that may be of the wrong type.
To ensure the code behaves as expected, it may be a good idea to add an explicit type cast.