FESCHatEpsilonReader.java

  1. /* Copyright 2002-2020 CS GROUP
  2.  * Licensed to CS GROUP (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.forces.gravity.potential;

  18. import java.io.BufferedReader;
  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.InputStreamReader;
  22. import java.nio.charset.StandardCharsets;
  23. import java.util.Map;
  24. import java.util.regex.Matcher;
  25. import java.util.regex.Pattern;

  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.SinCos;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;

  30. /** Reader for ocean tides files following the fes2004.dat format.
  31.  * @since 6.1
  32.  * @author Luc Maisonobe
  33.  */
  34. public class FESCHatEpsilonReader extends OceanTidesReader {

  35.     /** Default pattern for fields with unknown type (non-space characters). */
  36.     private static final String  UNKNOWN_TYPE_PATTERN = "\\S+";

  37.     /** Pattern for fields with integer type. */
  38.     private static final String  INTEGER_TYPE_PATTERN = "[-+]?\\p{Digit}+";

  39.     /** Pattern for fields with real type. */
  40.     private static final String  REAL_TYPE_PATTERN = "[-+]?(?:(?:\\p{Digit}+(?:\\.\\p{Digit}*)?)|(?:\\.\\p{Digit}+))(?:[eE][-+]?\\p{Digit}+)?";

  41.     /** Pattern for fields with Doodson number. */
  42.     private static final String  DOODSON_TYPE_PATTERN = "\\p{Digit}{2,3}[.,]\\p{Digit}{3}";

  43.     /** Pattern for regular data. */
  44.     private static final Pattern PATTERN = Pattern.compile("[.,]");

  45.     /** Sea water fensity. */
  46.     private static final double RHO   = 1025;

  47.     /** Gravitational constant (from IERS 2010, chapter 1). */
  48.     private static final double BIG_G = 6.67428e-11;

  49.     /** Earth mean gravity AT EQUATOR (from IERS 2010, chapter 1). */
  50.     private static final double GE    = 9.7803278;

  51.     /** Scale of the CHat parameters. */
  52.     private final double scaleCHat;

  53.     /** Scale of the epsilon parameters. */
  54.     private final double scaleEpsilon;

  55.     /** Load deformation coefficients for ocean tides. */
  56.     private final OceanLoadDeformationCoefficients oldc;

  57.     /** Map for astronomical amplitudes. */
  58.     private final Map<Integer, Double> astronomicalAmplitudes;

  59.     /** Simple constructor.
  60.      * @param supportedNames regular expression for supported files names
  61.      * @param scaleCHat scale of the CHat parameters
  62.      * @param scaleEpsilon scale of the epsilon parameters
  63.      * @param oldc load deformation coefficients for ocean tides
  64.      * @param astronomicalAmplitudes map for astronomical amplitudes
  65.      * @see AstronomicalAmplitudeReader#getAstronomicalAmplitudesMap()
  66.      */
  67.     public FESCHatEpsilonReader(final String supportedNames,
  68.                                 final double scaleCHat, final double scaleEpsilon,
  69.                                 final OceanLoadDeformationCoefficients oldc,
  70.                                 final Map<Integer, Double> astronomicalAmplitudes) {
  71.         super(supportedNames);
  72.         this.scaleCHat              = scaleCHat;
  73.         this.scaleEpsilon           = scaleEpsilon;
  74.         this.oldc                   = oldc;
  75.         this.astronomicalAmplitudes = astronomicalAmplitudes;
  76.     }

  77.     /** {@inheritDoc} */
  78.     @Override
  79.     public void loadData(final InputStream input, final String name)
  80.         throws IOException {

  81.         // FES ocean tides models have the following form:
  82.         //   Ocean tide model: FES2004 normalized model (fev. 2004) up to (100,100) in cm
  83.         //   (long period from FES2002 up to (50,50) + equilibrium Om1/Om2, atmospheric tide NOT included)
  84.         //   Doodson Darw  n   m    Csin+     Ccos+       Csin-     Ccos-       C+   eps+      C-   eps-
  85.         //    55.565 Om1   2   0 -0.540594  0.000000    0.000000  0.000000   0.5406 270.000 0.0000   0.000
  86.         //    55.575 Om2   2   0 -0.005218  0.000000    0.000000  0.000000   0.0052 270.000 0.0000   0.000
  87.         //    56.554 Sa    1   0  0.017233  0.000013    0.000000  0.000000   0.0172  89.957 0.0000   0.000
  88.         //    56.554 Sa    2   0 -0.046604 -0.000903    0.000000  0.000000   0.0466 268.890 0.0000   0.000
  89.         //    56.554 Sa    3   0 -0.000889  0.000049    0.000000  0.000000   0.0009 273.155 0.0000   0.000
  90.         final String[] fieldsPatterns = new String[] {
  91.             DOODSON_TYPE_PATTERN,
  92.             UNKNOWN_TYPE_PATTERN,
  93.             INTEGER_TYPE_PATTERN,
  94.             INTEGER_TYPE_PATTERN,
  95.             REAL_TYPE_PATTERN,
  96.             REAL_TYPE_PATTERN,
  97.             REAL_TYPE_PATTERN,
  98.             REAL_TYPE_PATTERN,
  99.             REAL_TYPE_PATTERN,
  100.             REAL_TYPE_PATTERN,
  101.             REAL_TYPE_PATTERN,
  102.             REAL_TYPE_PATTERN
  103.         };
  104.         final StringBuilder builder = new StringBuilder("^\\p{Space}*");
  105.         for (int i = 0; i < fieldsPatterns.length; ++i) {
  106.             builder.append("(");
  107.             builder.append(fieldsPatterns[i]);
  108.             builder.append(")");
  109.             builder.append((i < fieldsPatterns.length - 1) ? "\\p{Space}+" : "\\p{Space}*$");
  110.         }
  111.         final Pattern regularLinePattern = Pattern.compile(builder.toString());

  112.         final double commonFactor = 4 * FastMath.PI * BIG_G * RHO / GE;
  113.         final double[] kPrime = oldc.getCoefficients();

  114.         // parse the file
  115.         startParse(name);
  116.         try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
  117.             int lineNumber      = 0;
  118.             for (String line = r.readLine(); line != null; line = r.readLine()) {
  119.                 ++lineNumber;
  120.                 final Matcher regularMatcher = regularLinePattern.matcher(line);
  121.                 if (regularMatcher.matches()) {
  122.                     // we have found a regular data line

  123.                     // parse fields
  124.                     final int doodson = Integer.parseInt(PATTERN.matcher(regularMatcher.group(1)).replaceAll(""));
  125.                     final int n       = Integer.parseInt(regularMatcher.group(3));
  126.                     final int m       = Integer.parseInt(regularMatcher.group(4));

  127.                     if (canAdd(n, m)) {

  128.                         final double cHatPlus  = scaleCHat    * Double.parseDouble(regularMatcher.group(9));
  129.                         final double ePlus     = scaleEpsilon * Double.parseDouble(regularMatcher.group(10));
  130.                         final double cHatMinus = scaleCHat    * Double.parseDouble(regularMatcher.group(11));
  131.                         final double eMinus    = scaleEpsilon * Double.parseDouble(regularMatcher.group(12));

  132.                         // compute bias from table 6.6
  133.                         final double hf = astronomicalAmplitudes.getOrDefault(doodson, 0.0);
  134.                         final int cGamma = doodson / 100000;
  135.                         final double chiF;
  136.                         if (cGamma == 0) {
  137.                             chiF = hf > 0 ? FastMath.PI : 0.0;
  138.                         } else if (cGamma == 1) {
  139.                             chiF = hf > 0 ? 0.5 * FastMath.PI : -0.5 * FastMath.PI;
  140.                         } else if (cGamma == 2) {
  141.                             chiF = hf > 0 ? 0.0 : FastMath.PI;
  142.                         } else {
  143.                             chiF = 0;
  144.                         }

  145.                         // compute reference gravity coefficients by converting height coefficients
  146.                         // IERS conventions 2010, equation 6.21
  147.                         if (n >= kPrime.length) {
  148.                             throw new OrekitException(OrekitMessages.OCEAN_TIDE_LOAD_DEFORMATION_LIMITS,
  149.                                                       kPrime.length - 1, n, name);
  150.                         }
  151.                         final double termFactor = (1 + kPrime[n]) / (2 * n + 1);

  152.                         // an update on IERS conventions from 2012-08-10 states that for FES model:
  153.                         //      Note that, for zonal terms, FES2004 takes the approach to set
  154.                         //      the retrograde coefficients C-f,nO and S-f,n0 to zero and to double
  155.                         //      the prograde coefficients C+f,nO and S+f,n0. Therefore, after
  156.                         //      applying Equation (6.15), the ΔCn0 have the expected value but the
  157.                         //      ΔSn0 must be set to zero.
  158.                         // (see ftp://tai.bipm.org/iers/convupdt/chapter6/icc6.pdf)
  159.                         final SinCos scP    = FastMath.sinCos(ePlus  + chiF);
  160.                         final SinCos scM    = FastMath.sinCos(eMinus + chiF);
  161.                         final double cPlus  =                  commonFactor * termFactor * cHatPlus  * scP.sin();
  162.                         final double sPlus  =                  commonFactor * termFactor * cHatPlus  * scP.cos();
  163.                         final double cMinus = (m == 0) ? 0.0 : commonFactor * termFactor * cHatMinus * scM.sin();
  164.                         final double sMinus = (m == 0) ? 0.0 : commonFactor * termFactor * cHatMinus * scM.cos();

  165.                         // store parsed fields
  166.                         addWaveCoefficients(doodson, n, m, cPlus,  sPlus, cMinus, sMinus, lineNumber, line);

  167.                     }

  168.                 }
  169.             }
  170.         }
  171.         endParse();

  172.     }

  173. }