1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.troposphere;
18
19 import java.util.Collections;
20 import java.util.List;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.FieldSinCos;
25 import org.hipparchus.util.MathUtils;
26 import org.hipparchus.util.SinCos;
27 import org.orekit.bodies.FieldGeodeticPoint;
28 import org.orekit.bodies.GeodeticPoint;
29 import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
30 import org.orekit.models.earth.weather.PressureTemperatureHumidity;
31 import org.orekit.time.AbsoluteDate;
32 import org.orekit.time.FieldAbsoluteDate;
33 import org.orekit.utils.FieldTrackingCoordinates;
34 import org.orekit.utils.ParameterDriver;
35 import org.orekit.utils.TrackingCoordinates;
36
37
38
39
40
41
42
43
44
45
46
47 public class ModifiedHopfieldModel implements TroposphericModel {
48
49
50 private static final double HD0 = 40136.0;
51
52
53 private static final double HD1 = 148.72;
54
55
56 private static final double T0 = 273.16;
57
58
59 private static final double HW0 = 11000.0;
60
61
62 private static final double ND = 77.64e-6;
63
64
65 private static final double NW1 = -12.96e-6;
66
67
68 private static final double NW2 = 0.371800;
69
70
71 private static final double RE = 6378137.0;
72
73
74
75 public ModifiedHopfieldModel() {
76
77 }
78
79
80 @Override
81 public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
82 final GeodeticPoint point,
83 final PressureTemperatureHumidity weather,
84 final double[] parameters, final AbsoluteDate date) {
85
86
87 final double zenithAngle = MathUtils.SEMI_PI - trackingCoordinates.getElevation();
88
89
90 final double hd = HD0 + HD1 * (weather.getTemperature() - T0);
91 final double nd = ND * TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()) /
92 weather.getTemperature();
93
94
95 final double hw = HW0;
96 final double nw = (NW1 + NW2 / weather.getTemperature()) / weather.getTemperature();
97
98 return new TroposphericDelay(delay(0.0, hd, nd),
99 delay(0.0, hw, nw),
100 delay(zenithAngle, hd, nd),
101 delay(zenithAngle, hw, nw));
102
103 }
104
105
106
107
108
109
110
111
112
113
114
115
116 @Override
117 public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
118 final FieldGeodeticPoint<T> point,
119 final FieldPressureTemperatureHumidity<T> weather,
120 final T[] parameters, final FieldAbsoluteDate<T> date) {
121
122
123 final T zenithAngle = trackingCoordinates.getElevation().negate().add(MathUtils.SEMI_PI);
124
125
126 final T hd = weather.getTemperature().subtract(T0).multiply(HD1).add(HD0);
127 final T nd = TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()).
128 multiply(ND).
129 divide(weather.getTemperature());
130
131
132 final T hw = date.getField().getZero().newInstance(HW0);
133 final T nw = weather.getTemperature().reciprocal().multiply(NW2).add(NW1).divide(weather.getTemperature());
134
135 return new FieldTroposphericDelay<>(delay(date.getField().getZero(), hd, nd),
136 delay(date.getField().getZero(), hw, nw),
137 delay(zenithAngle, hd, nd),
138 delay(zenithAngle, hw, nw));
139
140 }
141
142
143 @Override
144 public List<ParameterDriver> getParametersDrivers() {
145 return Collections.emptyList();
146 }
147
148
149
150
151
152
153
154 private double delay(final double zenithAngle, final double hi, final double ni) {
155
156
157 final SinCos scZ = FastMath.sinCos(zenithAngle);
158 final double rePhi = RE + hi;
159 final double reS = RE * scZ.sin();
160 final double reC = RE * scZ.cos();
161 final double ri = FastMath.sqrt(rePhi * rePhi - reS * reS) - reC;
162
163 final double ai = -scZ.cos() / hi;
164 final double bi = -scZ.sin() * scZ.sin() / (2 * hi * RE);
165 final double ai2 = ai * ai;
166 final double bi2 = bi * bi;
167
168 final double f1i = 1;
169 final double f2i = 4 * ai;
170 final double f3i = 6 * ai2 + 4 * bi;
171 final double f4i = 4 * ai * (ai2 + 3 * bi);
172 final double f5i = ai2 * ai2 + 12 * ai2 * bi + 6 * bi2;
173 final double f6i = 4 * ai * bi * (ai2 + 3 * bi);
174 final double f7i = bi2 * (6 * ai2 + 4 * bi);
175 final double f8i = 4 * ai * bi * bi2;
176 final double f9i = bi2 * bi2;
177
178 return ni * (ri * (f1i +
179 ri * (f2i / 2 +
180 ri * (f3i / 3 +
181 ri * (f4i / 4 +
182 ri * (f5i / 5 +
183 ri * (f6i / 6 +
184 ri * (f7i / 7 +
185 ri * (f8i / 8 +
186 ri * f9i / 9)))))))));
187
188 }
189
190
191
192
193
194
195
196
197 private <T extends CalculusFieldElement<T>> T delay(final T zenithAngle, final T hi, final T ni) {
198
199
200 final FieldSinCos<T> scZ = FastMath.sinCos(zenithAngle);
201 final T rePhi = hi.add(RE);
202 final T reS = scZ.sin().multiply(RE);
203 final T reC = scZ.cos().multiply(RE);
204 final T ri = FastMath.sqrt(rePhi.multiply(rePhi).subtract(reS.multiply(reS))).subtract(reC);
205
206 final T ai = scZ.cos().negate().divide(hi);
207 final T bi = scZ.sin().multiply(scZ.sin()).negate().divide(hi.add(hi).multiply(RE));
208 final T ai2 = ai.multiply(ai);
209 final T bi2 = bi.multiply(bi);
210
211 final T f1i = ai.getField().getOne();
212 final T f2i = ai.multiply(4);
213 final T f3i = ai2.multiply(6).add(bi.multiply(4));
214 final T f4i = ai.multiply(4).multiply(ai2.add(bi.multiply(3)));
215 final T f5i = ai2.multiply(ai2).add(ai2.multiply(12).multiply(bi)).add(bi2.multiply(6));
216 final T f6i = ai.multiply(4).multiply(bi).multiply(ai2.add(bi.multiply(3)));
217 final T f7i = bi2.multiply(ai2.multiply(6).add(bi.multiply(4)));
218 final T f8i = ai.multiply(4).multiply(bi).multiply(bi2);
219 final T f9i = bi2.multiply(bi2);
220
221 return ni.
222 multiply(ri.
223 multiply(f1i.
224 add(ri.
225 multiply(f2i.divide(2).
226 add(ri.
227 multiply(f3i.divide(3).
228 add(ri.
229 multiply(f4i.divide(4).
230 add(ri.
231 multiply(f5i.divide(5).
232 add(ri.
233 multiply(f6i.divide(6).
234 add(ri.
235 multiply(f7i.divide(7).
236 add(ri.
237 multiply(f8i.divide(8).
238 add(ri.multiply(f9i.divide(9)))))))))))))))))));
239
240 }
241
242 }