BoxAndSolarArraySpacecraft.java

  1. /* Copyright 2002-2017 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

  18. import java.util.ArrayList;
  19. import java.util.List;

  20. import org.hipparchus.Field;
  21. import org.hipparchus.RealFieldElement;
  22. import org.hipparchus.analysis.differentiation.DSFactory;
  23. import org.hipparchus.analysis.differentiation.DerivativeStructure;
  24. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  25. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  26. import org.hipparchus.geometry.euclidean.threed.Rotation;
  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.errors.OrekitMessages;
  33. import org.orekit.forces.drag.DragSensitive;
  34. import org.orekit.forces.radiation.RadiationSensitive;
  35. import org.orekit.frames.Frame;
  36. import org.orekit.time.AbsoluteDate;
  37. import org.orekit.time.FieldAbsoluteDate;
  38. import org.orekit.utils.PVCoordinatesProvider;
  39. import org.orekit.utils.ParameterDriver;

  40. /** Class representing the features of a classical satellite
  41.  * with a convex body shape and rotating flat solar arrays.
  42.  * <p>
  43.  * The body can be either a simple parallelepipedic box aligned with
  44.  * spacecraft axes or a set of facets defined by their area and normal vector.
  45.  * This should handle accurately most spacecraft shapes.
  46.  * </p>
  47.  * <p>
  48.  * The solar array rotation with respect to satellite body can be either
  49.  * the best lighting orientation (i.e. Sun exactly in solar array meridian
  50.  * plane defined by solar array rotation axis and solar array normal vector)
  51.  * or a rotation evolving linearly according to a start position and an
  52.  * angular rate (which can be set to 0 for non-rotating panels, which may
  53.  * occur in special modes or during contingencies).
  54.  * </p>
  55.  * <p>
  56.  * The lift component of the drag force can be optionally considered. It should
  57.  * probably only be used for reentry computation, with much denser atmosphere
  58.  * than in regular orbit propagation. The lift component is computed using a
  59.  * ratio of molecules that experience specular reflection instead of diffuse
  60.  * reflection (absorption followed by outgassing at negligible velocity).
  61.  * Without lift (i.e. when the lift ratio is set to 0), drag force is along
  62.  * atmosphere relative velocity. With lift (i.e. when the lift ratio is set to any
  63.  * value between 0 and 1), the drag force depends on both relative velocity direction
  64.  * and facets normal orientation. For a single panel, if the relative velocity is
  65.  * head-on (i.e. aligned with the panel normal), the force will be in the same
  66.  * direction with and without lift, but the magnitude with lift ratio set to 1.0 will
  67.  * be twice the magnitude with lift ratio set to 0.0 (because atmosphere molecules
  68.  * bounces backward at same velocity in case of specular reflection).
  69.  * </p>
  70.  * <p>
  71.  * This model does not take cast shadow between body and solar array into account.
  72.  * </p>
  73.  *
  74.  * @author Luc Maisonobe
  75.  * @author Pascal Parraud
  76.  */
  77. public class BoxAndSolarArraySpacecraft implements RadiationSensitive, DragSensitive {

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

  85.     /** Driver for drag coefficient parameter. */
  86.     private final ParameterDriver dragParameterDriver;

  87.     /** Driver for lift ratio parameter (may be null is lift is ignored). */
  88.     private final ParameterDriver liftParameterDriver;

  89.     /** Driver for radiation pressure absorption coefficient parameter. */
  90.     private final ParameterDriver absorptionParameterDriver;

  91.     /** Driver for radiation pressure reflection coefficient parameter. */
  92.     private final ParameterDriver reflectionParameterDriver;

  93.     /** Surface vectors for body facets. */
  94.     private final List<Facet> facets;

  95.     /** Solar array area (m²). */
  96.     private final double solarArrayArea;

  97.     /** Reference date for linear rotation angle (may be null). */
  98.     private final AbsoluteDate referenceDate;

  99.     /** Rotation rate for linear rotation angle. */
  100.     private final double rotationRate;

  101.     /** Solar array reference axis in spacecraft frame (may be null). */
  102.     private final Vector3D saX;

  103.     /** Solar array third axis in spacecraft frame (may be null). */
  104.     private final Vector3D saY;

  105.     /** Solar array rotation axis in spacecraft frame. */
  106.     private final Vector3D saZ;

  107.     /** Sun model. */
  108.     private final PVCoordinatesProvider sun;

  109.     /** Factory for the DerivativeStructure instances. */
  110.     private final DSFactory factory;

  111.     /** Build a spacecraft model with best lighting of solar array.
  112.      * <p>
  113.      * This constructor builds an instance that completely ignores lift
  114.      * in atmospheric drag (the value of lift coefficient is set to zero,
  115.      * and there are no {@link ParameterDriver drivers} to change it).
  116.      * </p>
  117.      * <p>
  118.      * Solar arrays orientation will be such that at each time the Sun direction
  119.      * will always be in the solar array meridian plane defined by solar array
  120.      * rotation axis and solar array normal vector.
  121.      * </p>
  122.      * @param xLength length of the body along its X axis (m)
  123.      * @param yLength length of the body along its Y axis (m)
  124.      * @param zLength length of the body along its Z axis (m)
  125.      * @param sun sun model
  126.      * @param solarArrayArea area of the solar array (m²)
  127.      * @param solarArrayAxis solar array rotation axis in satellite frame
  128.      * @param dragCoeff drag coefficient (used only for drag)
  129.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  130.      * (used only for radiation pressure)
  131.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  132.      * (used only for radiation pressure)
  133.      */
  134.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  135.                                       final double zLength,
  136.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  137.                                       final Vector3D solarArrayAxis,
  138.                                       final double dragCoeff,
  139.                                       final double absorptionCoeff,
  140.                                       final double reflectionCoeff) {
  141.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  142.              dragCoeff, 0.0, false,
  143.              absorptionCoeff, reflectionCoeff);
  144.     }

  145.     /** Build a spacecraft model with best lighting of solar array.
  146.      * <p>
  147.      * Solar arrays orientation will be such that at each time the Sun direction
  148.      * will always be in the solar array meridian plane defined by solar array
  149.      * rotation axis and solar array normal vector.
  150.      * </p>
  151.      * @param xLength length of the body along its X axis (m)
  152.      * @param yLength length of the body along its Y axis (m)
  153.      * @param zLength length of the body along its Z axis (m)
  154.      * @param sun sun model
  155.      * @param solarArrayArea area of the solar array (m²)
  156.      * @param solarArrayAxis solar array rotation axis in satellite frame
  157.      * @param dragCoeff drag coefficient (used only for drag)
  158.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  159.      * that will experience specular reflection when hitting spacecraft instead
  160.      * of experiencing diffuse reflection, hence producing lift)
  161.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  162.      * (used only for radiation pressure)
  163.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  164.      * (used only for radiation pressure)
  165.      * @since 9.0
  166.      */
  167.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  168.                                       final double zLength,
  169.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  170.                                       final Vector3D solarArrayAxis,
  171.                                       final double dragCoeff, final double liftRatio,
  172.                                       final double absorptionCoeff,
  173.                                       final double reflectionCoeff) {
  174.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  175.              dragCoeff, liftRatio,
  176.              absorptionCoeff, reflectionCoeff);
  177.     }

  178.     /** Build a spacecraft model with best lighting of solar array.
  179.      * <p>
  180.      * The spacecraft body is described by an array of surface vectors. Each facet of
  181.      * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
  182.      * and whose norm is the surface area in m².
  183.      * </p>
  184.      * <p>
  185.      * Solar arrays orientation will be such that at each time the Sun direction
  186.      * will always be in the solar array meridian plane defined by solar array
  187.      * rotation axis and solar array normal vector.
  188.      * </p>
  189.      * @param facets body facets (only the facets with strictly positive area will be stored)
  190.      * @param sun sun model
  191.      * @param solarArrayArea area of the solar array (m²)
  192.      * @param solarArrayAxis solar array rotation axis in satellite frame
  193.      * @param dragCoeff drag coefficient (used only for drag)
  194.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  195.      * (used only for radiation pressure)
  196.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  197.      * (used only for radiation pressure)
  198.      */
  199.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  200.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  201.                                       final Vector3D solarArrayAxis,
  202.                                       final double dragCoeff,
  203.                                       final double absorptionCoeff,
  204.                                       final double reflectionCoeff) {
  205.         this(facets, sun, solarArrayArea, solarArrayAxis,
  206.              dragCoeff, 0.0, false,
  207.              absorptionCoeff, reflectionCoeff);
  208.     }

  209.     /** Build a spacecraft model with best lighting of solar array.
  210.      * <p>
  211.      * The spacecraft body is described by an array of surface vectors. Each facet of
  212.      * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
  213.      * and whose norm is the surface area in m².
  214.      * </p>
  215.      * <p>
  216.      * Solar arrays orientation will be such that at each time the Sun direction
  217.      * will always be in the solar array meridian plane defined by solar array
  218.      * rotation axis and solar array normal vector.
  219.      * </p>
  220.      * @param facets body facets (only the facets with strictly positive area will be stored)
  221.      * @param sun sun model
  222.      * @param solarArrayArea area of the solar array (m²)
  223.      * @param solarArrayAxis solar array rotation axis in satellite frame
  224.      * @param dragCoeff drag coefficient (used only for drag)
  225.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  226.      * that will experience specular reflection when hitting spacecraft instead
  227.      * of experiencing diffuse reflection, hence producing lift)
  228.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  229.      * (used only for radiation pressure)
  230.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  231.      * (used only for radiation pressure)
  232.      * @since 9.0
  233.      */
  234.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  235.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  236.                                       final Vector3D solarArrayAxis,
  237.                                       final double dragCoeff, final double liftRatio,
  238.                                       final double absorptionCoeff,
  239.                                       final double reflectionCoeff) {
  240.         this(facets, sun, solarArrayArea, solarArrayAxis,
  241.              dragCoeff, liftRatio, true,
  242.              absorptionCoeff, reflectionCoeff);
  243.     }

  244.     /** Build a spacecraft model with best lighting of solar array.
  245.      * <p>
  246.      * The spacecraft body is described by an array of surface vectors. Each facet of
  247.      * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
  248.      * and whose norm is the surface area in m².
  249.      * </p>
  250.      * <p>
  251.      * Solar arrays orientation will be such that at each time the Sun direction
  252.      * will always be in the solar array meridian plane defined by solar array
  253.      * rotation axis and solar array normal vector.
  254.      * </p>
  255.      * @param facets body facets (only the facets with strictly positive area will be stored)
  256.      * @param sun sun model
  257.      * @param solarArrayArea area of the solar array (m²)
  258.      * @param solarArrayAxis solar array rotation axis in satellite frame
  259.      * @param dragCoeff drag coefficient (used only for drag)
  260.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  261.      * that will experience specular reflection when hitting spacecraft instead
  262.      * of experiencing diffuse reflection, hence producing lift)
  263.      * @param useLift if true, lift should be used, otherwise it is completely ignored
  264.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  265.      * (used only for radiation pressure)
  266.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  267.      * (used only for radiation pressure)
  268.      * @since 9.0
  269.      */
  270.     private BoxAndSolarArraySpacecraft(final Facet[] facets,
  271.                                        final PVCoordinatesProvider sun, final double solarArrayArea,
  272.                                        final Vector3D solarArrayAxis,
  273.                                        final double dragCoeff, final double liftRatio, final boolean useLift,
  274.                                        final double absorptionCoeff,
  275.                                        final double reflectionCoeff) {

  276.         // drag
  277.         dragParameterDriver = buildDragParameterDriver(dragCoeff);
  278.         liftParameterDriver = useLift ? buildLiftParameterDriver(liftRatio) : null;

  279.         // radiation pressure
  280.         absorptionParameterDriver = buildAbsorptionParameterDriver(absorptionCoeff);
  281.         reflectionParameterDriver = buildReflectionParameterDriver(reflectionCoeff);

  282.         factory = new DSFactory(1, 1);

  283.         this.facets = filter(facets);

  284.         this.sun            = sun;
  285.         this.solarArrayArea = solarArrayArea;
  286.         this.referenceDate  = null;
  287.         this.rotationRate   = 0;

  288.         this.saZ = solarArrayAxis.normalize();
  289.         this.saY = null;
  290.         this.saX = null;

  291.     }

  292.     /** Build a spacecraft model with linear rotation of solar array.
  293.      * <p>
  294.      * Solar arrays orientation will be a regular rotation from the
  295.      * reference orientation at reference date and using a constant
  296.      * rotation rate.
  297.      * </p>
  298.      * @param xLength length of the body along its X axis (m)
  299.      * @param yLength length of the body along its Y axis (m)
  300.      * @param zLength length of the body along its Z axis (m)
  301.      * @param sun sun model
  302.      * @param solarArrayArea area of the solar array (m²)
  303.      * @param solarArrayAxis solar array rotation axis in satellite frame
  304.      * @param referenceDate reference date for the solar array rotation
  305.      * @param referenceNormal direction of the solar array normal at reference date
  306.      * in spacecraft frame
  307.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  308.      * @param dragCoeff drag coefficient (used only for drag)
  309.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  310.      * (used only for radiation pressure)
  311.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  312.      * (used only for radiation pressure)
  313.      */
  314.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  315.                                       final double zLength,
  316.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  317.                                       final Vector3D solarArrayAxis,
  318.                                       final AbsoluteDate referenceDate,
  319.                                       final Vector3D referenceNormal,
  320.                                       final double rotationRate,
  321.                                       final double dragCoeff,
  322.                                       final double absorptionCoeff,
  323.                                       final double reflectionCoeff) {
  324.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  325.              referenceDate, referenceNormal, rotationRate,
  326.              dragCoeff, 0.0, false,
  327.              absorptionCoeff, reflectionCoeff);
  328.     }

  329.     /** Build a spacecraft model with linear rotation of solar array.
  330.      * <p>
  331.      * Solar arrays orientation will be a regular rotation from the
  332.      * reference orientation at reference date and using a constant
  333.      * rotation rate.
  334.      * </p>
  335.      * @param xLength length of the body along its X axis (m)
  336.      * @param yLength length of the body along its Y axis (m)
  337.      * @param zLength length of the body along its Z axis (m)
  338.      * @param sun sun model
  339.      * @param solarArrayArea area of the solar array (m²)
  340.      * @param solarArrayAxis solar array rotation axis in satellite frame
  341.      * @param referenceDate reference date for the solar array rotation
  342.      * @param referenceNormal direction of the solar array normal at reference date
  343.      * in spacecraft frame
  344.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  345.      * @param dragCoeff drag coefficient (used only for drag)
  346.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  347.      * that will experience specular reflection when hitting spacecraft instead
  348.      * of experiencing diffuse reflection, hence producing lift)
  349.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  350.      * (used only for radiation pressure)
  351.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  352.      * (used only for radiation pressure)
  353.      * @since 9.0
  354.      */
  355.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  356.                                       final double zLength,
  357.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  358.                                       final Vector3D solarArrayAxis,
  359.                                       final AbsoluteDate referenceDate,
  360.                                       final Vector3D referenceNormal,
  361.                                       final double rotationRate,
  362.                                       final double dragCoeff, final double liftRatio,
  363.                                       final double absorptionCoeff,
  364.                                       final double reflectionCoeff) {
  365.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  366.              referenceDate, referenceNormal, rotationRate,
  367.              dragCoeff, liftRatio, true,
  368.              absorptionCoeff, reflectionCoeff);
  369.     }

  370.     /** Build a spacecraft model with linear rotation of solar array.
  371.      * <p>
  372.      * The spacecraft body is described by an array of surface vectors. Each facet of
  373.      * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
  374.      * and whose norm is the surface area in m².
  375.      * </p>
  376.      * <p>
  377.      * Solar arrays orientation will be a regular rotation from the
  378.      * reference orientation at reference date and using a constant
  379.      * rotation rate.
  380.      * </p>
  381.      * @param facets body facets (only the facets with strictly positive area will be stored)
  382.      * @param sun sun model
  383.      * @param solarArrayArea area of the solar array (m²)
  384.      * @param solarArrayAxis solar array rotation axis in satellite frame
  385.      * @param referenceDate reference date for the solar array rotation
  386.      * @param referenceNormal direction of the solar array normal at reference date
  387.      * in spacecraft frame
  388.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  389.      * @param dragCoeff drag coefficient (used only for drag)
  390.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  391.      * (used only for radiation pressure)
  392.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  393.      * (used only for radiation pressure)
  394.      */
  395.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  396.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  397.                                       final Vector3D solarArrayAxis,
  398.                                       final AbsoluteDate referenceDate,
  399.                                       final Vector3D referenceNormal,
  400.                                       final double rotationRate,
  401.                                       final double dragCoeff,
  402.                                       final double absorptionCoeff,
  403.                                       final double reflectionCoeff) {
  404.         this(facets, sun, solarArrayArea, solarArrayAxis, referenceDate, referenceNormal, rotationRate,
  405.              dragCoeff, 0.0, false,
  406.              absorptionCoeff, reflectionCoeff);
  407.     }

  408.     /** Build a spacecraft model with linear rotation of solar array.
  409.      * <p>
  410.      * The spacecraft body is described by an array of surface vectors. Each facet of
  411.      * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
  412.      * and whose norm is the surface area in m².
  413.      * </p>
  414.      * <p>
  415.      * Solar arrays orientation will be a regular rotation from the
  416.      * reference orientation at reference date and using a constant
  417.      * rotation rate.
  418.      * </p>
  419.      * @param facets body facets (only the facets with strictly positive area will be stored)
  420.      * @param sun sun model
  421.      * @param solarArrayArea area of the solar array (m²)
  422.      * @param solarArrayAxis solar array rotation axis in satellite frame
  423.      * @param referenceDate reference date for the solar array rotation
  424.      * @param referenceNormal direction of the solar array normal at reference date
  425.      * in spacecraft frame
  426.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  427.      * @param dragCoeff drag coefficient (used only for drag)
  428.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  429.      * that will experience specular reflection when hitting spacecraft instead
  430.      * of experiencing diffuse reflection, hence producing lift)
  431.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  432.      * (used only for radiation pressure)
  433.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  434.      * (used only for radiation pressure)
  435.      * @since 9.0
  436.      */
  437.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  438.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  439.                                       final Vector3D solarArrayAxis,
  440.                                       final AbsoluteDate referenceDate,
  441.                                       final Vector3D referenceNormal,
  442.                                       final double rotationRate,
  443.                                       final double dragCoeff, final double liftRatio,
  444.                                       final double absorptionCoeff,
  445.                                       final double reflectionCoeff) {
  446.         this(facets, sun, solarArrayArea, solarArrayAxis, referenceDate, referenceNormal, rotationRate,
  447.              dragCoeff, liftRatio, true,
  448.              absorptionCoeff, reflectionCoeff);
  449.     }

  450.     /** Build a spacecraft model with linear rotation of solar array.
  451.      * <p>
  452.      * The spacecraft body is described by an array of surface vectors. Each facet of
  453.      * the body is described by a vector normal to the facet (pointing outward of the spacecraft)
  454.      * and whose norm is the surface area in m².
  455.      * </p>
  456.      * <p>
  457.      * Solar arrays orientation will be a regular rotation from the
  458.      * reference orientation at reference date and using a constant
  459.      * rotation rate.
  460.      * </p>
  461.      * @param facets body facets (only the facets with strictly positive area will be stored)
  462.      * @param sun sun model
  463.      * @param solarArrayArea area of the solar array (m²)
  464.      * @param solarArrayAxis solar array rotation axis in satellite frame
  465.      * @param referenceDate reference date for the solar array rotation
  466.      * @param referenceNormal direction of the solar array normal at reference date
  467.      * in spacecraft frame
  468.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  469.      * @param dragCoeff drag coefficient (used only for drag)
  470.      * @param liftRatio lift ratio (proportion between 0 and 1 of atmosphere modecules
  471.      * that will experience specular reflection when hitting spacecraft instead
  472.      * of experiencing diffuse reflection, hence producing lift)
  473.      * @param useLift if true, lift should be used, otherwise it is completely ignored
  474.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  475.      * (used only for radiation pressure)
  476.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  477.      * (used only for radiation pressure)
  478.      * @since 9.0
  479.      */
  480.     private BoxAndSolarArraySpacecraft(final Facet[] facets,
  481.                                        final PVCoordinatesProvider sun, final double solarArrayArea,
  482.                                        final Vector3D solarArrayAxis,
  483.                                        final AbsoluteDate referenceDate,
  484.                                        final Vector3D referenceNormal,
  485.                                        final double rotationRate,
  486.                                        final double dragCoeff, final double liftRatio, final boolean useLift,
  487.                                        final double absorptionCoeff,
  488.                                        final double reflectionCoeff) {

  489.         // drag
  490.         dragParameterDriver = buildDragParameterDriver(dragCoeff);
  491.         liftParameterDriver = useLift ? buildLiftParameterDriver(liftRatio) : null;

  492.         // radiation pressure
  493.         absorptionParameterDriver = buildAbsorptionParameterDriver(absorptionCoeff);
  494.         reflectionParameterDriver = buildReflectionParameterDriver(reflectionCoeff);

  495.         factory = new DSFactory(1, 1);


  496.         this.facets = filter(facets.clone());

  497.         this.sun            = sun;
  498.         this.solarArrayArea = solarArrayArea;
  499.         this.referenceDate  = referenceDate;
  500.         this.rotationRate   = rotationRate;

  501.         this.saZ = solarArrayAxis.normalize();
  502.         this.saY = Vector3D.crossProduct(saZ, referenceNormal).normalize();
  503.         this.saX = Vector3D.crossProduct(saY, saZ);

  504.     }

  505.     /** Build the parameter driver for drag coefficient.
  506.      * @param coeff drag coefficient
  507.      * @return parameter driver for drag coefficient
  508.      * @since 9.0
  509.      */
  510.     private ParameterDriver buildDragParameterDriver(final double coeff) {
  511.         try {
  512.             return new ParameterDriver(DragSensitive.DRAG_COEFFICIENT,
  513.                                        coeff, SCALE, 0.0, Double.POSITIVE_INFINITY);
  514.         } catch (OrekitException oe) {
  515.             // this should never happen
  516.             throw new OrekitInternalError(oe);
  517.         }
  518.     }

  519.     /** Build the parameter driver for lift coefficient.
  520.      * @param coeff lift coefficient
  521.      * @return parameter driver for lift coefficient
  522.      * @since 9.0
  523.      */
  524.     private ParameterDriver buildLiftParameterDriver(final double coeff) {
  525.         try {
  526.             return new ParameterDriver(DragSensitive.LIFT_RATIO, coeff, SCALE, 0.0, 1.0);
  527.         } catch (OrekitException oe) {
  528.             // this should never happen
  529.             throw new OrekitInternalError(oe);
  530.         }
  531.     }

  532.     /** Build the parameter driver for absorption coefficient.
  533.      * @param coeff absorption coefficient
  534.      * @return parameter driver for absorption coefficient
  535.      * @since 9.0
  536.      */
  537.     private ParameterDriver buildAbsorptionParameterDriver(final double coeff) {
  538.         try {
  539.             return new ParameterDriver(RadiationSensitive.ABSORPTION_COEFFICIENT, coeff, SCALE, 0.0, 1.0);
  540.         } catch (OrekitException oe) {
  541.             // this should never happen
  542.             throw new OrekitInternalError(oe);
  543.         }
  544.     }

  545.     /** Build the parameter driver for reflection coefficient.
  546.      * @param coeff absorption coefficient
  547.      * @return parameter driver for reflection coefficient
  548.      * @since 9.0
  549.      */
  550.     private ParameterDriver buildReflectionParameterDriver(final double coeff) {
  551.         try {
  552.             return new ParameterDriver(RadiationSensitive.REFLECTION_COEFFICIENT, coeff, SCALE, 0.0, 1.0);
  553.         } catch (OrekitException oe) {
  554.             // this should never happen
  555.             throw new OrekitInternalError(oe);
  556.         }
  557.     }

  558.     /** {@inheritDoc} */
  559.     @Override
  560.     public ParameterDriver[] getDragParametersDrivers() {
  561.         return liftParameterDriver == null ?
  562.                new ParameterDriver[] {
  563.                    dragParameterDriver
  564.                } : new ParameterDriver[] {
  565.                    dragParameterDriver, liftParameterDriver
  566.                };
  567.     }

  568.     /** {@inheritDoc} */
  569.     @Override
  570.     public ParameterDriver[] getRadiationParametersDrivers() {
  571.         return new ParameterDriver[] {
  572.             absorptionParameterDriver, reflectionParameterDriver
  573.         };
  574.     }

  575.     /** Get solar array normal in spacecraft frame.
  576.      * @param date current date
  577.      * @param frame inertial reference frame for state (both orbit and attitude)
  578.      * @param position position of spacecraft in reference frame
  579.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  580.      * @return solar array normal in spacecraft frame
  581.      * @exception OrekitException if sun direction cannot be computed in best lighting
  582.      * configuration
  583.      */
  584.     public synchronized Vector3D getNormal(final AbsoluteDate date, final Frame frame,
  585.                                            final Vector3D position, final Rotation rotation)
  586.         throws OrekitException {

  587.         if (referenceDate != null) {
  588.             // use a simple rotation at fixed rate
  589.             final double alpha = rotationRate * date.durationFrom(referenceDate);
  590.             return new Vector3D(FastMath.cos(alpha), saX, FastMath.sin(alpha), saY);
  591.         }

  592.         // compute orientation for best lighting
  593.         final Vector3D sunInert = sun.getPVCoordinates(date, frame).getPosition().subtract(position).normalize();
  594.         final Vector3D sunSpacecraft = rotation.applyTo(sunInert);
  595.         final double d = Vector3D.dotProduct(sunSpacecraft, saZ);
  596.         final double f = 1 - d * d;
  597.         if (f < Precision.EPSILON) {
  598.             // extremely rare case: the sun is along solar array rotation axis
  599.             // (there will not be much output power ...)
  600.             // we set up an arbitrary normal
  601.             return saZ.orthogonal();
  602.         }

  603.         final double s = 1.0 / FastMath.sqrt(f);
  604.         return new Vector3D(s, sunSpacecraft, -s * d, saZ);

  605.     }

  606.     /** Get solar array normal in spacecraft frame.
  607.      * @param date current date
  608.      * @param frame inertial reference frame for state (both orbit and attitude)
  609.      * @param position position of spacecraft in reference frame
  610.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  611.      * @return solar array normal in spacecraft frame
  612.      * @param <T> type of the field elements
  613.      * @exception OrekitException if sun direction cannot be computed in best lighting
  614.      * configuration
  615.      */
  616.     public synchronized <T extends RealFieldElement<T>> FieldVector3D<T> getNormal(final FieldAbsoluteDate<T> date,
  617.                                                                                    final Frame frame,
  618.                                                                                    final FieldVector3D<T> position,
  619.                                                                                    final FieldRotation<T> rotation)
  620.         throws OrekitException {

  621.         if (referenceDate != null) {
  622.             // use a simple rotation at fixed rate
  623.             final T alpha = date.durationFrom(referenceDate).multiply(rotationRate);
  624.             return new FieldVector3D<>(alpha.cos(), saX, alpha.sin(), saY);
  625.         }

  626.         // compute orientation for best lighting
  627.         final FieldVector3D<T> sunInert = position.subtract(sun.getPVCoordinates(date.toAbsoluteDate(), frame).getPosition()).negate().normalize();
  628.         final FieldVector3D<T> sunSpacecraft = rotation.applyTo(sunInert);
  629.         final T d = FieldVector3D.dotProduct(sunSpacecraft, saZ);
  630.         final T f = d.multiply(d).subtract(1).negate();
  631.         if (f.getReal() < Precision.EPSILON) {
  632.             // extremely rare case: the sun is along solar array rotation axis
  633.             // (there will not be much output power ...)
  634.             // we set up an arbitrary normal
  635.             return new FieldVector3D<>(f.getField(), saZ.orthogonal());
  636.         }

  637.         final T s = f.sqrt().reciprocal();
  638.         return new FieldVector3D<>(s, sunSpacecraft,
  639.                                    s.multiply(d).negate(), new FieldVector3D<>(date.getField(), saZ));

  640.     }

  641.     /** Get solar array normal in spacecraft frame.
  642.      * @param date current date
  643.      * @param frame inertial reference frame for state (both orbit and attitude)
  644.      * @param position position of spacecraft in reference frame
  645.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  646.      * @return solar array normal in spacecraft frame
  647.      * @exception OrekitException if sun direction cannot be computed in best lighting
  648.      * configuration
  649.      */
  650.     public synchronized FieldVector3D<DerivativeStructure> getNormal(final AbsoluteDate date, final Frame frame,
  651.                                                                      final FieldVector3D<DerivativeStructure> position,
  652.                                                                      final FieldRotation<DerivativeStructure> rotation)
  653.         throws OrekitException {

  654.         final DerivativeStructure zero = position.getX().getField().getZero();

  655.         if (referenceDate != null) {
  656.             // use a simple rotation at fixed rate
  657.             final DerivativeStructure alpha = zero.add(rotationRate * date.durationFrom(referenceDate));
  658.             return new FieldVector3D<>(alpha.cos(), saX, alpha.sin(), saY);
  659.         }

  660.         // compute orientation for best lighting
  661.         final FieldVector3D<DerivativeStructure> sunInert =
  662.                 position.subtract(sun.getPVCoordinates(date, frame).getPosition()).negate().normalize();
  663.         final FieldVector3D<DerivativeStructure> sunSpacecraft = rotation.applyTo(sunInert);
  664.         final DerivativeStructure d = FieldVector3D.dotProduct(sunSpacecraft, saZ);
  665.         final DerivativeStructure f = d.multiply(d).subtract(1).negate();
  666.         if (f.getValue() < Precision.EPSILON) {
  667.             // extremely rare case: the sun is along solar array rotation axis
  668.             // (there will not be much output power ...)
  669.             // we set up an arbitrary normal
  670.             return new FieldVector3D<>(position.getX().getField(), saZ.orthogonal());
  671.         }

  672.         final DerivativeStructure s = f.sqrt().reciprocal();
  673.         return new FieldVector3D<>(s, sunSpacecraft,
  674.                                    s.multiply(d).negate(), new FieldVector3D<>(zero.getField(), saZ));

  675.     }


  676.     /** {@inheritDoc} */
  677.     @Override
  678.     public Vector3D dragAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
  679.                                      final Rotation rotation, final double mass,
  680.                                      final double density, final Vector3D relativeVelocity,
  681.                                      final double[] parameters)
  682.         throws OrekitException {

  683.         final double dragCoeff = parameters[0];
  684.         final double liftRatio = liftParameterDriver == null ? 0.0 : parameters[1];

  685.         // relative velocity in spacecraft frame
  686.         final double   vNorm2 = relativeVelocity.getNormSq();
  687.         final double   vNorm  = FastMath.sqrt(vNorm2);
  688.         final Vector3D vDir   = rotation.applyTo(relativeVelocity.scalarMultiply(1.0 / vNorm));
  689.         final double   coeff  = density * dragCoeff * vNorm2 / (2.0 * mass);
  690.         final double   oMr    = 1 - liftRatio;

  691.         // solar array contribution
  692.         final Vector3D frontNormal = getNormal(date, frame, position, rotation);
  693.         final double   s           = coeff * solarArrayArea * Vector3D.dotProduct(frontNormal, vDir);
  694.         Vector3D acceleration = new Vector3D(oMr * FastMath.abs(s), vDir,
  695.                                              liftRatio * s * 2,     frontNormal);

  696.         // body facets contribution
  697.         for (final Facet facet : facets) {
  698.             final double dot = Vector3D.dotProduct(facet.getNormal(), vDir);
  699.             if (dot < 0) {
  700.                 // the facet intercepts the incoming flux
  701.                 final double f = coeff * facet.getArea() * dot;
  702.                 acceleration = new Vector3D(1,                     acceleration,
  703.                                             oMr * FastMath.abs(f), vDir,
  704.                                             liftRatio * f * 2,     facet.getNormal());
  705.             }
  706.         }

  707.         // convert back to inertial frame
  708.         return rotation.applyInverseTo(acceleration);

  709.     }

  710.     /** {@inheritDoc} */
  711.     @Override
  712.     public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame,
  713.                                                                final Vector3D position, final Rotation rotation,
  714.                                                                final double mass, final  double density,
  715.                                                                final Vector3D relativeVelocity,
  716.                                                                final double[] parameters,
  717.                                                                final String paramName)
  718.         throws OrekitException {

  719.         final DerivativeStructure dragCoeffDS;
  720.         final DerivativeStructure liftRatioDS;
  721.         final DerivativeStructure oMrDS;
  722.         final Field<DerivativeStructure> field = factory.getDerivativeField();
  723.         if (dragParameterDriver.getName().equals(paramName)) {
  724.             final double liftRatio = liftParameterDriver == null ? 0.0 : parameters[1];
  725.             dragCoeffDS = factory.variable(0, parameters[0]);
  726.             liftRatioDS = factory.constant(liftRatio);
  727.             oMrDS       = factory.constant(1 - liftRatio);
  728.         } else if (liftParameterDriver != null && liftParameterDriver.getName().equals(paramName)) {
  729.             dragCoeffDS = factory.constant(parameters[0]);
  730.             liftRatioDS = factory.variable(0, parameters[1]);
  731.             oMrDS       = liftRatioDS.negate().add(1);
  732.         } else {
  733.             if (liftParameterDriver == null) {
  734.                 throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName,
  735.                                           dragParameterDriver.getName());
  736.             } else {
  737.                 throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName,
  738.                                           dragParameterDriver.getName(), liftParameterDriver.getName());
  739.             }
  740.         }

  741.         // relative velocity in spacecraft frame
  742.         final double                             vNorm2 = relativeVelocity.getNormSq();
  743.         final double                             vNorm  = FastMath.sqrt(vNorm2);
  744.         final FieldVector3D<DerivativeStructure> vDir   = new FieldVector3D<>(field,
  745.                                                                               rotation.applyTo(relativeVelocity.scalarMultiply(1.0 / vNorm)));
  746.         final DerivativeStructure                coeff  = dragCoeffDS.multiply(0.5 * density * vNorm2 / mass);

  747.         // solar array facet contribution
  748.         final FieldVector3D<DerivativeStructure> frontNormal = new FieldVector3D<>(field,
  749.                                                                                    getNormal(date, frame, position, rotation));
  750.         final DerivativeStructure                s           = coeff.
  751.                                                                multiply(solarArrayArea).
  752.                                                                multiply(FieldVector3D.dotProduct(frontNormal, vDir));
  753.         FieldVector3D<DerivativeStructure> acceleration = new FieldVector3D<>(s.abs().multiply(oMrDS), vDir,
  754.                                                                               s.multiply(liftRatioDS).multiply(2), frontNormal);

  755.         // body facets contribution
  756.         for (final Facet facet : facets) {
  757.             final DerivativeStructure dot = FieldVector3D.dotProduct(facet.getNormal(), vDir);
  758.             if (dot.getValue() < 0) {
  759.                 // the facet intercepts the incoming flux
  760.                 final DerivativeStructure f = coeff.multiply(facet.getArea()).multiply(dot);
  761.                 acceleration = new FieldVector3D<>(field.getOne(),          acceleration,
  762.                                                    f.abs().multiply(oMrDS), vDir,
  763.                                                    f.multiply(liftRatioDS).multiply(2), new FieldVector3D<>(field, facet.getNormal()));
  764.             }
  765.         }

  766.         // convert back to inertial frame
  767.         return new FieldRotation<>(field, rotation).applyInverseTo(acceleration);

  768.     }

  769.     /** {@inheritDoc} */
  770.     @Override
  771.     public Vector3D radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
  772.                                                   final Rotation rotation, final double mass, final Vector3D flux,
  773.                                                   final double[] parameters)
  774.         throws OrekitException {

  775.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  776.             // null illumination (we are probably in umbra)
  777.             return Vector3D.ZERO;
  778.         }

  779.         // radiation flux in spacecraft frame
  780.         final Vector3D fluxSat = rotation.applyTo(flux);

  781.         // solar array contribution
  782.         Vector3D normal = getNormal(date, frame, position, rotation);
  783.         double dot = Vector3D.dotProduct(normal, fluxSat);
  784.         if (dot > 0) {
  785.             // the solar array is illuminated backward,
  786.             // fix signs to compute contribution correctly
  787.             dot   = -dot;
  788.             normal = normal.negate();
  789.         }
  790.         Vector3D force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot, parameters);

  791.         // body facets contribution
  792.         for (final Facet bodyFacet : facets) {
  793.             normal = bodyFacet.getNormal();
  794.             dot = Vector3D.dotProduct(normal, fluxSat);
  795.             if (dot < 0) {
  796.                 // the facet intercepts the incoming flux
  797.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot, parameters));
  798.             }
  799.         }

  800.         // convert to inertial frame
  801.         return rotation.applyInverseTo(new Vector3D(1.0 / mass, force));

  802.     }

  803.     /** {@inheritDoc} */
  804.     @Override
  805.     public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame,
  806.                                                                             final Vector3D position, final Rotation rotation,
  807.                                                                             final double mass, final Vector3D flux,
  808.                                                                             final double[] parameters,
  809.                                                                             final String paramName)
  810.         throws OrekitException {

  811.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  812.             // null illumination (we are probably in umbra)
  813.             return FieldVector3D.getZero(factory.getDerivativeField());
  814.         }

  815.         final DerivativeStructure absorptionCoeffDS;
  816.         final DerivativeStructure specularReflectionCoeffDS;
  817.         if (ABSORPTION_COEFFICIENT.equals(paramName)) {
  818.             absorptionCoeffDS         = factory.variable(0, parameters[0]);
  819.             specularReflectionCoeffDS = factory.constant(parameters[1]);
  820.         } else if (REFLECTION_COEFFICIENT.equals(paramName)) {
  821.             absorptionCoeffDS         = factory.constant(parameters[0]);
  822.             specularReflectionCoeffDS = factory.variable(0, parameters[1]);
  823.         } else {
  824.             throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName,
  825.                                       ABSORPTION_COEFFICIENT + ", " + REFLECTION_COEFFICIENT);
  826.         }
  827.         final DerivativeStructure diffuseReflectionCoeffDS =
  828.                 absorptionCoeffDS.add(specularReflectionCoeffDS).subtract(1).negate();

  829.         // radiation flux in spacecraft frame
  830.         final Vector3D fluxSat = rotation.applyTo(flux);

  831.         // solar array contribution
  832.         Vector3D normal = getNormal(date, frame, position, rotation);
  833.         double dot = Vector3D.dotProduct(normal, fluxSat);
  834.         if (dot > 0) {
  835.             // the solar array is illuminated backward,
  836.             // fix signs to compute contribution correctly
  837.             dot   = -dot;
  838.             normal = normal.negate();
  839.         }
  840.         FieldVector3D<DerivativeStructure> force =
  841.                 facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot,
  842.                                            specularReflectionCoeffDS, diffuseReflectionCoeffDS);

  843.         // body facets contribution
  844.         for (final Facet bodyFacet : facets) {
  845.             normal = bodyFacet.getNormal();
  846.             dot = Vector3D.dotProduct(normal, fluxSat);
  847.             if (dot < 0) {
  848.                 // the facet intercepts the incoming flux
  849.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot,
  850.                                                              specularReflectionCoeffDS, diffuseReflectionCoeffDS));
  851.             }
  852.         }

  853.         // convert to inertial
  854.         return FieldRotation.applyInverseTo(rotation, new FieldVector3D<>(1.0 / mass, force));

  855.     }

  856.     /** {@inheritDoc} */
  857.     @Override
  858.     public <T extends RealFieldElement<T>> FieldVector3D<T>
  859.         dragAcceleration(final FieldAbsoluteDate<T> date, final Frame frame,
  860.                          final FieldVector3D<T> position, final FieldRotation<T> rotation,
  861.                          final T mass, final  T density, final FieldVector3D<T> relativeVelocity,
  862.                          final T[] parameters)
  863.         throws OrekitException {

  864.         final T dragCoeff = parameters[0];
  865.         final T liftRatio = liftParameterDriver == null ? dragCoeff.getField().getZero() : parameters[1];

  866.         // relative velocity in spacecraft frame
  867.         final T                vNorm2 = relativeVelocity.getNormSq();
  868.         final T                vNorm  = vNorm2.sqrt();
  869.         final FieldVector3D<T> vDir   = rotation.applyTo(relativeVelocity.scalarMultiply(vNorm.reciprocal()));
  870.         final T                coeff  = density.multiply(0.5).multiply(dragCoeff).multiply(vNorm2).divide(mass);
  871.         final T                oMr    = liftRatio.negate().add(1);

  872.         // solar array facet contribution
  873.         final FieldVector3D<T> frontNormal = getNormal(date, frame, position, rotation);
  874.         final T                s           = coeff.
  875.                                              multiply(solarArrayArea).
  876.                                              multiply(FieldVector3D.dotProduct(frontNormal, vDir));
  877.         FieldVector3D<T> acceleration = new FieldVector3D<>(s.abs().multiply(oMr), vDir,
  878.                                                             s.multiply(liftRatio).multiply(2), frontNormal);

  879.         // body facets contribution
  880.         final Field<T> field = coeff.getField();
  881.         for (final Facet facet : facets) {
  882.             final T dot = FieldVector3D.dotProduct(facet.getNormal(), vDir);
  883.             if (dot.getReal() < 0) {
  884.                 // the facet intercepts the incoming flux
  885.                 final T f = coeff.multiply(facet.getArea()).multiply(dot);
  886.                 acceleration = new FieldVector3D<>(field.getOne(),        acceleration,
  887.                                                    f.abs().multiply(oMr), vDir,
  888.                                                    f.multiply(liftRatio).multiply(2), new FieldVector3D<>(field, facet.getNormal()));
  889.             }
  890.         }

  891.         // convert back to inertial frame
  892.         return rotation.applyInverseTo(acceleration);

  893.     }

  894.     /** {@inheritDoc} */
  895.     @Override
  896.     public <T extends RealFieldElement<T>> FieldVector3D<T>
  897.         radiationPressureAcceleration(final FieldAbsoluteDate<T> date, final Frame frame,
  898.                                       final FieldVector3D<T> position,
  899.                                       final FieldRotation<T> rotation, final T mass,
  900.                                       final FieldVector3D<T> flux,
  901.                                       final T[] parameters)
  902.         throws OrekitException {

  903.         if (flux.getNormSq().getReal() < Precision.SAFE_MIN) {
  904.             // null illumination (we are probably in umbra)
  905.             return FieldVector3D.getZero(date.getField());
  906.         }

  907.         // radiation flux in spacecraft frame
  908.         final FieldVector3D<T> fluxSat = rotation.applyTo(flux);

  909.         // solar array contribution
  910.         FieldVector3D<T> normal = getNormal(date, frame, position, rotation);
  911.         T dot = FieldVector3D.dotProduct(normal, fluxSat);
  912.         if (dot.getReal() > 0) {
  913.             // the solar array is illuminated backward,
  914.             // fix signs to compute contribution correctly
  915.             dot    = dot.negate();
  916.             normal = normal.negate();
  917.         }
  918.         FieldVector3D<T> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot, parameters);

  919.         // body facets contribution
  920.         for (final Facet bodyFacet : facets) {
  921.             normal = new FieldVector3D<>(date.getField(), bodyFacet.getNormal());
  922.             dot = FieldVector3D.dotProduct(fluxSat, normal);
  923.             if (dot.getReal() < 0) {
  924.                 // the facet intercepts the incoming flux
  925.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot, parameters));
  926.             }
  927.         }

  928.         // convert to inertial frame
  929.         return rotation.applyInverseTo(new FieldVector3D<>(mass.reciprocal(), force));

  930.     }

  931.     /** Compute contribution of one facet to force.
  932.      * <p>This method implements equation 8-44 from David A. Vallado's
  933.      * Fundamentals of Astrodynamics and Applications, third edition,
  934.      * 2007, Microcosm Press.</p>
  935.      * @param normal facet normal
  936.      * @param area facet area
  937.      * @param fluxSat radiation pressure flux in spacecraft frame
  938.      * @param dot dot product of facet and fluxSat (must be negative)
  939.      * @param parameters values of the force model parameters
  940.      * @return contribution of the facet to force in spacecraft frame
  941.      */
  942.     private Vector3D facetRadiationAcceleration(final Vector3D normal, final double area, final Vector3D fluxSat,
  943.                                                 final double dot, final double[] parameters) {

  944.         final double absorptionCoeff         = parameters[0];
  945.         final double specularReflectionCoeff = parameters[1];
  946.         final double diffuseReflectionCoeff  = 1 - (absorptionCoeff + specularReflectionCoeff);

  947.         final double psr  = fluxSat.getNorm();

  948.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  949.         // cos (phi) = -dot / (psr * area)
  950.         // n         = facet / area
  951.         // s         = -fluxSat / psr
  952.         final double cN = 2 * area * dot * (diffuseReflectionCoeff / 3 - specularReflectionCoeff * dot / psr);
  953.         final double cS = (area * dot / psr) * (specularReflectionCoeff - 1);
  954.         return new Vector3D(cN, normal, cS, fluxSat);

  955.     }

  956.     /** Compute contribution of one facet to force.
  957.      * <p>This method implements equation 8-44 from David A. Vallado's
  958.      * Fundamentals of Astrodynamics and Applications, third edition,
  959.      * 2007, Microcosm Press.</p>
  960.      * @param normal facet normal
  961.      * @param area facet area
  962.      * @param fluxSat radiation pressure flux in spacecraft frame
  963.      * @param dot dot product of facet and fluxSat (must be negative)
  964.      * @param parameters values of the force model parameters
  965.      * @param <T> type of the field elements
  966.      * @return contribution of the facet to force in spacecraft frame
  967.      */
  968.     private <T extends RealFieldElement<T>> FieldVector3D<T>
  969.         facetRadiationAcceleration(final FieldVector3D<T> normal, final double area, final FieldVector3D<T> fluxSat,
  970.                                    final T dot, final T[] parameters) {

  971.         final T absorptionCoeff         = parameters[0];
  972.         final T specularReflectionCoeff = parameters[1];
  973.         final T diffuseReflectionCoeff  = absorptionCoeff.add(specularReflectionCoeff).negate().add(1);

  974.         final T psr  = fluxSat.getNorm();

  975.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  976.         // cos (phi) = -dot / (psr * area)
  977.         // n         = facet / area
  978.         // s         = -fluxSat / psr
  979.         final T cN = dot.multiply(-2 * area).multiply(dot.multiply(specularReflectionCoeff).divide(psr).subtract(diffuseReflectionCoeff.divide(3)));
  980.         final T cS = dot.multiply(area).multiply(specularReflectionCoeff.subtract(1)).divide(psr);
  981.         return new FieldVector3D<>(cN, normal, cS, fluxSat);

  982.     }

  983.     /** Compute contribution of one facet to force.
  984.      * <p>This method implements equation 8-44 from David A. Vallado's
  985.      * Fundamentals of Astrodynamics and Applications, third edition,
  986.      * 2007, Microcosm Press.</p>
  987.      * @param normal facet normal
  988.      * @param area facet area
  989.      * @param fluxSat radiation pressure flux in spacecraft frame
  990.      * @param dot dot product of facet and fluxSat (must be negative)
  991.      * @param specularReflectionCoeffDS specular reflection coefficient
  992.      * @param diffuseReflectionCoeffDS diffuse reflection coefficient
  993.      * @return contribution of the facet to force in spacecraft frame
  994.      */
  995.     private FieldVector3D<DerivativeStructure> facetRadiationAcceleration(final Vector3D normal, final double area,
  996.                                                                           final Vector3D fluxSat, final double dot,
  997.                                                                           final DerivativeStructure specularReflectionCoeffDS,
  998.                                                                           final DerivativeStructure diffuseReflectionCoeffDS) {
  999.         final double psr  = fluxSat.getNorm();

  1000.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  1001.         // cos (phi) = -dot / (psr * area)
  1002.         // n         = facet / area
  1003.         // s         = -fluxSat / psr
  1004.         final DerivativeStructure cN =
  1005.                 diffuseReflectionCoeffDS.divide(3).subtract(specularReflectionCoeffDS.multiply(dot / psr)).multiply(2 * area * dot);
  1006.         final DerivativeStructure cS = specularReflectionCoeffDS.subtract(1).multiply(area * dot / psr);

  1007.         return new FieldVector3D<>(cN, normal, cS, fluxSat);

  1008.     }

  1009.     /** Class representing a single facet of a convex spacecraft body.
  1010.      * <p>Instance of this class are guaranteed to be immutable.</p>
  1011.      * @author Luc Maisonobe
  1012.      */
  1013.     public static class Facet {

  1014.         /** Unit Normal vector, pointing outward. */
  1015.         private final Vector3D normal;

  1016.         /** Area in m². */
  1017.         private final double area;

  1018.         /** Simple constructor.
  1019.          * @param normal vector normal to the facet, pointing outward (will be normalized)
  1020.          * @param area facet area in m²
  1021.          */
  1022.         public Facet(final Vector3D normal, final double area) {
  1023.             this.normal = normal.normalize();
  1024.             this.area   = area;
  1025.         }

  1026.         /** Get unit normal vector.
  1027.          * @return unit normal vector
  1028.          */
  1029.         public Vector3D getNormal() {
  1030.             return normal;
  1031.         }

  1032.         /** Get facet area.
  1033.          * @return facet area
  1034.          */
  1035.         public double getArea() {
  1036.             return area;
  1037.         }

  1038.     }

  1039.     /** Build the surface vectors for body facets of a simple parallelepipedic box.
  1040.      * @param xLength length of the body along its X axis (m)
  1041.      * @param yLength length of the body along its Y axis (m)
  1042.      * @param zLength length of the body along its Z axis (m)
  1043.      * @return surface vectors array
  1044.      */
  1045.     private static Facet[] simpleBoxFacets(final double xLength, final double yLength, final double zLength) {
  1046.         return new Facet[] {
  1047.             new Facet(Vector3D.MINUS_I, yLength * zLength),
  1048.             new Facet(Vector3D.PLUS_I,  yLength * zLength),
  1049.             new Facet(Vector3D.MINUS_J, xLength * zLength),
  1050.             new Facet(Vector3D.PLUS_J,  xLength * zLength),
  1051.             new Facet(Vector3D.MINUS_K, xLength * yLength),
  1052.             new Facet(Vector3D.PLUS_K,  xLength * yLength)
  1053.         };
  1054.     }

  1055.     /** Filter out zero area facets.
  1056.      * @param facets original facets (may include zero area facets)
  1057.      * @return filtered array
  1058.      */
  1059.     private static List<Facet> filter(final Facet[] facets) {
  1060.         final List<Facet> filtered = new ArrayList<Facet>(facets.length);
  1061.         for (Facet facet : facets) {
  1062.             if (facet.getArea() > 0) {
  1063.                 filtered.add(facet);
  1064.             }
  1065.         }
  1066.         return filtered;
  1067.     }

  1068. }