AbstractGrid.java

  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.troposphere.iturp834;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.analysis.interpolation.GridAxis;
  21. import org.hipparchus.util.FastMath;
  22. import org.orekit.bodies.FieldGeodeticPoint;
  23. import org.orekit.bodies.GeodeticPoint;
  24. import org.orekit.errors.OrekitException;
  25. import org.orekit.errors.OrekitMessages;
  26. import org.orekit.utils.units.Unit;

  28. import java.io.BufferedReader;
  29. import java.io.IOException;
  30. import java.io.InputStream;
  31. import java.io.InputStreamReader;
  32. import java.nio.charset.StandardCharsets;
  33. import java.util.regex.Pattern;

  35. /** Container for grid data.
  36.  * @author Luc Maisonobe
  37.  * @since 13.0
  38.  */
  39. abstract class AbstractGrid {

  41.     /** ITU-R P.834 data resources directory. */
  42.     private static final String ITU_R_P_834 = "/assets/org/orekit/ITU-R-P.834/";

  44.     /** Pattern for splitting fields. */
  45.     private static final Pattern SPLITTER = Pattern.compile("\\s+");

  47.     /** Minimum latitude (degrees). */
  48.     private static final double MIN_LAT = -90.0;

  50.     /** Maximum latitude (degrees). */
  51.     private static final double MAX_LAT =  90.0;

  53.     /** Grid step in latitude (degrees). */
  54.     private static final double STEP_LAT =   1.5;

  56.     /** Grid step in latitude (radians). */
  57.     private static final double STEP_LAT_RAD = FastMath.toRadians(STEP_LAT);

  59.     /** Minimum longitude (degrees). */
  60.     private static final double MIN_LON = -180.0;

  62.     /** Maximum longitude (degrees). */
  63.     private static final double MAX_LON =  180.0;

  65.     /** Grid step in longitude (degrees). */
  66.     private static final double STEP_LON =   1.5;

  68.     /** Grid step in longitude (radians). */
  69.     private static final double STEP_LON_RAD = FastMath.toRadians(STEP_LON);

  71.     /** Latitude grid axis. */
  72.     private final GridAxis latitudeAxis;

  74.     /** Longitude grid axis. */
  75.     private final GridAxis longitudeAxis;

  77.     /** Simple constructor.
  78.      */
  79.     protected AbstractGrid() {
  80.         latitudeAxis  = buildAxis(MIN_LAT, MAX_LAT, STEP_LAT);
  81.         longitudeAxis = buildAxis(MIN_LON, MAX_LON, STEP_LON);
  82.     }

  84.     /** Get latitude axis.
  85.      * @return latitude axis
  86.      */
  87.     public GridAxis getLatitudeAxis() {
  88.         return latitudeAxis;
  89.     }

  91.     /** Get longitude axis.
  92.      * @return longitude axis
  93.      */
  94.     public GridAxis getLongitudeAxis() {
  95.         return longitudeAxis;
  96.     }

  98.     /** Get cell size in latitude.
  99.      * @return cell size in latitude
  100.      */
  101.     public double getSizeLat() {
  102.         return STEP_LAT_RAD;
  103.     }

  105.     /** Get cell size in longitude.
  106.      * @return cell size in longitude
  107.      */
  108.     public double getSizeLon() {
  109.         return STEP_LON_RAD;
  110.     }

  112.     /** Get one raw cell.
  113.      * @param location point location on Earth
  114.      * @param rawData raww grid data
  115.      * @return raw cell
  116.      */
  117.     protected GridCell getRawCell(final GeodeticPoint location, final double[][] rawData) {

  119.         // locate the point
  120.         final int    southIndex    = latitudeAxis.interpolationIndex(location.getLatitude());
  121.         final double southLatitude = latitudeAxis.node(southIndex);
  122.         final int    westIndex     = longitudeAxis.interpolationIndex(location.getLongitude());
  123.         final double westLongitude = longitudeAxis.node(westIndex);

  125.         // build the cell
  126.         return new GridCell(location.getLatitude() - southLatitude,
  127.                             location.getLongitude() - westLongitude,
  128.                             STEP_LAT_RAD, STEP_LON_RAD,
  129.                             rawData[southIndex + 1][westIndex],
  130.                             rawData[southIndex][westIndex],
  131.                             rawData[southIndex][westIndex + 1],
  132.                             rawData[southIndex + 1][westIndex + 1]);

  134.     }

  136.     /** Get one raw cell.
  137.      * @param <T> type of the field elements
  138.      * @param location point location on Earth
  139.      * @param rawData raww grid data
  140.      * @return raw cell
  141.      */
  142.     protected <T extends CalculusFieldElement<T>> FieldGridCell<T> getRawCell(final FieldGeodeticPoint<T> location,
  143.                                                                               final double[][] rawData) {

  145.         // locate the point
  146.         final int    southIndex    = latitudeAxis.interpolationIndex(location.getLatitude().getReal());
  147.         final double southLatitude = latitudeAxis.node(southIndex);
  148.         final int    westIndex     = longitudeAxis.interpolationIndex(location.getLongitude().getReal());
  149.         final double westLongitude = longitudeAxis.node(westIndex);

  151.         // build the cell
  152.         final T zero = location.getAltitude().getField().getZero();
  153.         return new FieldGridCell<>(location.getLatitude().subtract(southLatitude),
  154.                                    location.getLongitude().subtract(westLongitude),
  155.                                    STEP_LAT_RAD, STEP_LON_RAD,
  156.                                    zero.newInstance(rawData[southIndex + 1][westIndex]),
  157.                                    zero.newInstance(rawData[southIndex][westIndex]),
  158.                                    zero.newInstance(rawData[southIndex][westIndex + 1]),
  159.                                    zero.newInstance(rawData[southIndex + 1][westIndex + 1]));

  161.     }

  163.     /** Get one cell.
  164.      * @param location     point location on Earth
  165.      * @param secondOfYear second of year
  166.      * @return cell at location
  167.      */
  168.     public abstract GridCell getCell(GeodeticPoint location, double secondOfYear);

  170.     /** Get one cell.
  171.      * @param <T>          type of the field elements
  172.      * @param location     point location on Earth
  173.      * @param secondOfYear second of year
  174.      * @return cell at location
  175.      */
  176.     public abstract <T extends CalculusFieldElement<T>> FieldGridCell<T> getCell(FieldGeodeticPoint<T> location,
  177.                                                                                  T secondOfYear);

  179.     /** Build a grid axis for interpolating within a table.
  180.      * @param min min angle in degrees (included)
  181.      * @param max max angle in degrees (included)
  182.      * @param step step between points
  183.      * @return grid axis
  184.      */
  185.     private static GridAxis buildAxis(final double min, final double max, final double step) {
  186.         final double[] grid = new double[(int) FastMath.rint((max - min) / step) + 1];
  187.         for (int i = 0; i < grid.length; i++) {
  188.             grid[i] = FastMath.toRadians(min + i * step);
  189.         }
  190.         return new GridAxis(grid, 2);
  191.     }

  193.     /** Parse interpolation table from a resource file.
  194.      * @param unit unit of values in resource file
  195.      * @param name name of the resource to parse
  196.      * @return parsed interpolation function
  197.      */
  198.     protected double[][] parse(final Unit unit, final String name) {

  200.         // parse the file
  201.         final double[][] values = new double[latitudeAxis.size()][longitudeAxis.size()];
  202.         try (InputStream       is     = ITURP834WeatherParametersProvider.class.getResourceAsStream(ITU_R_P_834 + name);
  203.              InputStreamReader isr    = is  == null ? null : new InputStreamReader(is, StandardCharsets.UTF_8);
  204.              BufferedReader    reader = isr == null ? null : new BufferedReader(isr)) {
  205.             if (reader == null) {
  206.                 // this should never happen with embedded data
  207.                 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, name);
  208.             }
  209.             for (int row = 0; row < latitudeAxis.size(); ++row) {

  211.                 // the ITU-files are sorted in decreasing latitude order, from +90° to -90°…
  212.                 final int latitudeIndex = latitudeAxis.size() - 1 - row;

  214.                 final String   line   = reader.readLine();
  215.                 final String[] fields = line == null ? new String[0] : SPLITTER.split(line.trim());
  216.                 if (fields.length != longitudeAxis.size()) {
  217.                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, row + 1, name, line);
  218.                 }

  220.                 // distribute points in longitude
  221.                 for (int col = 0; col < longitudeAxis.size(); ++col) {
  222.                     // files are between 0° and 360° in longitude, with last column (360°) equal to first column (0°)
  223.                     // our tables, on the other hand, use canonical longitudes between -180° and +180°
  224.                     // we have to redistribute indices
  225.                     // col =   0 → base longitude =   0.0° → fixed longitude =    0.0° → longitudeIndex = 120
  226.                     // col =   1 → base longitude =   1.5° → fixed longitude =    1.5° → longitudeIndex = 121
  227.                     // …
  228.                     // col = 119 → base longitude = 178.5° → fixed longitude =  178.5° → longitudeIndex = 239
  229.                     // col = 120 → base longitude = 180.0° → fixed longitude =  180.0° → longitudeIndex = 240
  230.                     // col = 121 → base longitude = 181.5° → fixed longitude = -178.5° → longitudeIndex =   1
  231.                     // …
  232.                     // col = 239 → base longitude = 358.5° → fixed longitude =   -1.5° → longitudeIndex = 119
  233.                     // col = 240 → base longitude = 360.0° → fixed longitude =    0.0° → longitudeIndex = 120
  234.                     final int longitudeIndex = col < 121 ? col + 120 : col - 120;
  235.                     values[latitudeIndex][longitudeIndex] = unit.toSI(Double.parseDouble(fields[col]));
  236.                 }

  238.                 // the loop above stored longitude 180° at index 240, but longitude -180° is missing
  239.                 values[latitudeIndex][0] = values[latitudeIndex][longitudeAxis.size() - 1];

  241.             }
  242.         } catch (IOException ioe) {
  243.             // this should never happen with the embedded data
  244.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, ioe);
  245.         }

  247.         // build the interpolating function
  248.         return values;

  250.     }

  252. }
Created with JaCoCo 0.8.12.202403310830