ModipGrid.java

/* Copyright 2002-2025 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.ionosphere.nequick;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.orekit.data.DataSource;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;

import java.io.BufferedReader;
import java.io.IOException;
import java.io.Reader;
import java.util.regex.Pattern;

/** Modified Dip Latitude grid.
 * <p>
 * The modip grid allows to estimate modip μ [deg] at a given point (φ,λ)
 * by interpolation of the relevant values contained in the support file.
 * </p>
 * <p>
 * The data is parsed from files sampling Earth in latitude and longitude. In
 * the case of Galileo-specific modip file, sampling step is 5° in latitude and
 * 10° in longitude, and extra rows/columns have been added to wrap around
 * anti-meridian and longitude and also "past" pole (adding extra points with
 * 180° longitude phasing offset) in latitude. In the case of ITU original
 * NeQuick 2 model, sampling step is 1° in latitude and 2° in longitude, without
 * any extra rows/columns. We therefore add the extra rows/columns upon parsing
 * the ITU file.
 * </p>
 * @author Bryan Cazabonne
 * @author Luc Maisonobe
 * @since 13.0
 */
class ModipGrid {

    /** Pattern for delimiting regular expressions. */
    private static final Pattern SEPARATOR = Pattern.compile("\\s+");

    /** Number of cells of modip grid in longitude (without wrapping). */
    private final int nbCellsLon;

    /** Cell size in longitude. */
    private final double sizeLon;

    /** Number of cells of modip grid in latitude (without wrapping). */
    private final int nbCellsLat;

    /** Cell size in latitude. */
    private final double sizeLat;

    /** Modip grid. */
    private final double[][] grid;

    /** Build a new modip grid.
     * @param nbCellsLon number of cells of modip grid in longitude (without wrapping)
     * @param nbCellsLat number of cells of modip grid in latitude (without wrapping)
     * @param source source of the grid file
     * @param wrappingAlreadyIncluded indicator for already included modip grid wrapping in resource file
     */
    ModipGrid(final int nbCellsLon, final int nbCellsLat, final DataSource source,
              final boolean wrappingAlreadyIncluded) {
        this.nbCellsLon = nbCellsLon;
        this.sizeLon    = MathUtils.TWO_PI / nbCellsLon;
        this.nbCellsLat = nbCellsLat;
        this.sizeLat    = FastMath.PI / nbCellsLat;
        this.grid       = new double[nbCellsLat + 3][nbCellsLon + 3];
        loadData(source, !wrappingAlreadyIncluded);
    }

    /** Compute modip for a location.
     * @param latitude latitude
     * @param longitude longitude
     * @return modip at specified location, in degrees
     */
    public double computeMODIP(final double latitude, final double longitude) {

        // index in latitude (note that Galileo document uses x for interpolation in latitude)
        // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
        final double x  = (latitude + MathUtils.SEMI_PI) / sizeLat + 1;
        if (x < 1) {
            return -90;
        }
        if (x >= nbCellsLat + 1) {
            return 90;
        }
        final int    i      = (int) FastMath.floor(x);
        final double deltaX = x - i;

        // index in longitude (note that Galileo document uses y for interpolation in longitude)
        // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
        double y  = (longitude + FastMath.PI) / sizeLon + 1;
        while (y < 1) {
            y += nbCellsLon;
        }
        while (y >= nbCellsLon + 1) {
            y -= nbCellsLon;
        }
        final int    j      = (int) FastMath.floor(y);
        final double deltaY = y - j;

        // zi coefficients (Eq. 12 and 13)
        final double z1 = interpolate(grid[i - 1][j - 1], grid[i][j - 1], grid[i + 1][j - 1], grid[i + 2][j - 1], deltaX);
        final double z2 = interpolate(grid[i - 1][j],     grid[i][j],     grid[i + 1][j],     grid[i + 2][j],     deltaX);
        final double z3 = interpolate(grid[i - 1][j + 1], grid[i][j + 1], grid[i + 1][j + 1], grid[i + 2][j + 1], deltaX);
        final double z4 = interpolate(grid[i - 1][j + 2], grid[i][j + 2], grid[i + 1][j + 2], grid[i + 2][j + 2], deltaX);

        // Modip (Ref Eq. 16)
        return interpolate(z1, z2, z3, z4, deltaY);

    }

    /**
     * This method provides a third order interpolation function
     * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
     *
     * @param z1 z1 coefficient
     * @param z2 z2 coefficient
     * @param z3 z3 coefficient
     * @param z4 z4 coefficient
     * @param x position
     * @return a third order interpolation
     */
    private double interpolate(final double z1, final double z2, final double z3, final double z4,
                               final double x) {

        if (x < 5e-11) {
            return z2;
        }

        final double delta = 2.0 * x - 1.0;
        final double g1 = z3 + z2;
        final double g2 = z3 - z2;
        final double g3 = z4 + z1;
        final double g4 = (z4 - z1) / 3.0;
        final double a0 = 9.0 * g1 - g3;
        final double a1 = 9.0 * g2 - g4;
        final double a2 = g3 - g1;
        final double a3 = g4 - g2;
        return 0.0625 * (a0 + delta * (a1 + delta * (a2 + delta * a3)));

    }

    /** Compute modip for a location.
     * @param <T> type of the field elements
     * @param latitude latitude
     * @param longitude longitude
     * @return modip at specified location
     */
    public <T extends CalculusFieldElement<T>> T computeMODIP(final T latitude, final T longitude) {

        // index in latitude (note that Galileo document uses x for interpolation in latitude)
        // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
        final T x  = latitude.add(MathUtils.SEMI_PI).divide(sizeLat).add(1);
        if (x.getReal() < 1) {
            return latitude.newInstance(-90);
        }
        if (x.getReal() >= nbCellsLat + 1) {
            return latitude.newInstance(90);
        }
        final int i      = (int) FastMath.floor(x.getReal());
        final T   deltaX = x.subtract(i);

        // index in longitude (note that Galileo document uses y for interpolation in longitude)
        // ensuring we have two points before (or at) index and two points after index for 3rd order interpolation
        T y  = longitude.add(FastMath.PI).divide(sizeLon).add(1);
        while (y.getReal() < 1) {
            y = y.add(nbCellsLon);
        }
        while (y.getReal() >= nbCellsLon + 1) {
            y = y.subtract(nbCellsLon);
        }
        final int j      = (int) FastMath.floor(y.getReal());
        final T   deltaY = y.subtract(j);

        // zi coefficients (Eq. 12 and 13)
        final T z1 = interpolate(x.newInstance(grid[i - 1][j - 1]), x.newInstance(grid[i][j - 1]),
                                 x.newInstance(grid[i + 1][j - 1]), x.newInstance(grid[i + 2][j - 1]),
                                 deltaX);
        final T z2 = interpolate(x.newInstance(grid[i - 1][j]),     x.newInstance(grid[i][j]),
                                 x.newInstance(grid[i + 1][j]),     x.newInstance(grid[i + 2][j]),
                                 deltaX);
        final T z3 = interpolate(x.newInstance(grid[i - 1][j + 1]), x.newInstance(grid[i][j + 1]),
                                 x.newInstance(grid[i + 1][j + 1]), x.newInstance(grid[i + 2][j + 1]),
                                 deltaX);
        final T z4 = interpolate(x.newInstance(grid[i - 1][j + 2]), x.newInstance(grid[i][j + 2]),
                                 x.newInstance(grid[i + 1][j + 2]), x.newInstance(grid[i + 2][j + 2]),
                                 deltaX);

        // Modip (Ref Eq. 16)
        return interpolate(z1, z2, z3, z4, deltaY);

    }

    /**
     * This method provides a third order interpolation function
     * as recommended in the reference document (Ref Eq. 128 to Eq. 138)
     *
     * @param <T> type of the field elements
     * @param z1 z1 coefficient
     * @param z2 z2 coefficient
     * @param z3 z3 coefficient
     * @param z4 z4 coefficient
     * @param x position
     * @return a third order interpolation
     */
    private <T extends CalculusFieldElement<T>> T interpolate(final T z1, final T z2, final T z3, final T z4,
                                                              final T x) {

        if (x.getReal() < 5e-11) {
            return z2;
        }

        final T delta = x.multiply(2.0).subtract(1.0);
        final T g1    = z3.add(z2);
        final T g2    = z3.subtract(z2);
        final T g3    = z4.add(z1);
        final T g4    = z4.subtract(z1).divide(3.0);
        final T a0    = g1.multiply(9.0).subtract(g3);
        final T a1    = g2.multiply(9.0).subtract(g4);
        final T a2    = g3.subtract(g1);
        final T a3    = g4.subtract(g2);
        return delta.multiply(a3).add(a2).multiply(delta).add(a1).multiply(delta).add(a0).multiply(0.0625);

    }

    /**
     * Load data from a stream.
     *
     * @param source grid source
     * @param addWrapping if true, wrapping should be added to loaded data
     */
    private void loadData(final DataSource source, final boolean addWrapping) {
        // if we must add wrapping, we must keep some empty rows and columns that will be filled later on
        final int first = addWrapping ? 1 : 0;

        // Open stream and parse data
        int lineNumber = 0;
        String line = null;
        int row = first;
        try (Reader         r  = source.getOpener().openReaderOnce();
             BufferedReader br = new BufferedReader(r)) {

            for (line = br.readLine(); line != null; line = br.readLine()) {
                ++lineNumber;

                // Read grid data
                final String[] fields = SEPARATOR.split(line.trim());
                if (fields.length == (addWrapping ? grid[row].length - 2 : grid[row].length)) {
                    // this should be a regular data line (i.e. not a header line)
                    for (int i = 0; i < fields.length; i++) {
                        grid[row][first + i] = Double.parseDouble(fields[i]);
                    }
                    ++row;
                }

            }

        } catch (IOException | NumberFormatException e) {
            throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, lineNumber, source.getName(), line);
        }

        if (addWrapping) {

            // we must add rows phased 180° in longitude after poles, for the sake of interpolation
            final double[] extraNorth = grid[0];
            final double[] refNorth   = grid[2];
            final double[] extraSouth = grid[grid.length - 1];
            final double[] refSouth   = grid[grid.length - 3];
            final int      deltaPhase = nbCellsLon / 2;
            for (int j = 1; j <= deltaPhase; j++) {
                extraNorth[j]              = refNorth[j + deltaPhase];
                extraNorth[j + deltaPhase] = refNorth[j];
                extraSouth[j]              = refSouth[j + deltaPhase];
                extraSouth[j + deltaPhase] = refSouth[j];
            }
            extraNorth[nbCellsLon + 1] = extraNorth[1];
            extraSouth[nbCellsLon + 1] = extraSouth[1];

            // we must copy columns to wrap around Earth in longitude
            // the first three columns must be identical to the last three columns
            // the columns 1 and gi.length - 2 (i.e. anti-meridian) are already
            // identical in the original resource file
            for (final double[] gi : grid) {
                gi[0] = gi[gi.length - 3];
                gi[gi.length - 1] = gi[2];
            }

            row++;

        }

        // Throw an exception if modip grid was not loaded properly
        if (row != grid.length) {
            throw new OrekitException(OrekitMessages.MODIP_GRID_NOT_LOADED, source.getName());
        }

    }

}