1 /* Copyright 2011-2012 Space Applications Services
2 * Licensed to CS Communication & Systèmes (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
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.time.AbsoluteDate;
27 import org.orekit.time.FieldAbsoluteDate;
28 import org.orekit.utils.ParameterDriver;
29
30 /** The Marini-Murray tropospheric delay model for laser ranging.
31 *
32 * @see "Marini, J.W., and C.W. Murray, correction of Laser Range Tracking Data for
33 * Atmospheric Refraction at Elevations Above 10 degrees, X-591-73-351, NASA GSFC, 1973"
34 *
35 * @author Joris Olympio
36 */
37 public class MariniMurrayModel implements DiscreteTroposphericModel {
38
39 /** The temperature at the station, K. */
40 private double T0;
41
42 /** The atmospheric pressure, mbar. */
43 private double P0;
44
45 /** water vapor pressure at the laser site, mbar. */
46 private double e0;
47
48 /** Laser wavelength, micrometers. */
49 private double lambda;
50
51 /** Create a new Marini-Murray model for the troposphere using the given
52 * environmental conditions.
53 * @param t0 the temperature at the station, K
54 * @param p0 the atmospheric pressure at the station, mbar
55 * @param rh the humidity at the station, percent (50% -> 0.5)
56 * @param lambda laser wavelength (c/f), nm
57 */
58 public MariniMurrayModel(final double t0, final double p0, final double rh, final double lambda) {
59 this.T0 = t0;
60 this.P0 = p0;
61 this.e0 = getWaterVapor(rh);
62 this.lambda = lambda * 1e-3;
63 }
64
65 /** Create a new Marini-Murray model using a standard atmosphere model.
66 *
67 * <ul>
68 * <li>temperature: 20 degree Celsius</li>
69 * <li>pressure: 1013.25 mbar</li>
70 * <li>humidity: 50%</li>
71 * </ul>
72 *
73 * @param lambda laser wavelength (c/f), nm
74 *
75 * @return a Marini-Murray model with standard environmental values
76 */
77 public static MariniMurrayModel getStandardModel(final double lambda) {
78 return new MariniMurrayModel(273.15 + 20, 1013.25, 0.5, lambda);
79 }
80
81 /** {@inheritDoc} */
82 @Override
83 public double pathDelay(final double elevation, final GeodeticPoint point,
84 final double[] parameters, final AbsoluteDate date) {
85 final double A = 0.002357 * P0 + 0.000141 * e0;
86 final double K = 1.163 - 0.00968 * FastMath.cos(2 * point.getLatitude()) - 0.00104 * T0 + 0.00001435 * P0;
87 final double B = (1.084 * 1e-8) * P0 * T0 * K + (4.734 * 1e-8) * P0 * (P0 / T0) * (2 * K) / (3 * K - 1);
88 final double flambda = getLaserFrequencyParameter();
89
90 final double fsite = getSiteFunctionValue(point);
91
92 final double sinE = FastMath.sin(elevation);
93 final double dR = (flambda / fsite) * (A + B) / (sinE + B / ((A + B) * (sinE + 0.01)) );
94 return dR;
95 }
96
97 /** {@inheritDoc} */
98 @Override
99 public <T extends CalculusFieldElement<T>> T pathDelay(final T elevation, final FieldGeodeticPoint<T> point,
100 final T[] parameters, final FieldAbsoluteDate<T> date) {
101 final double A = 0.002357 * P0 + 0.000141 * e0;
102 final T K = FastMath.cos(point.getLatitude().multiply(2.)).multiply(0.00968).negate().add(1.163).subtract(0.00104 * T0).add(0.00001435 * P0);
103 final T B = K.multiply((1.084 * 1e-8) * P0 * T0).add(K.multiply(2.).multiply((4.734 * 1e-8) * P0 * (P0 / T0)).divide(K.multiply(3.).subtract(1.)));
104 final double flambda = getLaserFrequencyParameter();
105
106 final T fsite = getSiteFunctionValue(point);
107
108 final T sinE = FastMath.sin(elevation);
109 final T dR = fsite.divide(flambda).reciprocal().multiply(B.add(A)).divide(sinE.add(sinE.add(0.01).multiply(B.add(A)).divide(B).reciprocal()));
110 return dR;
111 }
112
113 /** {@inheritDoc} */
114 @Override
115 public List<ParameterDriver> getParametersDrivers() {
116 return Collections.emptyList();
117 }
118
119 /** Get the laser frequency parameter f(lambda).
120 * It is one for Ruby laser (lambda = 0.6943 micron)
121 * For infrared lasers, f(lambda) = 0.97966.
122 *
123 * @return the laser frequency parameter f(lambda).
124 */
125 private double getLaserFrequencyParameter() {
126 return 0.9650 + 0.0164 * FastMath.pow(lambda, -2) + 0.000228 * FastMath.pow(lambda, -4);
127 }
128
129 /** Get the laser frequency parameter f(lambda).
130 *
131 * @param point station location
132 * @return the laser frequency parameter f(lambda).
133 */
134 private double getSiteFunctionValue(final GeodeticPoint point) {
135 return 1. - 0.0026 * FastMath.cos(2 * point.getLatitude()) - 0.00031 * 0.001 * point.getAltitude();
136 }
137
138 /** Get the laser frequency parameter f(lambda).
139 *
140 * @param <T> type of the elements
141 * @param point station location
142 * @return the laser frequency parameter f(lambda).
143 */
144 private <T extends CalculusFieldElement<T>> T getSiteFunctionValue(final FieldGeodeticPoint<T> point) {
145 return FastMath.cos(point.getLatitude().multiply(2)).multiply(0.0026).add(point.getAltitude().multiply(0.001).multiply(0.00031)).negate().add(1.);
146 }
147
148 /** Get the water vapor.
149 * The water vapor model is the one of Giacomo and Davis as indicated in IERS TN 32, chap. 9.
150 *
151 * See: Giacomo, P., Equation for the dertermination of the density of moist air, Metrologia, V. 18, 1982
152 *
153 * @param rh relative humidity, in percent (50% -> 0.5).
154 * @return the water vapor, in mbar (1 mbar = 100 Pa).
155 */
156 private double getWaterVapor(final double rh) {
157
158 // saturation water vapor, equation (3) of reference paper, in mbar
159 // with amended 1991 values (see reference paper)
160 final double es = 0.01 * FastMath.exp((1.2378847 * 1e-5) * T0 * T0 -
161 (1.9121316 * 1e-2) * T0 +
162 33.93711047 -
163 (6.3431645 * 1e3) * 1. / T0);
164
165 // enhancement factor, equation (4) of reference paper
166 final double fw = 1.00062 + (3.14 * 1e-6) * P0 + (5.6 * 1e-7) * FastMath.pow(T0 - 273.15, 2);
167
168 final double e = rh * fw * es;
169 return e;
170 }
171
172 }