1   /* Copyright 2022-2025 Thales Alenia Space
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.weather;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.util.FastMath;
21  import org.orekit.bodies.FieldGeodeticPoint;
22  import org.orekit.bodies.GeodeticPoint;
23  import org.orekit.models.earth.weather.water.WaterVaporPressureProvider;
24  import org.orekit.time.AbsoluteDate;
25  import org.orekit.time.FieldAbsoluteDate;
26  
27  /** Converter for weather parameters that change with height.
28   * <p>
29   * Height variations correspond to equations 5.98, 5.99 and 5.100 from
30   * Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007
31   * </p>
32   * @author Luc Maisonobe
33   * @since 12.1
34   */
35  public class HeightDependentPressureTemperatureHumidityConverter {
36  
37      /** Water pressure provider for water vapor pressure. */
38      private final WaterVaporPressureProvider provider;
39  
40      /** Simple constructor.
41       * <p>
42       * Points outside of altitude range will be silently clipped back to range.
43       * </p>
44       * @param provider provider for water vapor pressure
45       */
46      public HeightDependentPressureTemperatureHumidityConverter(final WaterVaporPressureProvider provider) {
47          this.provider = provider;
48      }
49  
50      /** Convert weather parameters.
51       * @param pth0 weather at reference altitude
52       * @param h altitude at which weather is requested
53       * @return converted weather
54       */
55      public PressureTemperatureHumidity convert(final PressureTemperatureHumidity pth0,
56                                                 final double h) {
57  
58          // retrieve parameters at reference altitude
59          final double rh0 = provider.relativeHumidity(pth0.getPressure(), pth0.getTemperature(), pth0.getWaterVaporPressure());
60  
61          // compute changes due to altitude change
62          final double dh = h - pth0.getAltitude();
63          final double p  = pth0.getPressure() * FastMath.pow(1.0 - 2.26e-5 * dh, 5.225);
64          final double t  = pth0.getTemperature() - 6.5e-3 * dh;
65          final double rh = rh0 * FastMath.exp(-6.396e-4 * dh);
66  
67          return new PressureTemperatureHumidity(h, p, t, provider.waterVaporPressure(p, t, rh),
68                                                 pth0.getTm(), pth0.getLambda());
69  
70      }
71  
72      /** Convert weather parameters.
73       * @param <T> type of the elements
74       * @param pth0 weather at reference altitude
75       * @param h altitude at which weather is requested
76       * @return converted weather
77       */
78      public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T> convert(final FieldPressureTemperatureHumidity<T> pth0,
79                                                                                             final T h) {
80          // retrieve parameters at reference altitude
81          final T rh0 = provider.relativeHumidity(pth0.getPressure(), pth0.getTemperature(), pth0.getWaterVaporPressure());
82  
83          // compute changes due to altitude change
84          final T dh = h.subtract(pth0.getAltitude());
85          final T t  = pth0.getTemperature().subtract(dh.multiply(6.5e-3));
86          final T p  = pth0.getPressure().multiply(dh.multiply(2.26e-5).negate().add(1.0).pow(5.225));
87          final T rh = rh0.multiply(FastMath.exp(dh.multiply(-6.396e-4)));
88          return new FieldPressureTemperatureHumidity<>(h, p, t, provider.waterVaporPressure(p, t, rh),
89                                                        pth0.getTm(), pth0.getLambda());
90      }
91  
92      /** Generate a provider applying altitude dependency to fixed weather parameters.
93       * @param basePTH base weather parameters
94       * @return a provider that applies altitude dependency
95       * @since 13.0
96       */
97      public PressureTemperatureHumidityProvider getProvider(final PressureTemperatureHumidity basePTH) {
98          return new PressureTemperatureHumidityProvider() {
99  
100             /** {@inheritDoc} */
101             @Override
102             public PressureTemperatureHumidity getWeatherParameters(final GeodeticPoint location,
103                                                                     final AbsoluteDate date) {
104                 return convert(basePTH, location.getAltitude());
105             }
106 
107             /** {@inheritDoc} */
108             @Override
109             public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T>
110                 getWeatherParameters(final FieldGeodeticPoint<T> location, final FieldAbsoluteDate<T> date) {
111                 return convert(new FieldPressureTemperatureHumidity<>(date.getField(), basePTH),
112                                location.getAltitude());
113             }
114         };
115     }
116 
117 }