/* Copyright 2002-2024 CS GROUP
    * Licensed to CS GROUP (CS) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * CS licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *   http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.orekit.models.earth.weather;
  
  import java.io.BufferedInputStream;
  import java.io.IOException;
  import java.io.InputStream;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.util.FastMath;
  import org.orekit.bodies.FieldGeodeticPoint;
  import org.orekit.bodies.GeodeticPoint;
  import org.orekit.data.DataSource;
  import org.orekit.models.earth.troposphere.AzimuthalGradientCoefficients;
  import org.orekit.models.earth.troposphere.AzimuthalGradientProvider;
  import org.orekit.models.earth.troposphere.FieldAzimuthalGradientCoefficients;
  import org.orekit.models.earth.troposphere.FieldViennaACoefficients;
  import org.orekit.models.earth.troposphere.TroposphericModelUtils;
  import org.orekit.models.earth.troposphere.ViennaACoefficients;
  import org.orekit.models.earth.troposphere.ViennaAProvider;
  import org.orekit.time.AbsoluteDate;
  import org.orekit.time.FieldAbsoluteDate;
  import org.orekit.time.TimeScale;
  import org.orekit.utils.Constants;
  
  /** Base class for Global Pressure and Temperature 2, 2w and 3 models.
   * These models are empirical models that provide the temperature, the pressure and the water vapor pressure
   * of a site depending its latitude and  longitude. These models also {@link ViennaACoefficients provide}
   * the a<sub>h</sub> and a<sub>w</sub> coefficients for Vienna models.
   * <p>
   * The requisite coefficients for the computation of the weather parameters are provided by the
   * Department of Geodesy and Geoinformation of the Vienna University. They are based on an
   * external grid file like "gpt2_1.grd" (1° x 1°), "gpt2_5.grd" (5° x 5°), "gpt2_1w.grd" (1° x 1°),
   * "gpt2_5w.grd" (5° x 5°), "gpt3_1.grd" (1° x 1°), or "gpt3_5.grd" (5° x 5°) available at:
   * <a href="https://vmf.geo.tuwien.ac.at/codes/"> link</a>
   * </p>
   * <p>
   * A bilinear interpolation is performed in order to obtained the correct values of the weather parameters.
   * </p>
   * <p>
   * The format is always the same, with and example shown below for the pressure and the temperature.
   * The "GPT2w" model (w stands for wet) also provide humidity parameters and the "GPT3" model also
   * provides horizontal gradient, so the number of columns vary depending on the model.
   * <p>
   * Example:
   * </p>
   * <pre>
   * %  lat    lon   p:a0    A1   B1   A2   B2  T:a0    A1   B1   A2   B2
   *   87.5    2.5 101421    21  409 -217 -122 259.2 -13.2 -6.1  2.6  0.3
   *   87.5    7.5 101416    21  411 -213 -120 259.3 -13.1 -6.1  2.6  0.3
   *   87.5   12.5 101411    22  413 -209 -118 259.3 -13.1 -6.1  2.6  0.3
   *   87.5   17.5 101407    23  415 -205 -116 259.4 -13.0 -6.1  2.6  0.3
   *   ...
   * </pre>
   *
   * @see "K. Lagler, M. Schindelegger, J. Böhm, H. Krasna, T. Nilsson (2013),
   * GPT2: empirical slant delay model for radio space geodetic techniques. Geophys
   * Res Lett 40(6):1069–1073. doi:10.1002/grl.50288"
   *
   * @author Bryan Cazabonne
   * @author Luc Maisonobe
   * @since 12.1
   */
  public class AbstractGlobalPressureTemperature
      implements ViennaAProvider, AzimuthalGradientProvider, PressureTemperatureHumidityProvider {
  
      /** Standard gravity constant [m/s²]. */
      private static final double G = Constants.G0_STANDARD_GRAVITY;
  
      /** Ideal gas constant for dry air [J/kg/K]. */
      private static final double R = 287.0;
  
      /** Loaded grid. */
      private final Grid grid;
  
      /** UTC time scale. */
      private final TimeScale utc;
  
      /**
       * Constructor with source of GPTn auxiliary data given by user.
       *
       * @param source grid data source
       * @param utc UTC time scale.
       * @param expected expected seasonal models types
       * @exception IOException if grid data cannot be read
      */
     protected AbstractGlobalPressureTemperature(final DataSource source, final TimeScale utc,
                                                 final SeasonalModelType... expected)
         throws IOException {
         this.utc = utc;
 
         // load the grid data
         try (InputStream         is     = source.getOpener().openStreamOnce();
              BufferedInputStream bis    = new BufferedInputStream(is)) {
             final GptNParser     parser = new GptNParser(expected);
             parser.loadData(bis, source.getName());
             grid = parser.getGrid();
         }
 
     }
 
     /**
      * Constructor with already loaded grid.
      *
      * @param grid loaded grid
      * @param utc UTC time scale.
      * @deprecated as of 12.1 used only by {@link GlobalPressureTemperature2Model}
      */
     @Deprecated
     protected AbstractGlobalPressureTemperature(final Grid grid, final TimeScale utc) {
         this.grid = grid;
         this.utc  = utc;
     }
 
     /** {@inheritDoc} */
     @Override
     public ViennaACoefficients getA(final GeodeticPoint location, final AbsoluteDate date) {
 
         // set up interpolation parameters
         final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();
 
         // ah and aw coefficients
         return new ViennaACoefficients(interpolator.interpolate(e -> e.getModel(SeasonalModelType.AH).evaluate(dayOfYear)) * 0.001,
                                        interpolator.interpolate(e -> e.getModel(SeasonalModelType.AW).evaluate(dayOfYear)) * 0.001);
 
     }
 
     /** {@inheritDoc} */
     @Override
     public PressureTemperatureHumidity getWeatherParamerers(final GeodeticPoint location, final AbsoluteDate date) {
 
         // set up interpolation parameters
         final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();
 
         // Corrected height (can be negative)
         final double undu            = interpolator.interpolate(e -> e.getUndulation());
         final double correctedheight = location.getAltitude() - undu - interpolator.interpolate(e -> e.getHs());
 
         // Temperature gradient [K/m]
         final double dTdH = interpolator.interpolate(e -> e.getModel(SeasonalModelType.DT).evaluate(dayOfYear)) * 0.001;
 
         // Specific humidity
         final double qv = interpolator.interpolate(e -> e.getModel(SeasonalModelType.QV).evaluate(dayOfYear)) * 0.001;
 
         // For the computation of the temperature and the pressure, we use
         // the standard ICAO atmosphere formulas.
 
         // Temperature [K]
         final double t0 = interpolator.interpolate(e -> e.getModel(SeasonalModelType.TEMPERATURE).evaluate(dayOfYear));
         final double temperature = t0 + dTdH * correctedheight;
 
         // Pressure [hPa]
         final double p0       = interpolator.interpolate(e -> e.getModel(SeasonalModelType.PRESSURE).evaluate(dayOfYear));
         final double exponent = G / (dTdH * R);
         final double pressure = p0 * FastMath.pow(1 - (dTdH / t0) * correctedheight, exponent) * 0.01;
 
         // Water vapor pressure [hPa]
         final double e0 = qv * pressure / (0.622 + 0.378 * qv);
 
         // mean temperature weighted with water vapor pressure
         final double tm = grid.hasModels(SeasonalModelType.TM) ?
                           interpolator.interpolate(e -> e.getModel(SeasonalModelType.TM).evaluate(dayOfYear)) :
                           Double.NaN;
 
         // water vapor decrease factor
         final double lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
                               interpolator.interpolate(e -> e.getModel(SeasonalModelType.LAMBDA).evaluate(dayOfYear)) :
                               Double.NaN;
 
         return new PressureTemperatureHumidity(location.getAltitude(),
                                                TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
                                                temperature,
                                                TroposphericModelUtils.HECTO_PASCAL.toSI(e0),
                                                tm,
                                                lambda);
 
     }
 
     /** {@inheritDoc} */
     @Override
     public AzimuthalGradientCoefficients getGradientCoefficients(final GeodeticPoint location, final AbsoluteDate date) {
 
         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
             // set up interpolation parameters
             final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
             final int              dayOfYear    = date.getComponents(utc).getDate().getDayOfYear();
 
             return new AzimuthalGradientCoefficients(interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_H).evaluate(dayOfYear)),
                                                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_H).evaluate(dayOfYear)),
                                                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_W).evaluate(dayOfYear)),
                                                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_W).evaluate(dayOfYear)));
         } else {
             return null;
         }
 
     }
 
     /** {@inheritDoc} */
     @Override
     public <T extends CalculusFieldElement<T>> FieldViennaACoefficients<T> getA(final FieldGeodeticPoint<T> location,
                                                                                 final FieldAbsoluteDate<T> date) {
 
         // set up interpolation parameters
         final FieldCellInterpolator<T> interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();
 
         // ah and aw coefficients
         return new FieldViennaACoefficients<>(interpolator.interpolate(e -> e.getModel(SeasonalModelType.AH).evaluate(dayOfYear)).multiply(0.001),
                                               interpolator.interpolate(e -> e.getModel(SeasonalModelType.AW).evaluate(dayOfYear)).multiply(0.001));
 
     }
 
     /** {@inheritDoc} */
     @Override
     public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T> getWeatherParamerers(final FieldGeodeticPoint<T> location,
                                                                                                         final FieldAbsoluteDate<T> date) {
 
         // set up interpolation parameters
         final FieldCellInterpolator<T> interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
         final int dayOfYear = date.getComponents(utc).getDate().getDayOfYear();
 
         // Corrected height (can be negative)
         final T undu            = interpolator.interpolate(e -> e.getUndulation());
         final T correctedheight = location.getAltitude().subtract(undu).subtract(interpolator.interpolate(e -> e.getHs()));
 
         // Temperature gradient [K/m]
         final T dTdH = interpolator.interpolate(e -> e.getModel(SeasonalModelType.DT).evaluate(dayOfYear)).multiply(0.001);
 
         // Specific humidity
         final T qv = interpolator.interpolate(e -> e.getModel(SeasonalModelType.QV).evaluate(dayOfYear)).multiply(0.001);
 
         // For the computation of the temperature and the pressure, we use
         // the standard ICAO atmosphere formulas.
 
         // Temperature [K]
         final T t0 = interpolator.interpolate(e -> e.getModel(SeasonalModelType.TEMPERATURE).evaluate(dayOfYear));
         final T temperature = correctedheight.multiply(dTdH).add(t0);
 
         // Pressure [hPa]
         final T p0       = interpolator.interpolate(e -> e.getModel(SeasonalModelType.PRESSURE).evaluate(dayOfYear));
         final T exponent = dTdH.multiply(R).reciprocal().multiply(G);
         final T pressure = FastMath.pow(correctedheight.multiply(dTdH.negate().divide(t0)).add(1), exponent).multiply(p0.multiply(0.01));
 
         // Water vapor pressure [hPa]
         final T e0 = pressure.multiply(qv.divide(qv.multiply(0.378).add(0.622 )));
 
         // mean temperature weighted with water vapor pressure
         final T tm = grid.hasModels(SeasonalModelType.TM) ?
                      interpolator.interpolate(e -> e.getModel(SeasonalModelType.TM).evaluate(dayOfYear)) :
                      date.getField().getZero().newInstance(Double.NaN);
 
         // water vapor decrease factor
         final T lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
                          interpolator.interpolate(e -> e.getModel(SeasonalModelType.LAMBDA).evaluate(dayOfYear)) :
                          date.getField().getZero().newInstance(Double.NaN);
 
         return new FieldPressureTemperatureHumidity<>(location.getAltitude(),
                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
                                                       temperature,
                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(e0),
                                                       tm,
                                                       lambda);
 
     }
 
     /** {@inheritDoc} */
     @Override
     public <T extends CalculusFieldElement<T>> FieldAzimuthalGradientCoefficients<T> getGradientCoefficients(final FieldGeodeticPoint<T> location,
                                                                                                              final FieldAbsoluteDate<T> date) {
 
         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
             // set up interpolation parameters
             final FieldCellInterpolator<T> interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude());
             final int                      dayOfYear    = date.getComponents(utc).getDate().getDayOfYear();
 
             return new FieldAzimuthalGradientCoefficients<>(interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_H).evaluate(dayOfYear)),
                                                             interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_H).evaluate(dayOfYear)),
                                                             interpolator.interpolate(e -> e.getModel(SeasonalModelType.GN_W).evaluate(dayOfYear)),
                                                             interpolator.interpolate(e -> e.getModel(SeasonalModelType.GE_W).evaluate(dayOfYear)));
         } else {
             return null;
         }
 
     }
 
 }