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