/* Copyright 2002-2017 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.forces.gravity.potential;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.text.ParseException;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.errors.OrekitParseException;
import org.orekit.time.DateComponents;
import org.orekit.utils.Constants;
/** Reader for the ICGEM gravity field format.
*
* <p>This format is used to describe the gravity field of EIGEN models
* published by the GFZ Potsdam since 2004. It is described in Franz
* Barthelmes and Christoph Förste paper: "the ICGEM-format".
* The 2006-02-28 version of this paper can be found <a
* href="http://op.gfz-potsdam.de/grace/results/grav/g005_ICGEM-Format.pdf">here</a>
* and the 2011-06-07 version of this paper can be found <a
* href="http://icgem.gfz-potsdam.de/ICGEM/documents/ICGEM-Format-2011.pdf">here</a>.
* These versions differ in time-dependent coefficients, which are linear-only prior
* to 2011 (up to eigen-5 model) and have also harmonic effects after that date
* (starting with eigen-6 model). Both versions are supported by the class.</p>
* <p>
* This reader uses relaxed check on the gravity constant key so any key ending
* in gravity_constant is accepted and not only earth_gravity_constant as specified
* in the previous documents. This allows to read also non Earth gravity fields
* as found in <a href="http://icgem.gfz-potsdam.de/ICGEM/ModelstabBodies.html">ICGEM
* - Gravity Field Models of other Celestial Bodies</a> page to be read.
* </p>
* <p>
* In order to simplify implementation, some design choices have been made: the
* reference date and the periods of harmonic pulsations are stored globally and
* not on a per-coefficient basis. This has the following implications:
* </p>
* <ul>
* <li>
* all time-stamped coefficients must share the same reference date, otherwise
* an error will be triggered during parsing,
* </li>
* <li>
* in order to avoid too large memory and CPU consumption, only a few different
* periods should appear in the file,
* </li>
* <li>
* only one occurrence of each coefficient may appear in the file, otherwise
* an error will be triggered during parsing. Multiple occurrences with different
* time stamps are forbidden (both because they correspond to a duplicated entry
* and because they define two different reference dates as per previous design
* choice).
* </li>
* </ul>
*
* <p> The proper way to use this class is to call the {@link GravityFieldFactory}
* which will determine which reader to use with the selected gravity field file.</p>
*
* @see GravityFieldFactory
* @author Luc Maisonobe
*/
public class ICGEMFormatReader extends PotentialCoefficientsReader {
/** Product type. */
private static final String PRODUCT_TYPE = "product_type";
/** Gravity field product type. */
private static final String GRAVITY_FIELD = "gravity_field";
/** Gravity constant marker. */
private static final String GRAVITY_CONSTANT = "gravity_constant";
/** Reference radius. */
private static final String REFERENCE_RADIUS = "radius";
/** Max degree. */
private static final String MAX_DEGREE = "max_degree";
/** Tide system indicator. */
private static final String TIDE_SYSTEM_INDICATOR = "tide_system";
/** Indicator value for zero-tide system. */
private static final String ZERO_TIDE = "zero_tide";
/** Indicator value for tide-free system. */
private static final String TIDE_FREE = "tide_free";
/** Indicator value for unknown tide system. */
private static final String TIDE_UNKNOWN = "unknown";
/** Normalization indicator. */
private static final String NORMALIZATION_INDICATOR = "norm";
/** Indicator value for normalized coefficients. */
private static final String NORMALIZED = "fully_normalized";
/** Indicator value for un-normalized coefficients. */
private static final String UNNORMALIZED = "unnormalized";
/** End of header marker. */
private static final String END_OF_HEADER = "end_of_head";
/** Gravity field coefficient. */
private static final String GFC = "gfc";
/** Time stamped gravity field coefficient. */
private static final String GFCT = "gfct";
/** Gravity field coefficient first time derivative. */
private static final String DOT = "dot";
/** Gravity field coefficient trend. */
private static final String TRND = "trnd";
/** Gravity field coefficient sine amplitude. */
private static final String ASIN = "asin";
/** Gravity field coefficient cosine amplitude. */
private static final String ACOS = "acos";
/** Tide system. */
private TideSystem tideSystem;
/** Indicator for normalized coefficients. */
private boolean normalized;
/** Reference date. */
private DateComponents referenceDate;
/** Secular trend of the cosine coefficients. */
private final List<List<Double>> cTrend;
/** Secular trend of the sine coefficients. */
private final List<List<Double>> sTrend;
/** Cosine part of the cosine coefficients pulsation. */
private final Map<Double, List<List<Double>>> cCos;
/** Sine part of the cosine coefficients pulsation. */
private final Map<Double, List<List<Double>>> cSin;
/** Cosine part of the sine coefficients pulsation. */
private final Map<Double, List<List<Double>>> sCos;
/** Sine part of the sine coefficients pulsation. */
private final Map<Double, List<List<Double>>> sSin;
/** Simple constructor.
* @param supportedNames regular expression for supported files names
* @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
*/
public ICGEMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
super(supportedNames, missingCoefficientsAllowed);
referenceDate = null;
cTrend = new ArrayList<List<Double>>();
sTrend = new ArrayList<List<Double>>();
cCos = new HashMap<Double, List<List<Double>>>();
cSin = new HashMap<Double, List<List<Double>>>();
sCos = new HashMap<Double, List<List<Double>>>();
sSin = new HashMap<Double, List<List<Double>>>();
}
/** {@inheritDoc} */
public void loadData(final InputStream input, final String name)
throws IOException, ParseException, OrekitException {
// reset the indicator before loading any data
setReadComplete(false);
referenceDate = null;
cTrend.clear();
sTrend.clear();
cCos.clear();
cSin.clear();
sCos.clear();
sSin.clear();
// by default, the field is normalized with unknown tide system
// (will be overridden later if non-default)
normalized = true;
tideSystem = TideSystem.UNKNOWN;
final BufferedReader r = new BufferedReader(new InputStreamReader(input, "UTF-8"));
boolean inHeader = true;
double[][] c = null;
double[][] s = null;
boolean okCoeffs = false;
int lineNumber = 0;
for (String line = r.readLine(); line != null; line = r.readLine()) {
try {
++lineNumber;
if (line.trim().length() == 0) {
continue;
}
final String[] tab = line.split("\\s+");
if (inHeader) {
if ((tab.length == 2) && PRODUCT_TYPE.equals(tab[0])) {
if (!GRAVITY_FIELD.equals(tab[1])) {
throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, line);
}
} else if ((tab.length == 2) && tab[0].endsWith(GRAVITY_CONSTANT)) {
setMu(parseDouble(tab[1]));
} else if ((tab.length == 2) && REFERENCE_RADIUS.equals(tab[0])) {
setAe(parseDouble(tab[1]));
} else if ((tab.length == 2) && MAX_DEGREE.equals(tab[0])) {
final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
final int order = FastMath.min(getMaxParseOrder(), degree);
c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
} else if ((tab.length == 2) && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
if (ZERO_TIDE.equals(tab[1])) {
tideSystem = TideSystem.ZERO_TIDE;
} else if (TIDE_FREE.equals(tab[1])) {
tideSystem = TideSystem.TIDE_FREE;
} else if (TIDE_UNKNOWN.equals(tab[1])) {
tideSystem = TideSystem.UNKNOWN;
} else {
throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, line);
}
} else if ((tab.length == 2) && NORMALIZATION_INDICATOR.equals(tab[0])) {
if (NORMALIZED.equals(tab[1])) {
normalized = true;
} else if (UNNORMALIZED.equals(tab[1])) {
normalized = false;
} else {
throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, line);
}
} else if ((tab.length == 2) && END_OF_HEADER.equals(tab[0])) {
inHeader = false;
}
} else {
if ((tab.length == 7 && GFC.equals(tab[0])) || (tab.length == 8 && GFCT.equals(tab[0]))) {
final int i = Integer.parseInt(tab[1]);
final int j = Integer.parseInt(tab[2]);
if (i < c.length && j < c[i].length) {
parseCoefficient(tab[3], c, i, j, "C", name);
parseCoefficient(tab[4], s, i, j, "S", name);
okCoeffs = true;
if (tab.length == 8) {
// check the reference date (format yyyymmdd)
final DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
Integer.parseInt(tab[7].substring(4, 6)),
Integer.parseInt(tab[7].substring(6, 8)));
if (referenceDate == null) {
// first reference found, store it
referenceDate = localRef;
} else if (!referenceDate.equals(localRef)) {
throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
referenceDate, localRef, name);
}
}
}
} else if (tab.length == 7 && (DOT.equals(tab[0]) || TRND.equals(tab[0]))) {
final int i = Integer.parseInt(tab[1]);
final int j = Integer.parseInt(tab[2]);
if (i < c.length && j < c[i].length) {
// store the secular trend coefficients
extendListOfLists(cTrend, i, j, 0.0);
extendListOfLists(sTrend, i, j, 0.0);
parseCoefficient(tab[3], cTrend, i, j, "Ctrend", name);
parseCoefficient(tab[4], sTrend, i, j, "Strend", name);
}
} else if (tab.length == 8 && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {
final int i = Integer.parseInt(tab[1]);
final int j = Integer.parseInt(tab[2]);
if (i < c.length && j < c[i].length) {
// extract arrays corresponding to period
final Double period = Double.valueOf(tab[7]);
if (!cCos.containsKey(period)) {
cCos.put(period, new ArrayList<List<Double>>());
cSin.put(period, new ArrayList<List<Double>>());
sCos.put(period, new ArrayList<List<Double>>());
sSin.put(period, new ArrayList<List<Double>>());
}
final List<List<Double>> cCosPeriod = cCos.get(period);
final List<List<Double>> cSinPeriod = cSin.get(period);
final List<List<Double>> sCosPeriod = sCos.get(period);
final List<List<Double>> sSinPeriod = sSin.get(period);
// store the pulsation coefficients
extendListOfLists(cCosPeriod, i, j, 0.0);
extendListOfLists(cSinPeriod, i, j, 0.0);
extendListOfLists(sCosPeriod, i, j, 0.0);
extendListOfLists(sSinPeriod, i, j, 0.0);
if (ACOS.equals(tab[0])) {
parseCoefficient(tab[3], cCosPeriod, i, j, "Ccos", name);
parseCoefficient(tab[4], sCosPeriod, i, j, "SCos", name);
} else {
parseCoefficient(tab[3], cSinPeriod, i, j, "Csin", name);
parseCoefficient(tab[4], sSinPeriod, i, j, "Ssin", name);
}
}
} else {
throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, line);
}
}
} catch (NumberFormatException nfe) {
final OrekitParseException pe = new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, line);
pe.initCause(nfe);
throw pe;
}
}
if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
// ensure at least the (0, 0) element is properly set
if (Precision.equals(c[0][0], 0.0, 0)) {
c[0][0] = 1.0;
}
}
if (Double.isNaN(getAe()) || Double.isNaN(getMu()) || !okCoeffs) {
String loaderName = getClass().getName();
loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
name, loaderName);
}
setRawCoefficients(normalized, c, s, name);
setTideSystem(tideSystem);
setReadComplete(true);
}
/** Get a provider for read spherical harmonics coefficients.
* <p>
* ICGEM fields do include time-dependent parts which are taken into account
* in the returned provider.
* </p>
* @param wantNormalized if true, the provider will provide normalized coefficients,
* otherwise it will provide un-normalized coefficients
* @param degree maximal degree
* @param order maximal order
* @return a new provider
* @exception OrekitException if the requested maximal degree or order exceeds the
* available degree or order or if no gravity field has read yet
* @since 6.0
*/
public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
final int degree, final int order)
throws OrekitException {
RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
if (cTrend.isEmpty() && cCos.isEmpty()) {
// there are no time-dependent coefficients
return provider;
}
if (!cTrend.isEmpty()) {
// add the secular trend layer
final double[][] cArrayTrend = toArray(cTrend);
final double[][] sArrayTrend = toArray(sTrend);
rescale(1.0 / Constants.JULIAN_YEAR, normalized, cArrayTrend, sArrayTrend, wantNormalized, cArrayTrend, sArrayTrend);
provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArrayTrend, sArrayTrend);
}
for (final Map.Entry<Double, List<List<Double>>> entry : cCos.entrySet()) {
final double period = entry.getKey();
// add the pulsating layer for the current period
final double[][] cArrayCos = toArray(cCos.get(period));
final double[][] sArrayCos = toArray(sCos.get(period));
final double[][] cArraySin = toArray(cSin.get(period));
final double[][] sArraySin = toArray(sSin.get(period));
rescale(1.0, normalized, cArrayCos, sArrayCos, wantNormalized, cArrayCos, sArrayCos);
rescale(1.0, normalized, cArraySin, sArraySin, wantNormalized, cArraySin, sArraySin);
provider = new PulsatingSphericalHarmonics(provider, period * Constants.JULIAN_YEAR,
cArrayCos, cArraySin, sArrayCos, sArraySin);
}
return provider;
}
}