1 /* Copyright 2002-2019 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (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
19 import java.io.BufferedReader;
20 import java.io.IOException;
21 import java.io.InputStream;
22 import java.io.InputStreamReader;
23 import java.util.Map;
24 import java.util.regex.Matcher;
25 import java.util.regex.Pattern;
26
27 import org.hipparchus.util.FastMath;
28 import org.orekit.errors.OrekitException;
29 import org.orekit.errors.OrekitMessages;
30
31 /** Reader for ocean tides files following the fes2004.dat format.
32 * @since 6.1
33 * @author Luc Maisonobe
34 */
35 public class FESCHatEpsilonReader extends OceanTidesReader {
36
37 /** Default pattern for fields with unknown type (non-space characters). */
38 private static final String UNKNOWN_TYPE_PATTERN = "\\S+";
39
40 /** Pattern for fields with integer type. */
41 private static final String INTEGER_TYPE_PATTERN = "[-+]?\\p{Digit}+";
42
43 /** Pattern for fields with real type. */
44 private static final String REAL_TYPE_PATTERN = "[-+]?(?:(?:\\p{Digit}+(?:\\.\\p{Digit}*)?)|(?:\\.\\p{Digit}+))(?:[eE][-+]?\\p{Digit}+)?";
45
46 /** Pattern for fields with Doodson number. */
47 private static final String DOODSON_TYPE_PATTERN = "\\p{Digit}{2,3}[.,]\\p{Digit}{3}";
48
49 /** Sea water fensity. */
50 private static final double RHO = 1025;
51
52 /** Gravitational constant (from IERS 2010, chapter 1). */
53 private static final double BIG_G = 6.67428e-11;
54
55 /** Earth mean gravity AT EQUATOR (from IERS 2010, chapter 1). */
56 private static final double GE = 9.7803278;
57
58 /** Scale of the CHat parameters. */
59 private final double scaleCHat;
60
61 /** Scale of the epsilon parameters. */
62 private final double scaleEpsilon;
63
64 /** Load deformation coefficients for ocean tides. */
65 private final OceanLoadDeformationCoefficients oldc;
66
67 /** Map for astronomical amplitudes. */
68 private final Map<Integer, Double> astronomicalAmplitudes;
69
70 /** Simple constructor.
71 * @param supportedNames regular expression for supported files names
72 * @param scaleCHat scale of the CHat parameters
73 * @param scaleEpsilon scale of the epsilon parameters
74 * @param oldc load deformation coefficients for ocean tides
75 * @param astronomicalAmplitudes map for astronomical amplitudes
76 * @see AstronomicalAmplitudeReader#getAstronomicalAmplitudesMap()
77 */
78 public FESCHatEpsilonReader(final String supportedNames,
79 final double scaleCHat, final double scaleEpsilon,
80 final OceanLoadDeformationCoefficients oldc,
81 final Map<Integer, Double> astronomicalAmplitudes) {
82 super(supportedNames);
83 this.scaleCHat = scaleCHat;
84 this.scaleEpsilon = scaleEpsilon;
85 this.oldc = oldc;
86 this.astronomicalAmplitudes = astronomicalAmplitudes;
87 }
88
89 /** {@inheritDoc} */
90 @Override
91 public void loadData(final InputStream input, final String name)
92 throws IOException {
93
94 // FES ocean tides models have the following form:
95 // Ocean tide model: FES2004 normalized model (fev. 2004) up to (100,100) in cm
96 // (long period from FES2002 up to (50,50) + equilibrium Om1/Om2, atmospheric tide NOT included)
97 // Doodson Darw n m Csin+ Ccos+ Csin- Ccos- C+ eps+ C- eps-
98 // 55.565 Om1 2 0 -0.540594 0.000000 0.000000 0.000000 0.5406 270.000 0.0000 0.000
99 // 55.575 Om2 2 0 -0.005218 0.000000 0.000000 0.000000 0.0052 270.000 0.0000 0.000
100 // 56.554 Sa 1 0 0.017233 0.000013 0.000000 0.000000 0.0172 89.957 0.0000 0.000
101 // 56.554 Sa 2 0 -0.046604 -0.000903 0.000000 0.000000 0.0466 268.890 0.0000 0.000
102 // 56.554 Sa 3 0 -0.000889 0.000049 0.000000 0.000000 0.0009 273.155 0.0000 0.000
103 final String[] fieldsPatterns = new String[] {
104 DOODSON_TYPE_PATTERN,
105 UNKNOWN_TYPE_PATTERN,
106 INTEGER_TYPE_PATTERN,
107 INTEGER_TYPE_PATTERN,
108 REAL_TYPE_PATTERN,
109 REAL_TYPE_PATTERN,
110 REAL_TYPE_PATTERN,
111 REAL_TYPE_PATTERN,
112 REAL_TYPE_PATTERN,
113 REAL_TYPE_PATTERN,
114 REAL_TYPE_PATTERN,
115 REAL_TYPE_PATTERN
116 };
117 final StringBuilder builder = new StringBuilder("^\\p{Space}*");
118 for (int i = 0; i < fieldsPatterns.length; ++i) {
119 builder.append("(");
120 builder.append(fieldsPatterns[i]);
121 builder.append(")");
122 builder.append((i < fieldsPatterns.length - 1) ? "\\p{Space}+" : "\\p{Space}*$");
123 }
124 final Pattern regularLinePattern = Pattern.compile(builder.toString());
125
126 final double commonFactor = 4 * FastMath.PI * BIG_G * RHO / GE;
127 final double[] kPrime = oldc.getCoefficients();
128
129 // parse the file
130 startParse(name);
131 final BufferedReader r = new BufferedReader(new InputStreamReader(input, "UTF-8"));
132 int lineNumber = 0;
133 for (String line = r.readLine(); line != null; line = r.readLine()) {
134 ++lineNumber;
135 final Matcher regularMatcher = regularLinePattern.matcher(line);
136 if (regularMatcher.matches()) {
137 // we have found a regular data line
138
139 // parse fields
140 final int doodson = Integer.parseInt(regularMatcher.group(1).replaceAll("[.,]", ""));
141 final int n = Integer.parseInt(regularMatcher.group(3));
142 final int m = Integer.parseInt(regularMatcher.group(4));
143
144 if (canAdd(n, m)) {
145
146 final double cHatPlus = scaleCHat * Double.parseDouble(regularMatcher.group(9));
147 final double ePlus = scaleEpsilon * Double.parseDouble(regularMatcher.group(10));
148 final double cHatMinus = scaleCHat * Double.parseDouble(regularMatcher.group(11));
149 final double eMinus = scaleEpsilon * Double.parseDouble(regularMatcher.group(12));
150
151 // compute bias from table 6.6
152 final double hf = astronomicalAmplitudes.containsKey(doodson) ? astronomicalAmplitudes.get(doodson) : 0.0;
153 final int cGamma = doodson / 100000;
154 final double chiF;
155 if (cGamma == 0) {
156 chiF = hf > 0 ? FastMath.PI : 0.0;
157 } else if (cGamma == 1) {
158 chiF = hf > 0 ? 0.5 * FastMath.PI : -0.5 * FastMath.PI;
159 } else if (cGamma == 2) {
160 chiF = hf > 0 ? 0.0 : FastMath.PI;
161 } else {
162 chiF = 0;
163 }
164
165 // compute reference gravity coefficients by converting height coefficients
166 // IERS conventions 2010, equation 6.21
167 if (n >= kPrime.length) {
168 throw new OrekitException(OrekitMessages.OCEAN_TIDE_LOAD_DEFORMATION_LIMITS,
169 kPrime.length - 1, n, name);
170 }
171 final double termFactor = (1 + kPrime[n]) / (2 * n + 1);
172
173 // an update on IERS conventions from 2012-08-10 states that for FES model:
174 // Note that, for zonal terms, FES2004 takes the approach to set
175 // the retrograde coefficients C-f,nO and S-f,n0 to zero and to double
176 // the prograde coefficients C+f,nO and S+f,n0. Therefore, after
177 // applying Equation (6.15), the ΔCn0 have the expected value but the
178 // ΔSn0 must be set to zero.
179 // (see ftp://tai.bipm.org/iers/convupdt/chapter6/icc6.pdf)
180 final double cPlus = commonFactor * termFactor * cHatPlus * FastMath.sin(ePlus + chiF);
181 final double sPlus = commonFactor * termFactor * cHatPlus * FastMath.cos(ePlus + chiF);
182 final double cMinus = (m == 0) ? 0.0 : commonFactor * termFactor * cHatMinus * FastMath.sin(eMinus + chiF);
183 final double sMinus = (m == 0) ? 0.0 : commonFactor * termFactor * cHatMinus * FastMath.cos(eMinus + chiF);
184
185 // store parsed fields
186 addWaveCoefficients(doodson, n, m, cPlus, sPlus, cMinus, sMinus, lineNumber, line);
187
188 }
189
190 }
191 }
192 endParse();
193
194 }
195
196 }