BoxAndSolarArraySpacecraft.java

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

  19. import java.util.ArrayList;
  20. import java.util.Collections;
  21. import java.util.List;
  22. import java.util.stream.Collectors;

  24. import org.hipparchus.CalculusFieldElement;
  25. import org.hipparchus.Field;
  26. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  27. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.Precision;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitInternalError;
  32. import org.orekit.forces.drag.DragSensitive;
  33. import org.orekit.forces.radiation.RadiationSensitive;
  34. import org.orekit.propagation.FieldSpacecraftState;
  35. import org.orekit.propagation.SpacecraftState;
  36. import org.orekit.utils.ExtendedPositionProvider;
  37. import org.orekit.utils.ParameterDriver;

  39. /** Class representing the features of a classical satellite with a convex body shape.
  40.  * <p>
  41.  * The body can be either a simple parallelepipedic box aligned with
  42.  * spacecraft axes or a set of panels defined by their area and normal vector.
  43.  * Some panels may be moving to model solar arrays (or antennas that could
  44.  * point anywhere). This should handle accurately most spacecraft shapes. This model
  45.  * does not take cast shadows into account.
  46.  * </p>
  47.  * <p>
  48.  * The lift component of the drag force can be optionally considered. It should
  49.  * probably only be used for reentry computation, with much denser atmosphere
  50.  * than in regular orbit propagation. The lift component is computed using a
  51.  * ratio of molecules that experience specular reflection instead of diffuse
  52.  * reflection (absorption followed by outgassing at negligible velocity).
  53.  * Without lift (i.e. when the lift ratio is set to 0), drag force is along
  54.  * atmosphere relative velocity. With lift (i.e. when the lift ratio is set to any
  55.  * value between 0 and 1), the drag force depends on both relative velocity direction
  56.  * and panels normal orientation. For a single panel, if the relative velocity is
  57.  * head-on (i.e. aligned with the panel normal), the force will be in the same
  58.  * direction with and without lift, but the magnitude with lift ratio set to 1.0 will
  59.  * be twice the magnitude with lift ratio set to 0.0 (because atmosphere molecules
  60.  * bounces backward at same velocity in case of specular reflection).
  61.  * </p>
  62.  * <p>
  63.  * Each {@link Panel panel} has its own set of radiation and drag coefficients. In
  64.  * orbit determination context, it would not be possible to estimate each panel
  65.  * individually, therefore {@link #getDragParametersDrivers()} returns a single
  66.  * {@link ParameterDriver parameter driver} representing a {@link DragSensitive#GLOBAL_DRAG_FACTOR
  67.  * global drag multiplicative factor} that applies to all panels drag coefficients
  68.  * and the {@link #getRadiationParametersDrivers()} returns a single
  69.  * {@link ParameterDriver parameter driver} representing a
  70.  * {@link RadiationSensitive#GLOBAL_RADIATION_FACTOR global radiation multiplicative factor}
  71.  * that applies to all panels radiation coefficients.
  72.  * </p>
  73.  *
  74.  * @author Luc Maisonobe
  75.  * @author Pascal Parraud
  76.  */
  77. public class BoxAndSolarArraySpacecraft implements RadiationSensitive, DragSensitive {

  79.     /** Parameters scaling factor.
  80.      * <p>
  81.      * We use a power of 2 to avoid numeric noise introduction
  82.      * in the multiplications/divisions sequences.
  83.      * </p>
  84.      */
  85.     private final double SCALE = FastMath.scalb(1.0, -3);

  87.     /** Driver for drag multiplicative factor parameter. */
  88.     private final ParameterDriver dragFactorParameterDriver;

  90.     /** Driver for radiation pressure multiplicative factor parameter. */
  91.     private final ParameterDriver radiationFactorParameterDriver;

  93.     /** Panels composing the spacecraft. */
  94.     private final List<Panel> panels;

  96.     /** Build a spacecraft model.
  97.      * @param panels panels composing the body, solar arrays and antennas
  98.      * (only the panels with strictly positive area will be stored)
  99.      * @since 12.0
  100.      */
  101.     public BoxAndSolarArraySpacecraft(final List<Panel> panels) {

  103.         try {
  104.             dragFactorParameterDriver      = new ParameterDriver(DragSensitive.GLOBAL_DRAG_FACTOR,
  105.                                                                  1.0, SCALE, 0.0, Double.POSITIVE_INFINITY);
  106.             radiationFactorParameterDriver = new ParameterDriver(RadiationSensitive.GLOBAL_RADIATION_FACTOR,
  107.                                                                  1.0, SCALE, 0.0, Double.POSITIVE_INFINITY);
  108.         } catch (OrekitException oe) {
  109.             // this should never happen
  110.             throw new OrekitInternalError(oe);
  111.         }

  113.         // remove spurious panels
  114.         this.panels = panels.stream().filter(p -> p.getArea() > 0).collect(Collectors.toList());

  116.     }

  118.     /** Build a spacecraft model with best lighting of solar array.
  119.      * <p>
  120.      * Solar arrays orientation will be such that at each time the Sun direction
  121.      * will always be in the solar array meridian plane defined by solar array
  122.      * rotation axis and solar array normal vector.
  123.      * </p>
  124.      * @param xLength length of the body along its X axis (m)
  125.      * @param yLength length of the body along its Y axis (m)
  126.      * @param zLength length of the body along its Z axis (m)
  127.      * @param sun sun model
  128.      * @param solarArrayArea area of the solar array (m²)
  129.      * @param solarArrayAxis solar array rotation axis in satellite frame
  130.      * @param dragCoeff drag coefficient (used only for drag)
  131.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  132.      * that will experience specular reflection when hitting spacecraft instead
  133.      * of experiencing diffuse reflection, hence producing lift)
  134.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  135.      * (used only for radiation pressure)
  136.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  137.      * (used only for radiation pressure)
  138.      * @since 12.0
  139.      */
  140.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength, final double zLength,
  141.                                       final ExtendedPositionProvider sun,
  142.                                       final double solarArrayArea, final Vector3D solarArrayAxis,
  143.                                       final double dragCoeff, final double liftRatio,
  144.                                       final double absorptionCoeff, final double reflectionCoeff) {
  145.         this(buildPanels(xLength, yLength, zLength,
  146.                          sun, solarArrayArea, solarArrayAxis,
  147.                          dragCoeff, liftRatio, absorptionCoeff, reflectionCoeff));
  148.     }

  150.     /** Get the panels composing the body.
  151.      * @return unmodifiable view of the panels composing the body
  152.      * @since 12.0
  153.      */
  154.     public List<Panel> getPanels() {
  155.         return Collections.unmodifiableList(panels);
  156.     }

  158.     /** {@inheritDoc} */
  159.     @Override
  160.     public List<ParameterDriver> getDragParametersDrivers() {
  161.         return Collections.singletonList(dragFactorParameterDriver);
  162.     }

  164.     /** {@inheritDoc} */
  165.     @Override
  166.     public List<ParameterDriver> getRadiationParametersDrivers() {
  167.         return Collections.singletonList(radiationFactorParameterDriver);
  168.     }

  170.     /** {@inheritDoc} */
  171.     @Override
  172.     public Vector3D dragAcceleration(final SpacecraftState state,
  173.                                      final double density, final Vector3D relativeVelocity,
  174.                                      final double[] parameters) {

  176.         final double dragFactor = parameters[0];

  178.         // relative velocity in spacecraft frame
  179.         final double   vNorm2 = relativeVelocity.getNormSq();
  180.         final double   vNorm  = FastMath.sqrt(vNorm2);
  181.         final Vector3D vDir   = state.getAttitude().getRotation().applyTo(relativeVelocity.scalarMultiply(1.0 / vNorm));
  182.         final double   coeff  = density * dragFactor * vNorm2 / (2.0 * state.getMass());

  184.         // panels contribution
  185.         Vector3D acceleration = Vector3D.ZERO;
  186.         for (final Panel panel : panels) {
  187.             Vector3D normal = panel.getNormal(state);
  188.             double dot = Vector3D.dotProduct(normal, vDir);
  189.             if (panel.isDoubleSided() && dot > 0) {
  190.                 // the flux comes from the back side
  191.                 normal = normal.negate();
  192.                 dot    = -dot;
  193.             }
  194.             if (dot < 0) {
  195.                 // the panel intercepts the incoming flux
  196.                 final double f         = coeff * panel.getDrag() * panel.getArea() * dot;
  197.                 final double liftRatio = panel.getLiftRatio();
  198.                 acceleration = new Vector3D(1,                                 acceleration,
  199.                                             (1 - liftRatio) * FastMath.abs(f), vDir,
  200.                                             liftRatio * f * 2,                 normal);
  201.             }
  202.         }

  204.         // convert back to inertial frame
  205.         return state.getAttitude().getRotation().applyInverseTo(acceleration);

  207.     }

  209.     /** {@inheritDoc} */
  210.     @Override
  211.     public <T extends CalculusFieldElement<T>> FieldVector3D<T>
  212.         dragAcceleration(final FieldSpacecraftState<T> state,
  213.                          final  T density, final FieldVector3D<T> relativeVelocity,
  214.                          final T[] parameters) {

  216.         final Field<T> field = state.getDate().getField();
  217.         final T dragFactor = parameters[0];

  219.         // relative velocity in spacecraft frame
  220.         final T                vNorm2 = relativeVelocity.getNormSq();
  221.         final T                vNorm  = FastMath.sqrt(vNorm2);
  222.         final FieldVector3D<T> vDir   = state.getAttitude().getRotation().applyTo(relativeVelocity.scalarMultiply(vNorm.reciprocal()));
  223.         final T                coeff  = density.multiply(dragFactor).multiply(vNorm2).divide(state.getMass().multiply(2.0));

  225.         // panels contribution
  226.         FieldVector3D<T> acceleration = FieldVector3D.getZero(field);
  227.         for (final Panel panel : panels) {
  228.             FieldVector3D<T> normal = panel.getNormal(state);
  229.             T dot = FieldVector3D.dotProduct(normal, vDir);
  230.             if (panel.isDoubleSided() && dot.getReal() > 0) {
  231.                 // the flux comes from the back side
  232.                 normal = normal.negate();
  233.                 dot    = dot.negate();
  234.             }
  235.             if (panel.isDoubleSided() || dot.getReal() < 0) {
  236.                 // the panel intercepts the incoming flux
  237.                 final T      f         = coeff.multiply(panel.getDrag() * panel.getArea()).multiply(dot);
  238.                 final double liftRatio = panel.getLiftRatio();
  239.                 acceleration = new FieldVector3D<>(field.getOne(),                         acceleration,
  240.                                                   FastMath.abs(f).multiply(1 - liftRatio), vDir,
  241.                                                   f.multiply(2 * liftRatio),               normal);
  242.             }
  243.         }

  245.         // convert back to inertial frame
  246.         return state.getAttitude().getRotation().applyInverseTo(acceleration);

  248.     }

  250.     /** {@inheritDoc} */
  251.     @Override
  252.     public Vector3D radiationPressureAcceleration(final SpacecraftState state,
  253.                                                   final Vector3D flux,
  254.                                                   final double[] parameters) {

  256.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  257.             // null illumination (we are probably in umbra)
  258.             return Vector3D.ZERO;
  259.         }

  261.         // radiation flux in spacecraft frame
  262.         final double   radiationFactor = parameters[0];
  263.         final Vector3D fluxSat         = state.getAttitude().getRotation().applyTo(flux).
  264.                                          scalarMultiply(radiationFactor);

  266.         // panels contribution
  267.         Vector3D force = Vector3D.ZERO;
  268.         for (final Panel panel : panels) {
  269.             Vector3D normal = panel.getNormal(state);
  270.             double dot = Vector3D.dotProduct(normal, fluxSat);
  271.             if (panel.isDoubleSided() && dot > 0) {
  272.                 // the flux comes from the back side
  273.                 normal = normal.negate();
  274.                 dot    = -dot;
  275.             }
  276.             if (dot < 0) {
  277.                 // the panel intercepts the incoming flux

  279.                 final double absorptionCoeff         = panel.getAbsorption();
  280.                 final double specularReflectionCoeff = panel.getReflection();
  281.                 final double diffuseReflectionCoeff  = 1 - (absorptionCoeff + specularReflectionCoeff);
  282.                 final double psr                     = fluxSat.getNorm();

  284.                 // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  285.                 // cos (phi) = -dot / (psr * area)
  286.                 // n         = panel / area
  287.                 // s         = -fluxSat / psr
  288.                 final double cN = 2 * panel.getArea() * dot * (diffuseReflectionCoeff / 3 - specularReflectionCoeff * dot / psr);
  289.                 final double cS = (panel.getArea() * dot / psr) * (specularReflectionCoeff - 1);
  290.                 force = new Vector3D(1, force, cN, normal, cS, fluxSat);

  292.             }
  293.         }

  295.         // convert to inertial frame
  296.         return state.getAttitude().getRotation().applyInverseTo(new Vector3D(1.0 / state.getMass(), force));

  298.     }

  300.     /** {@inheritDoc}
  301.      * <p>This method implements equation 8-44 from David A. Vallado's
  302.      * Fundamentals of Astrodynamics and Applications, third edition,
  303.      * 2007, Microcosm Press.</p>
  304.      */
  305.     @Override
  306.     public <T extends CalculusFieldElement<T>> FieldVector3D<T>
  307.         radiationPressureAcceleration(final FieldSpacecraftState<T> state,
  308.                                       final FieldVector3D<T> flux,
  309.                                       final T[] parameters) {

  311.         final Field<T> field = state.getDate().getField();
  312.         if (flux.getNormSq().getReal() < Precision.SAFE_MIN) {
  313.             // null illumination (we are probably in umbra)
  314.             return FieldVector3D.getZero(field);
  315.         }

  317.         // radiation flux in spacecraft frame
  318.         final T                radiationFactor = parameters[0];
  319.         final FieldVector3D<T> fluxSat         = state.getAttitude().getRotation().applyTo(flux).
  320.                                                  scalarMultiply(radiationFactor);

  322.         // panels contribution
  323.         FieldVector3D<T> force = FieldVector3D.getZero(field);
  324.         for (final Panel panel : panels) {
  325.             FieldVector3D<T> normal = panel.getNormal(state);
  326.             T dot = FieldVector3D.dotProduct(normal, fluxSat);
  327.             if (panel.isDoubleSided() && dot.getReal() > 0) {
  328.                 // the flux comes from the back side
  329.                 normal = normal.negate();
  330.                 dot    = dot.negate();
  331.             }
  332.             if (dot.getReal() < 0) {
  333.                 // the panel intercepts the incoming flux

  335.                 final double absorptionCoeff         = panel.getAbsorption();
  336.                 final double specularReflectionCoeff = panel.getReflection();
  337.                 final double diffuseReflectionCoeff  = 1 - (absorptionCoeff + specularReflectionCoeff);
  338.                 final T      psr                     = fluxSat.getNorm();

  340.                 // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  341.                 // cos (phi) = -dot / (psr * area)
  342.                 // n         = panel / area
  343.                 // s         = -fluxSat / psr
  344.                 final T cN = dot.multiply(-2 * panel.getArea()).multiply(dot.multiply(specularReflectionCoeff).divide(psr).subtract(diffuseReflectionCoeff / 3));
  345.                 final T cS = dot.multiply(panel.getArea()).multiply(specularReflectionCoeff - 1).divide(psr);
  346.                 force = new FieldVector3D<>(field.getOne(), force, cN, normal, cS, fluxSat);
  347.             }
  348.         }

  350.         // convert to inertial frame
  351.         return state.getAttitude().getRotation().applyInverseTo(new FieldVector3D<>(state.getMass().reciprocal(), force));

  353.     }

  355.     /** Build the panels of a simple parallelepipedic box.
  356.      * @param xLength length of the body along its X axis (m)
  357.      * @param yLength length of the body along its Y axis (m)
  358.      * @param zLength length of the body along its Z axis (m)
  359.      * @param drag drag coefficient
  360.      * @param liftRatio drag lift ratio (proportion between 0 and 1 of atmosphere modecules
  361.      * that will experience specular reflection when hitting spacecraft instead
  362.      * of experiencing diffuse reflection, hence producing lift)
  363.      * @param absorption radiation pressure absorption coefficient (between 0 and 1)
  364.      * @param reflection radiation pressure specular reflection coefficient (between 0 and 1)
  365.      * @return surface vectors array
  366.      * @since 12.0
  367.      */
  368.     public static List<Panel> buildBox(final double xLength, final double yLength, final double zLength,
  369.                                        final double drag, final double liftRatio,
  370.                                        final double absorption, final double reflection) {

  372.         final List<Panel> panels = new ArrayList<>(6);

  374.         // spacecraft body, composed of single-sided panels
  375.         panels.add(new FixedPanel(Vector3D.MINUS_I, yLength * zLength, false, drag, liftRatio, absorption, reflection));
  376.         panels.add(new FixedPanel(Vector3D.PLUS_I,  yLength * zLength, false, drag, liftRatio, absorption, reflection));
  377.         panels.add(new FixedPanel(Vector3D.MINUS_J, xLength * zLength, false, drag, liftRatio, absorption, reflection));
  378.         panels.add(new FixedPanel(Vector3D.PLUS_J,  xLength * zLength, false, drag, liftRatio, absorption, reflection));
  379.         panels.add(new FixedPanel(Vector3D.MINUS_K, xLength * yLength, false, drag, liftRatio, absorption, reflection));
  380.         panels.add(new FixedPanel(Vector3D.PLUS_K,  xLength * yLength, false, drag, liftRatio, absorption, reflection));

  382.         return panels;

  384.     }

  386.     /** Build the panels of a simple parallelepiped box plus one solar array panel.
  387.      * @param xLength length of the body along its X axis (m)
  388.      * @param yLength length of the body along its Y axis (m)
  389.      * @param zLength length of the body along its Z axis (m)
  390.      * @param sun sun model
  391.      * @param solarArrayArea area of the solar array (m²)
  392.      * @param solarArrayAxis solar array rotation axis in satellite frame
  393.      * @param drag drag coefficient
  394.      * @param liftRatio drag lift ratio (proportion between 0 and 1 of atmosphere modecules
  395.      * that will experience specular reflection when hitting spacecraft instead
  396.      * of experiencing diffuse reflection, hence producing lift)
  397.      * @param absorption radiation pressure absorption coefficient (between 0 and 1)
  398.      * @param reflection radiation pressure specular reflection coefficient (between 0 and 1)
  399.      * @return surface vectors array
  400.      * @since 12.0
  401.      */
  402.     public static List<Panel> buildPanels(final double xLength, final double yLength, final double zLength,
  403.                                           final ExtendedPositionProvider sun,
  404.                                           final double solarArrayArea, final Vector3D solarArrayAxis,
  405.                                           final double drag, final double liftRatio,
  406.                                           final double absorption, final double reflection) {

  408.         // spacecraft body
  409.         final List<Panel> panels = buildBox(xLength, yLength, zLength, drag, liftRatio, absorption, reflection);

  411.         // solar array
  412.         panels.add(new PointingPanel(solarArrayAxis, sun, solarArrayArea, drag, liftRatio, absorption, reflection));

  414.         return panels;

  416.     }

  418. }
Created with JaCoCo 0.8.12.202403310830