KnockeRediffusedForceModel.java

  1. /* Copyright 2002-2024 CS GROUP
  2.  * Licensed to CS GROUP (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.forces.radiation;

  18. import java.util.List;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  21. import org.hipparchus.analysis.polynomials.PolynomialsUtils;
  22. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  23. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  24. import org.hipparchus.geometry.euclidean.threed.Rotation;
  25. import org.hipparchus.geometry.euclidean.threed.RotationConvention;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.FieldSinCos;
  29. import org.hipparchus.util.MathUtils;
  30. import org.hipparchus.util.SinCos;
  31. import org.orekit.annotation.DefaultDataContext;
  32. import org.orekit.bodies.OneAxisEllipsoid;
  33. import org.orekit.data.DataContext;
  34. import org.orekit.forces.ForceModel;
  35. import org.orekit.frames.FieldStaticTransform;
  36. import org.orekit.frames.Frame;
  37. import org.orekit.frames.StaticTransform;
  38. import org.orekit.propagation.FieldSpacecraftState;
  39. import org.orekit.propagation.SpacecraftState;
  40. import org.orekit.time.AbsoluteDate;
  41. import org.orekit.time.FieldAbsoluteDate;
  42. import org.orekit.time.TimeScale;
  43. import org.orekit.utils.Constants;
  44. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  45. import org.orekit.utils.ParameterDriver;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  115.     /** General constructor.
  116.      * @param sun Sun model
  117.      * @param spacecraft the object physical and geometrical information
  118.      * @param equatorialRadius the Earth equatorial radius in m
  119.      * @param angularResolution angular resolution in rad
  120.      * @param utc the UTC time scale to define reference epoch
  121.      */
  122.     public KnockeRediffusedForceModel (final ExtendedPVCoordinatesProvider sun,
  123.                                        final RadiationSensitive spacecraft,
  124.                                        final double equatorialRadius,
  125.                                        final double angularResolution,
  126.                                        final TimeScale utc) {
  127.         this.sun               = sun;
  128.         this.spacecraft        = spacecraft;
  129.         this.equatorialRadius  = equatorialRadius;
  130.         this.angularResolution = angularResolution;
  131.         this.referenceEpoch    = new AbsoluteDate(1981, 12, 22, 0, 0, 0.0, utc);
  132.     }


  133.     /** {@inheritDoc} */
  134.     @Override
  135.     public boolean dependsOnPositionOnly() {
  136.         return false;
  137.     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

  273.     }


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

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

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

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

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

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

  297.     }

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

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

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

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

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

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

  322.     }


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

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

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

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

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

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

  346.     }

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

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

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


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

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

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


  368.     /** Compute elementary rediffused flux on satellite.
  369.      * @param state the current spacecraft state
  370.      * @param elementCenter the position of the considered area center
  371.      * @param sunPosition the position of the Sun in the spacecraft frame
  372.      * @param earth the Earth model
  373.      * @param elementArea the area of the current element
  374.      * @return the rediffused flux from considered element on the spacecraft
  375.      */
  376.     public Vector3D computeElementaryFlux(final SpacecraftState state,
  377.                                           final Vector3D elementCenter,
  378.                                           final Vector3D sunPosition,
  379.                                           final OneAxisEllipsoid earth,
  380.                                           final double elementArea) {

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

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

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

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

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

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

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

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

  397.             // Initialize albedo
  398.             double a = 0.0;

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

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

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

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

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

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

  416.         } else {

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

  420.     }


  421.     /** Compute elementary rediffused flux on satellite.
  422.      * @param state the current spacecraft state
  423.      * @param elementCenter the position of the considered area center
  424.      * @param sunPosition the position of the Sun in the spacecraft frame
  425.      * @param earth the Earth model
  426.      * @param elementArea the area of the current element
  427.      * @param <T> extends CalculusFieldElement
  428.      * @return the rediffused flux from considered element on the spacecraft
  429.      */
  430.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> computeElementaryFlux(final FieldSpacecraftState<T> state,
  431.                                                                                       final FieldVector3D<T> elementCenter,
  432.                                                                                       final FieldVector3D<T> sunPosition,
  433.                                                                                       final OneAxisEllipsoid earth,
  434.                                                                                       final T elementArea) {

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

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

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

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

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

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

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

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

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

  453.             // Initialize albedo
  454.             T a = zero;

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

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

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

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

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

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

  472.         } else {

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

  476.     }

  477. }