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.orekit.bodies.FieldGeodeticPoint;
25 import org.orekit.bodies.GeodeticPoint;
26 import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
27 import org.orekit.models.earth.weather.PressureTemperatureHumidity;
28 import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
29 import org.orekit.time.AbsoluteDate;
30 import org.orekit.time.FieldAbsoluteDate;
31 import org.orekit.utils.FieldTrackingCoordinates;
32 import org.orekit.utils.ParameterDriver;
33 import org.orekit.utils.TrackingCoordinates;
34 import org.orekit.utils.units.Unit;
35 import org.orekit.utils.units.UnitsConverter;
36
37
38
39
40
41
42
43
44
45
46 public class MariniMurray implements TroposphericModel {
47
48
49 private final double fLambda;
50
51
52
53
54 private final PressureTemperatureHumidityProvider pthProvider;
55
56
57
58
59
60
61
62
63
64 public MariniMurray(final double lambda, final Unit lambdaUnits, final PressureTemperatureHumidityProvider pthProvider) {
65
66 this.pthProvider = pthProvider;
67
68
69 final double lambdaMicrometer = new UnitsConverter(lambdaUnits, TroposphericModelUtils.MICRO_M).convert(lambda);
70 final double l2 = lambdaMicrometer * lambdaMicrometer;
71 this.fLambda = 0.9650 + (0.0164 + 0.000228 / l2) / l2;
72
73 }
74
75
76 @Override
77 public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
78 final double[] parameters, final AbsoluteDate date) {
79
80 final PressureTemperatureHumidity weather = pthProvider.getWeatherParameters(point, date);
81 final double p = weather.getPressure();
82 final double t = weather.getTemperature();
83 final double e = weather.getWaterVaporPressure();
84
85
86 final double Ah = 0.00002357 * p;
87 final double Aw = 0.00000141 * e;
88 final double K = 1.163 - 0.00968 * FastMath.cos(2 * point.getLatitude()) - 0.00104 * t + 0.0000001435 * p;
89 final double B = 1.084e-10 * p * t * K + 4.734e-12 * p * (p / t) * (2 * K) / (3 * K - 1);
90 final double flambda = getLaserFrequencyParameter();
91
92 final double fsite = getSiteFunctionValue(point);
93
94 final double sinE = FastMath.sin(trackingCoordinates.getElevation());
95 final double totalZenith = (flambda / fsite) * (Ah + Aw + B) / (1.0 + B / ((Ah + Aw + B) * (1.0 + 0.01)));
96 final double totalElev = (flambda / fsite) * (Ah + Aw + B) / (sinE + B / ((Ah + Aw + B) * (sinE + 0.01)));
97 final double hydrostaticZenith = (flambda / fsite) * (Ah + B) / (1.0 + B / ((Ah + B) * (1.0 + 0.01)));
98 final double hydrostaticElev = (flambda / fsite) * (Ah + B) / (sinE + B / ((Ah + B) * (sinE + 0.01)));
99 return new TroposphericDelay(hydrostaticZenith, totalZenith - hydrostaticZenith,
100 hydrostaticElev, totalElev - hydrostaticElev);
101 }
102
103
104 @Override
105 public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
106 final FieldGeodeticPoint<T> point,
107 final T[] parameters, final FieldAbsoluteDate<T> date) {
108
109 final FieldPressureTemperatureHumidity<T> weather = pthProvider.getWeatherParameters(point, date);
110 final T p = weather.getPressure();
111 final T t = weather.getTemperature();
112 final T e = weather.getWaterVaporPressure();
113
114
115 final T Ah = p.multiply(0.00002357);
116 final T Aw = e.multiply(0.00000141);
117 final T K = FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00968).negate().
118 add(1.163).
119 subtract(t.multiply(0.00104)).
120 add(p.multiply(0.0000001435));
121 final T B = K.multiply(t.multiply(p).multiply(1.084e-10 )).
122 add(K.multiply(2.).multiply(p.multiply(p).divide(t).multiply(4.734e-12)).divide(K.multiply(3.).subtract(1.)));
123 final double flambda = getLaserFrequencyParameter();
124
125 final T fsite = getSiteFunctionValue(point);
126
127 final T sinE = FastMath.sin(trackingCoordinates.getElevation());
128 final T one = date.getField().getOne();
129 final T totalZenith = fsite.divide(flambda).reciprocal().
130 multiply(B.add(Ah).add(Aw)).
131 divide(one.add(one.add(0.01).multiply(B.add(Ah).add(Aw)).divide(B).reciprocal()));
132 final T totalElev = fsite.divide(flambda).reciprocal().
133 multiply(B.add(Ah).add(Aw)).
134 divide(sinE.add(sinE.add(0.01).multiply(B.add(Ah).add(Aw)).divide(B).reciprocal()));
135 final T hydrostaticZenith = fsite.divide(flambda).reciprocal().
136 multiply(B.add(Ah)).
137 divide(one.add(one.add(0.01).multiply(B.add(Ah)).divide(B).reciprocal()));
138 final T hydrostaticElev = fsite.divide(flambda).reciprocal().
139 multiply(B.add(Ah)).
140 divide(sinE.add(sinE.add(0.01).multiply(B.add(Ah)).divide(B).reciprocal()));
141 return new FieldTroposphericDelay<>(hydrostaticZenith, totalZenith.subtract(hydrostaticZenith),
142 hydrostaticElev, totalElev.subtract(hydrostaticElev));
143 }
144
145
146 @Override
147 public List<ParameterDriver> getParametersDrivers() {
148 return Collections.emptyList();
149 }
150
151
152
153
154
155
156
157 private double getLaserFrequencyParameter() {
158 return fLambda;
159 }
160
161
162
163
164
165
166 private double getSiteFunctionValue(final GeodeticPoint point) {
167 return 1. - 0.0026 * FastMath.cos(2 * point.getLatitude()) - 0.00031 * 0.001 * point.getAltitude();
168 }
169
170
171
172
173
174
175
176 private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
177 return FastMath.cos(point.getLatitude().multiply(2)).multiply(0.0026).add(point.getAltitude().multiply(0.001).multiply(0.00031)).negate().add(1.);
178 }
179
180 }