CssiSpaceWeatherData.java

/* Copyright 2020 Clément Jonglez
 * 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.
 * Clément Jonglez 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.atmosphere.data;

import java.util.List;
import java.util.stream.Collectors;

import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.AbstractSelfFeedingLoader;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.models.earth.atmosphere.DTM2000InputParameters;
import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherDataLoader.LineParameters;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeStamped;
import org.orekit.utils.Constants;
import org.orekit.utils.ImmutableTimeStampedCache;

/**
 * This class provides three-hourly and daily solar activity data needed by atmospheric
 * models: F107 solar flux, Ap and Kp indexes.
 * The {@link org.orekit.data.DataLoader} implementation and the parsing is handled by the class {@link CssiSpaceWeatherDataLoader}.
 * <p>
 * The data are retrieved through space weather files offered by AGI/CSSI on the AGI
 * <a href="ftp://ftp.agi.com/pub/DynamicEarthData/SpaceWeather-All-v1.2.txt">FTP</a> as
 * well as on the CelesTrack <a href="http://celestrak.com/SpaceData/">website</a>.
 * These files are updated several times a day by using several sources mentioned in the
 * <a href="http://celestrak.com/SpaceData/SpaceWx-format.php">Celestrak space
 * weather data documentation</a>.
 * </p>
 *
 * @author Clément Jonglez
 * @since 10.2
 */
public class CssiSpaceWeatherData extends AbstractSelfFeedingLoader
        implements DTM2000InputParameters, NRLMSISE00InputParameters {

    /** Default regular expression for supported names that works with all officially published files. */
    public static final String DEFAULT_SUPPORTED_NAMES = "^S(?:pace)?W(?:eather)?-(?:All)?.*\\.txt$";

    /** Serializable UID. */
    private static final long serialVersionUID = 4249411710645968978L;

    /** Size of the list. */
    private static final int N_NEIGHBORS = 2;

    /** Data set. */
    private final transient ImmutableTimeStampedCache<TimeStamped> data;

    /** UTC time scale. */
    private final TimeScale utc;

    /** First available date. */
    private final AbsoluteDate firstDate;

    /** Date of last data before the prediction starts. */
    private final AbsoluteDate lastObservedDate;

    /** Date of last daily prediction before the monthly prediction starts. */
    private final AbsoluteDate lastDailyPredictedDate;

    /** Last available date. */
    private final AbsoluteDate lastDate;

    /** Previous set of solar activity parameters. */
    private LineParameters previousParam;

    /** Current set of solar activity parameters. */
    private LineParameters nextParam;

    /**
     * Simple constructor. This constructor uses the default data context.
     * <p>
     * The original file names provided by AGI/CSSI are of the form:
     * SpaceWeather-All-v1.2.txt (AGI's ftp) or SW-Last5Years.txt (CelesTrak's website).
     * So a recommended regular expression for the supported names that works
     * with all published files is: {@link #DEFAULT_SUPPORTED_NAMES}.
     * </p>
     *
     * @param supportedNames regular expression for supported AGI/CSSI space weather files names
     */
    @DefaultDataContext
    public CssiSpaceWeatherData(final String supportedNames) {
        this(supportedNames, DataContext.getDefault().getDataProvidersManager(),
             DataContext.getDefault().getTimeScales().getUTC());
    }

    /**
     * Constructor that allows specifying the source of the CSSI space weather
     * file.
     *
     * @param supportedNames       regular expression for supported AGI/CSSI space weather files names
     * @param dataProvidersManager provides access to auxiliary data files.
     * @param utc                  UTC time scale.
     */
    public CssiSpaceWeatherData(final String supportedNames,
                                final DataProvidersManager dataProvidersManager,
                                final TimeScale utc) {
        super(supportedNames, dataProvidersManager);

        this.utc = utc;
        final CssiSpaceWeatherDataLoader loader =
            new CssiSpaceWeatherDataLoader(utc);
        this.feed(loader);
        data =
            new ImmutableTimeStampedCache<>(N_NEIGHBORS, loader.getDataSet());
        firstDate = loader.getMinDate();
        lastDate = loader.getMaxDate();
        lastObservedDate = loader.getLastObservedDate();
        lastDailyPredictedDate = loader.getLastDailyPredictedDate();
    }

    /** {@inheritDoc} */
    public AbsoluteDate getMinDate() {
        return firstDate;
    }

    /** {@inheritDoc} */
    public AbsoluteDate getMaxDate() {
        return lastDate;
    }

    /**
     * Find the data bracketing a specified date.
     *
     * @param date date to bracket
     */
    private void bracketDate(final AbsoluteDate date) {
        if (date.durationFrom(firstDate) < 0) {
            throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE,
                    date, firstDate, lastDate, firstDate.durationFrom(date));
        }
        if (date.durationFrom(lastDate) > 0) {
            throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER,
                    date, firstDate, lastDate, date.durationFrom(lastDate));
        }

        // don't search if the cached selection is fine
        if (previousParam != null && date.durationFrom(previousParam.getDate()) > 0 &&
                        date.durationFrom(nextParam.getDate()) <= 0) {
            return;
        }

        final List<TimeStamped> neigbors = data.getNeighbors(date).collect(Collectors.toList());
        previousParam = (LineParameters) neigbors.get(0);
        nextParam = (LineParameters) neigbors.get(1);
        if (previousParam.getDate().compareTo(date) > 0) { // TODO delete dead code
            /**
             * Throwing exception if neighbors are unbalanced because we are at the
             * beginning of the data set
             */
            throw new OrekitException(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE, date, firstDate, lastDate);
        }
    }

    /**
     * Performs a linear interpolation between two values The weights are computed
     * from the time delta between previous date, current date, next date.
     *
     * @param date          the current date
     * @param previousValue the value at previous date
     * @param nextValue     the value at next date
     * @return the value interpolated for the current date
     */
    private double getLinearInterpolation(final AbsoluteDate date, final double previousValue, final double nextValue) {
        // perform a linear interpolation
        final AbsoluteDate previousDate = previousParam.getDate();
        final AbsoluteDate currentDate = nextParam.getDate();
        final double dt = currentDate.durationFrom(previousDate);
        final double previousWeight = currentDate.durationFrom(date) / dt;
        final double nextWeight = date.durationFrom(previousDate) / dt;

        // returns the data interpolated at the date
        return previousValue * previousWeight + nextValue * nextWeight;
    }

    /** {@inheritDoc} */
    public double getInstantFlux(final AbsoluteDate date) {
        // Interpolating two neighboring daily fluxes
        // get the neighboring dates
        bracketDate(date);
        return getLinearInterpolation(date, previousParam.getF107Obs(), nextParam.getF107Obs());
    }

    /** {@inheritDoc} */
    public double getMeanFlux(final AbsoluteDate date) {
        return getAverageFlux(date);
    }

    /** {@inheritDoc} */
    public double getThreeHourlyKP(final AbsoluteDate date) {
        if (date.compareTo(lastObservedDate) <= 0) {
            /**
             * If observation data is available, it contains three-hourly data
             */
            bracketDate(date);
            final double hourOfDay = date.offsetFrom(previousParam.getDate(), utc) / 3600;
            int i_kp = (int) (hourOfDay / 3);
            if (i_kp >= 8) {
                /**
                 * hourOfDay can take the value 24.0 at midnight due to floating point precision
                 * when bracketing the dates or during a leap second because the hour of day is
                 * computed in UTC view
                 */
                i_kp = 7;
            }
            return previousParam.getThreeHourlyKp(i_kp);
        } else {
            /**
             * Only predictions are available, there are no three-hourly data
             */
            return get24HoursKp(date);
        }
    }

    /** {@inheritDoc} */
    public double get24HoursKp(final AbsoluteDate date) {
        if (date.compareTo(lastDailyPredictedDate) <= 0) {
            // Daily data is available, just taking the daily average
            bracketDate(date);
            return previousParam.getKpSum() / 8;
        } else {
            // Only monthly data is available, better interpolate between two months
            // get the neighboring dates
            bracketDate(date);
            return getLinearInterpolation(date, previousParam.getKpSum() / 8, nextParam.getKpSum() / 8);
        }
    }

    /** {@inheritDoc} */
    public double getDailyFlux(final AbsoluteDate date) {
        // Getting the value for the previous day
        return getDailyFluxOnDay(date.shiftedBy(-Constants.JULIAN_DAY));
    }

    /**
     * Gets the daily flux on the current day.
     *
     * @param date the current date
     * @return the daily F10.7 flux (observed)
     */
    private double getDailyFluxOnDay(final AbsoluteDate date) {
        if (date.compareTo(lastDailyPredictedDate) <= 0) {
            // Getting the value for the previous day
            bracketDate(date);
            return previousParam.getF107Obs();
        } else {
            // Only monthly data is available, better interpolate between two months
            // get the neighboring dates
            bracketDate(date);
            return getLinearInterpolation(date, previousParam.getF107Obs(), nextParam.getF107Obs());
        }
    }

    /** {@inheritDoc} */
    public double getAverageFlux(final AbsoluteDate date) {
        if (date.compareTo(lastDailyPredictedDate) <= 0) {
            bracketDate(date);
            return previousParam.getCtr81Obs();
        } else {
            // Only monthly data is available, better interpolate between two months
            // get the neighboring dates
            bracketDate(date);
            return getLinearInterpolation(date, previousParam.getCtr81Obs(), nextParam.getCtr81Obs());
        }
    }

    /** {@inheritDoc} */
    public double[] getAp(final AbsoluteDate date) {
        final double[] apArray = new double[7];
        apArray[0] = getDailyAp(date);
        apArray[1] = getThreeHourlyAp(date);
        apArray[2] = getThreeHourlyAp(date.shiftedBy(-3.0 * 3600.0));
        apArray[3] = getThreeHourlyAp(date.shiftedBy(-6.0 * 3600.0));
        apArray[4] = getThreeHourlyAp(date.shiftedBy(-9.0 * 3600.0));
        apArray[5] = get24HoursAverageAp(date.shiftedBy(-12.0 * 3600.0));
        apArray[6] = get24HoursAverageAp(date.shiftedBy(-36.0 * 3600.0));
        return apArray;
    }

    /**
     * Gets the value of the three-hourly Ap index for the given date.
     *
     * @param date the current date
     * @return the current three-hourly Ap index
     */
    private double getThreeHourlyAp(final AbsoluteDate date) {
        if (date.compareTo(lastObservedDate.shiftedBy(Constants.JULIAN_DAY)) < 0) {
            /**
             * If observation data is available, it contains three-hourly data.
             */
            bracketDate(date);
            final double hourOfDay = date.offsetFrom(previousParam.getDate(), utc) / 3600;
            int i_ap = (int) (hourOfDay / 3);
            if (i_ap >= 8) {
                /**
                 * hourOfDay can take the value 24.0 at midnight due to floating point precision
                 * when bracketing the dates or during a leap second because the hour of day is
                 * computed in UTC view
                 */
                i_ap = 7;
            }
            return previousParam.getThreeHourlyAp(i_ap);
        } else {
            /**
             * Only predictions are available, there are no three-hourly data
             */
            return getDailyAp(date);
        }
    }

    /**
     * Gets the running average of the 8 three-hourly Ap indices prior to current
     * time If three-hourly data is available, the result is different than
     * getDailyAp.
     *
     * @param date the current date
     * @return the 24 hours running average of the Ap index
     */
    private double get24HoursAverageAp(final AbsoluteDate date) {
        if (date.compareTo(lastDailyPredictedDate) <= 0) {
            // Computing running mean
            double apSum = 0.0;
            for (int i = 0; i < 8; i++) {
                apSum += getThreeHourlyAp(date.shiftedBy(-3.0 * 3600 * i));
            }
            return apSum / 8;
        } else {
            /**
             * Only monthly predictions are available, no need to compute the average from
             * three hourly data
             */
            return getDailyAp(date);
        }
    }

    /**
     * Get the daily Ap index for the given date.
     *
     * @param date the current date
     * @return the daily Ap index
     */
    private double getDailyAp(final AbsoluteDate date) {
        if (date.compareTo(lastDailyPredictedDate) <= 0) {
            // Daily data is available, just taking the daily average
            bracketDate(date);
            return previousParam.getApAvg();
        } else {
            // Only monthly data is available, better interpolate between two months
            // get the neighboring dates
            bracketDate(date);
            return getLinearInterpolation(date, previousParam.getApAvg(), nextParam.getApAvg());
        }
    }

    public String getSupportedNames() {
        return super.getSupportedNames();
    }
}