KnockeRediffusedForceModel.java

/* Copyright 2002-2024 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.forces.radiation;

import java.util.List;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.analysis.polynomials.PolynomialFunction;
import org.hipparchus.analysis.polynomials.PolynomialsUtils;
import org.hipparchus.geometry.euclidean.threed.FieldRotation;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.forces.ForceModel;
import org.orekit.frames.FieldStaticTransform;
import org.orekit.frames.Frame;
import org.orekit.frames.StaticTransform;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.utils.Constants;
import org.orekit.utils.ExtendedPVCoordinatesProvider;
import org.orekit.utils.ParameterDriver;

/** The Knocke Earth Albedo and IR emission force model.
 * <p>
 * This model is based on "EARTH RADIATION PRESSURE EFFECTS ON SATELLITES", 1988, by P. C. Knocke, J. C. Ries, and B. D. Tapley.
 * </p> <p>
 * This model represents the effects of radiation pressure coming from the Earth.
 * It considers Solar radiation which has been reflected by Earth (albedo) and Earth infrared emissions.
 * The planet is considered as a sphere and is divided into elementary areas.
 * Each elementary area is considered as a plane and emits radiation according to Lambert's law.
 * The flux the satellite receives is then equal to the sum of the elementary fluxes coming from Earth.
 * </p> <p>
 * The radiative model of the satellite, and its ability to diffuse, reflect  or absorb radiation is handled
 * by a {@link RadiationSensitive radiation sensitive model}.
 * </p> <p>
 * <b>Caution:</b> This model is only suitable for Earth. Using it with another central body is prone to error..
 * </p>
 *
 * @author Thomas Paulet
 * @since 10.3
 */
public class KnockeRediffusedForceModel implements ForceModel {

    /** 7Earth rotation around Sun pulsation in rad/sec. */
    private static final double EARTH_AROUND_SUN_PULSATION = MathUtils.TWO_PI / Constants.JULIAN_YEAR;

    /** Coefficient for solar irradiance computation. */
    private static final double ES_COEFF = 4.5606E-6;

    /** First coefficient for albedo computation. */
    private static final double A0 = 0.34;

    /** Second coefficient for albedo computation. */
    private static final double C0 = 0.;

    /** Third coefficient for albedo computation. */
    private static final double C1 = 0.10;

    /** Fourth coefficient for albedo computation. */
    private static final double C2 = 0.;

    /** Fifth coefficient for albedo computation. */
    private static final double A2 = 0.29;

    /** First coefficient for Earth emissivity computation. */
    private static final double E0 = 0.68;

    /** Second coefficient for Earth emissivity computation. */
    private static final double K0 = 0.;

    /** Third coefficient for Earth emissivity computation. */
    private static final double K1 = -0.07;

    /** Fourth coefficient for Earth emissivity computation. */
    private static final double K2 = 0.;

    /** Fifth coefficient for Earth emissivity computation. */
    private static final double E2 = -0.18;

    /** Sun model. */
    private final ExtendedPVCoordinatesProvider sun;

    /** Spacecraft. */
    private final RadiationSensitive spacecraft;

    /** Angular resolution for emissivity and albedo computation in rad. */
    private final double angularResolution;

    /** Earth equatorial radius in m. */
    private final double equatorialRadius;

    /** Reference date for periodic terms: December 22nd 1981.
     * Without more precision, the choice is to set it at midnight, UTC. */
    private final AbsoluteDate referenceEpoch;

    /** Default Constructor.
     * <p>This constructor uses the {@link DataContext#getDefault() default data context}</p>.
     * @param sun Sun model
     * @param spacecraft the object physical and geometrical information
     * @param equatorialRadius the Earth equatorial radius in m
     * @param angularResolution angular resolution in rad
     */
    @DefaultDataContext
    public KnockeRediffusedForceModel (final ExtendedPVCoordinatesProvider sun,
                                       final RadiationSensitive spacecraft,
                                       final double equatorialRadius,
                                       final double angularResolution) {

        this(sun, spacecraft, equatorialRadius, angularResolution, DataContext.getDefault().getTimeScales().getUTC());
    }

    /** General constructor.
     * @param sun Sun model
     * @param spacecraft the object physical and geometrical information
     * @param equatorialRadius the Earth equatorial radius in m
     * @param angularResolution angular resolution in rad
     * @param utc the UTC time scale to define reference epoch
     */
    public KnockeRediffusedForceModel (final ExtendedPVCoordinatesProvider sun,
                                       final RadiationSensitive spacecraft,
                                       final double equatorialRadius,
                                       final double angularResolution,
                                       final TimeScale utc) {
        this.sun               = sun;
        this.spacecraft        = spacecraft;
        this.equatorialRadius  = equatorialRadius;
        this.angularResolution = angularResolution;
        this.referenceEpoch    = new AbsoluteDate(1981, 12, 22, 0, 0, 0.0, utc);
    }


    /** {@inheritDoc} */
    @Override
    public boolean dependsOnPositionOnly() {
        return false;
    }

    /** {@inheritDoc} */
    @Override
    public Vector3D acceleration(final SpacecraftState s,
                                 final double[] parameters) {

        // Get date
        final AbsoluteDate date = s.getDate();

        // Get frame
        final Frame frame = s.getFrame();

        // Get satellite position
        final Vector3D satellitePosition = s.getPosition();

        // Get Sun position
        final Vector3D sunPosition = sun.getPosition(date, frame);

        // Get spherical Earth model
        final OneAxisEllipsoid earth = new OneAxisEllipsoid(equatorialRadius, 0.0, frame);

        // Project satellite on Earth as vector
        final Vector3D projectedToGround = satellitePosition.normalize().scalarMultiply(equatorialRadius);

        // Get elementary vector east for Earth browsing using rotations
        final Vector3D east = earth.transform(satellitePosition, frame, date).getEast();

        // Initialize rediffused flux with elementary flux coming from the circular area around the projected satellite
        final double centerArea = MathUtils.TWO_PI * equatorialRadius * equatorialRadius *
                                 (1.0 - FastMath.cos(angularResolution));
        Vector3D rediffusedFlux = computeElementaryFlux(s, projectedToGround, sunPosition, earth, centerArea);

        // Sectorize the part of Earth which is seen by the satellite into crown sectors with constant angular resolution
        for (double eastAxisOffset = 1.5 * angularResolution;
             eastAxisOffset < FastMath.asin(equatorialRadius / satellitePosition.getNorm());
             eastAxisOffset = eastAxisOffset + angularResolution) {

            // Build rotation transformations to get first crown elementary sector center
            final StaticTransform eastRotation = StaticTransform.of(date,
                                                          new Rotation(east, eastAxisOffset, RotationConvention.VECTOR_OPERATOR));

            // Get first elementary crown sector center
            final Vector3D firstCrownSectorCenter = eastRotation.transformPosition(projectedToGround);

            // Browse the entire crown
            for (double radialAxisOffset = 0.5 * angularResolution;
                 radialAxisOffset < MathUtils.TWO_PI;
                 radialAxisOffset = radialAxisOffset + angularResolution) {

                // Build rotation transformations to get elementary area center
                final StaticTransform radialRotation  = StaticTransform.of(date,
                                                                new Rotation(projectedToGround, radialAxisOffset, RotationConvention.VECTOR_OPERATOR));

                // Get current elementary crown sector center
                final Vector3D currentCenter = radialRotation.transformPosition(firstCrownSectorCenter);

                // Compute current elementary crown sector area, it results of the integration of an elementary crown sector
                // over the angular resolution
                final double sectorArea = equatorialRadius * equatorialRadius *
                                          2.0 * angularResolution * FastMath.sin(0.5 * angularResolution) * FastMath.sin(eastAxisOffset);

                // Add current sector contribution to total rediffused flux
                rediffusedFlux = rediffusedFlux.add(computeElementaryFlux(s, currentCenter, sunPosition, earth, sectorArea));
            }
        }

        return spacecraft.radiationPressureAcceleration(s, rediffusedFlux, parameters);
    }


    /** {@inheritDoc} */
    @Override
    public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
                                                                         final T[] parameters) {
        // Get date
        final FieldAbsoluteDate<T> date = s.getDate();

        // Get frame
        final Frame frame = s.getFrame();

        // Get zero
        final T zero = date.getField().getZero();

        // Get satellite position
        final FieldVector3D<T> satellitePosition = s.getPosition();

        // Get Sun position
        final FieldVector3D<T> sunPosition = sun.getPosition(date, frame);

        // Get spherical Earth model
        final OneAxisEllipsoid earth = new OneAxisEllipsoid(equatorialRadius, 0.0, frame);

        // Project satellite on Earth as vector
        final FieldVector3D<T> projectedToGround = satellitePosition.normalize().scalarMultiply(equatorialRadius);

        // Get elementary vector east for Earth browsing using rotations
        final FieldVector3D<T> east = earth.transform(satellitePosition, frame, date).getEast();

        // Initialize rediffused flux with elementary flux coming from the circular area around the projected satellite
        final T centerArea = zero.getPi().multiply(2.0).multiply(equatorialRadius).multiply(equatorialRadius).
                        multiply(1.0 - FastMath.cos(angularResolution));
        FieldVector3D<T> rediffusedFlux = computeElementaryFlux(s, projectedToGround, sunPosition, earth, centerArea);

        // Sectorize the part of Earth which is seen by the satellite into crown sectors with constant angular resolution
        for (double eastAxisOffset = 1.5 * angularResolution;
             eastAxisOffset < FastMath.asin(equatorialRadius / satellitePosition.getNorm().getReal());
             eastAxisOffset = eastAxisOffset + angularResolution) {

            // Build rotation transformations to get first crown elementary sector center
            final FieldStaticTransform<T> eastRotation = FieldStaticTransform.of(date,
                                                                        new FieldRotation<>(east,
                                                                                            zero.add(eastAxisOffset),
                                                                                            RotationConvention.VECTOR_OPERATOR));

            // Get first elementary crown sector center
            final FieldVector3D<T> firstCrownSectorCenter = eastRotation.transformPosition(projectedToGround);

            // Browse the entire crown
            for (double radialAxisOffset = 0.5 * angularResolution;
                 radialAxisOffset < MathUtils.TWO_PI;
                 radialAxisOffset = radialAxisOffset + angularResolution) {

                // Build rotation transformations to get elementary area center
                final FieldStaticTransform<T> radialRotation  = FieldStaticTransform.of(date,
                                                                               new FieldRotation<>(projectedToGround,
                                                                                                   zero.add(radialAxisOffset),
                                                                                                   RotationConvention.VECTOR_OPERATOR));

                // Get current elementary crown sector center
                final FieldVector3D<T> currentCenter = radialRotation.transformPosition(firstCrownSectorCenter);

                // Compute current elementary crown sector area, it results of the integration of an elementary crown sector
                // over the angular resolution
                final T sectorArea = zero.newInstance(equatorialRadius * equatorialRadius *
                        2.0 * angularResolution * FastMath.sin(0.5 * angularResolution) * FastMath.sin(eastAxisOffset));

                // Add current sector contribution to total rediffused flux
                rediffusedFlux = rediffusedFlux.add(computeElementaryFlux(s, currentCenter, sunPosition, earth, sectorArea));
            }
        }

        return spacecraft.radiationPressureAcceleration(s, rediffusedFlux, parameters);
    }


    /** {@inheritDoc} */
    @Override
    public List<ParameterDriver> getParametersDrivers() {
        return spacecraft.getRadiationParametersDrivers();
    }

    /** Compute Earth albedo.
     * Albedo value represents the fraction of solar radiative flux that is reflected by Earth.
     * Its value is in [0;1].
     * @param date the date
     * @param phi the equatorial latitude in rad
     * @return the albedo in [0;1]
     */
    public double computeAlbedo(final AbsoluteDate date, final double phi) {

        // Get duration since coefficient reference epoch
        final double deltaT = date.durationFrom(referenceEpoch);

        // Compute 1rst Legendre polynomial coeficient
        final SinCos sc = FastMath.sinCos(EARTH_AROUND_SUN_PULSATION * deltaT);
        final double A1 = C0 +
                          C1 * sc.cos() +
                          C2 * sc.sin();

        // Get 1rst and 2nd order Legendre polynomials
        final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
        final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);

        // Get latitude sinus
        final double sinPhi = FastMath.sin(phi);

        // Compute albedo
        return A0 +
               A1 * firstLegendrePolynomial.value(sinPhi) +
               A2 * secondLegendrePolynomial.value(sinPhi);

    }


    /** Compute Earth albedo.
     * Albedo value represents the fraction of solar radiative flux that is reflected by Earth.
     * Its value is in [0;1].
     * @param date the date
     * @param phi the equatorial latitude in rad
     * @param <T> extends CalculusFieldElement
     * @return the albedo in [0;1]
     */
    public <T extends CalculusFieldElement<T>> T computeAlbedo(final FieldAbsoluteDate<T> date, final T phi) {

        // Get duration since coefficient reference epoch
        final T deltaT = date.durationFrom(referenceEpoch);

        // Compute 1rst Legendre polynomial coeficient
        final FieldSinCos<T> sc = FastMath.sinCos(deltaT.multiply(EARTH_AROUND_SUN_PULSATION));
        final T A1 = sc.cos().multiply(C1).add(
                     sc.sin().multiply(C2)).add(C0);

        // Get 1rst and 2nd order Legendre polynomials
        final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
        final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);

        // Get latitude sinus
        final T sinPhi = FastMath.sin(phi);

        // Compute albedo
        return firstLegendrePolynomial.value(sinPhi).multiply(A1).add(
               secondLegendrePolynomial.value(sinPhi).multiply(A2)).add(A0);

    }

    /** Compute Earth emisivity.
     * Emissivity is used to compute the infrared flux that is emitted by Earth.
     * Its value is in [0;1].
     * @param date the date
     * @param phi the equatorial latitude in rad
     * @return the emissivity in [0;1]
     */
    public double computeEmissivity(final AbsoluteDate date, final double phi) {

        // Get duration since coefficient reference epoch
        final double deltaT = date.durationFrom(referenceEpoch);

        // Compute 1rst Legendre polynomial coefficient
        final SinCos sc = FastMath.sinCos(EARTH_AROUND_SUN_PULSATION * deltaT);
        final double E1 = K0 +
                          K1 * sc.cos() +
                          K2 * sc.sin();

        // Get 1rst and 2nd order Legendre polynomials
        final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
        final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);

        // Get latitude sinus
        final double sinPhi = FastMath.sin(phi);

        // Compute albedo
        return E0 +
               E1 * firstLegendrePolynomial.value(sinPhi) +
               E2 * secondLegendrePolynomial.value(sinPhi);

    }


    /** Compute Earth emisivity.
     * Emissivity is used to compute the infrared flux that is emitted by Earth.
     * Its value is in [0;1].
     * @param date the date
     * @param phi the equatorial latitude in rad
     * @param <T> extends CalculusFieldElement
     * @return the emissivity in [0;1]
     */
    public <T extends CalculusFieldElement<T>> T computeEmissivity(final FieldAbsoluteDate<T> date, final T phi) {

        // Get duration since coefficient reference epoch
        final T deltaT = date.durationFrom(referenceEpoch);

        // Compute 1rst Legendre polynomial coeficient
        final FieldSinCos<T> sc = FastMath.sinCos(deltaT.multiply(EARTH_AROUND_SUN_PULSATION));
        final T E1 = sc.cos().multiply(K1).add(
                     sc.sin().multiply(K2)).add(K0);

        // Get 1rst and 2nd order Legendre polynomials
        final PolynomialFunction firstLegendrePolynomial  = PolynomialsUtils.createLegendrePolynomial(1);
        final PolynomialFunction secondLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(2);

        // Get latitude sinus
        final T sinPhi = FastMath.sin(phi);

        // Compute albedo
        return firstLegendrePolynomial.value(sinPhi).multiply(E1).add(
               secondLegendrePolynomial.value(sinPhi).multiply(E2)).add(E0);

    }

    /** Compute total solar flux impacting Earth.
     * @param sunPosition the Sun position in an Earth centered frame
     * @return the total solar flux impacting Earth in J/m^3
     */
    public double computeSolarFlux(final Vector3D sunPosition) {

        // Compute Earth - Sun distance in UA
        final double earthSunDistance = sunPosition.getNorm() / Constants.JPL_SSD_ASTRONOMICAL_UNIT;

        // Compute Solar flux
        return ES_COEFF * Constants.SPEED_OF_LIGHT / (earthSunDistance * earthSunDistance);
    }


    /** Compute total solar flux impacting Earth.
     * @param sunPosition the Sun position in an Earth centered frame
     * @param <T> extends CalculusFieldElement
     * @return the total solar flux impacting Earth in J/m^3
     */
    public <T extends CalculusFieldElement<T>> T computeSolarFlux(final FieldVector3D<T> sunPosition) {

        // Compute Earth - Sun distance in UA
        final T earthSunDistance = sunPosition.getNorm().divide(Constants.JPL_SSD_ASTRONOMICAL_UNIT);

        // Compute Solar flux
        return earthSunDistance.multiply(earthSunDistance).reciprocal().multiply(ES_COEFF * Constants.SPEED_OF_LIGHT);
    }


    /** Compute elementary rediffused flux on satellite.
     * @param state the current spacecraft state
     * @param elementCenter the position of the considered area center
     * @param sunPosition the position of the Sun in the spacecraft frame
     * @param earth the Earth model
     * @param elementArea the area of the current element
     * @return the rediffused flux from considered element on the spacecraft
     */
    public Vector3D computeElementaryFlux(final SpacecraftState state,
                                          final Vector3D elementCenter,
                                          final Vector3D sunPosition,
                                          final OneAxisEllipsoid earth,
                                          final double elementArea) {

        // Get satellite position
        final Vector3D satellitePosition = state.getPosition();

        // Get current date
        final AbsoluteDate date = state.getDate();

        // Get frame
        final Frame frame = state.getFrame();

        // Get solar flux impacting Earth
        final double solarFlux = computeSolarFlux(sunPosition);

        // Get satellite viewing angle as seen from current elementary area
        final double alpha = Vector3D.angle(elementCenter, satellitePosition);

        // Check that satellite sees the current area
        if (FastMath.abs(alpha) < MathUtils.SEMI_PI) {

            // Get current elementary area center latitude
            final double currentLatitude = earth.transform(elementCenter, frame, date).getLatitude();

            // Compute Earth emissivity value
            final double e = computeEmissivity(date, currentLatitude);

            // Initialize albedo
            double a = 0.0;

            // Check if elementary area is in day light
            final double sunAngle = Vector3D.angle(elementCenter, sunPosition);

            if (FastMath.abs(sunAngle) < MathUtils.SEMI_PI) {
                // Elementary area is in day light, compute albedo value
                a = computeAlbedo(date, currentLatitude);
            }

            // Compute elementary area contribution to rediffused flux
            final double albedoAndIR = a * solarFlux * FastMath.cos(sunAngle) +
                                       e * solarFlux * 0.25;

            // Compute elementary area - satellite vector and distance
            final Vector3D r = satellitePosition.subtract(elementCenter);
            final double rNorm = r.getNorm();

            // Compute attenuated projected elemetary area vector
            final Vector3D projectedAreaVector = r.scalarMultiply(elementArea * FastMath.cos(alpha) /
                                                                 (FastMath.PI * rNorm * rNorm * rNorm));

            // Compute elementary radiation flux from current elementary area
            return projectedAreaVector.scalarMultiply(albedoAndIR / Constants.SPEED_OF_LIGHT);

        } else {

            // Spacecraft does not see the elementary area
            return new Vector3D(0.0, 0.0, 0.0);
        }

    }


    /** Compute elementary rediffused flux on satellite.
     * @param state the current spacecraft state
     * @param elementCenter the position of the considered area center
     * @param sunPosition the position of the Sun in the spacecraft frame
     * @param earth the Earth model
     * @param elementArea the area of the current element
     * @param <T> extends CalculusFieldElement
     * @return the rediffused flux from considered element on the spacecraft
     */
    public <T extends CalculusFieldElement<T>> FieldVector3D<T> computeElementaryFlux(final FieldSpacecraftState<T> state,
                                                                                      final FieldVector3D<T> elementCenter,
                                                                                      final FieldVector3D<T> sunPosition,
                                                                                      final OneAxisEllipsoid earth,
                                                                                      final T elementArea) {

        // Get satellite position
        final FieldVector3D<T> satellitePosition = state.getPosition();

        // Get current date
        final FieldAbsoluteDate<T> date = state.getDate();

        // Get frame
        final Frame frame = state.getFrame();

        // Get zero
        final T zero = date.getField().getZero();

        // Get solar flux impacting Earth
        final T solarFlux = computeSolarFlux(sunPosition);

        // Get satellite viewing angle as seen from current elementary area
        final T alpha = FieldVector3D.angle(elementCenter, satellitePosition);

        // Check that satellite sees the current area
        if (FastMath.abs(alpha).getReal() < MathUtils.SEMI_PI) {

            // Get current elementary area center latitude
            final T currentLatitude = earth.transform(elementCenter, frame, date).getLatitude();

            // Compute Earth emissivity value
            final T e = computeEmissivity(date, currentLatitude);

            // Initialize albedo
            T a = zero;

            // Check if elementary area is in day light
            final T sunAngle = FieldVector3D.angle(elementCenter, sunPosition);

            if (FastMath.abs(sunAngle).getReal() < MathUtils.SEMI_PI) {
                // Elementary area is in day light, compute albedo value
                a = computeAlbedo(date, currentLatitude);
            }

            // Compute elementary area contribution to rediffused flux
            final T albedoAndIR = a.multiply(solarFlux).multiply(FastMath.cos(sunAngle)).add(
                                  e.multiply(solarFlux).multiply(0.25));

            // Compute elementary area - satellite vector and distance
            final FieldVector3D<T> r = satellitePosition.subtract(elementCenter);
            final T rNorm = r.getNorm();

            // Compute attenuated projected elemetary area vector
            final FieldVector3D<T> projectedAreaVector = r.scalarMultiply(elementArea.multiply(FastMath.cos(alpha)).divide(
                                                                          rNorm.square().multiply(rNorm).multiply(zero.getPi())));

            // Compute elementary radiation flux from current elementary area
            return projectedAreaVector.scalarMultiply(albedoAndIR.divide(Constants.SPEED_OF_LIGHT));

        } else {

            // Spacecraft does not see the elementary area
            return new FieldVector3D<>(zero, zero, zero);
        }

    }

}