BoxAndSolarArraySpacecraft.java

  1. /* Copyright 2002-2016 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.analysis.differentiation.DerivativeStructure;
  21. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  22. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  23. import org.hipparchus.geometry.euclidean.threed.Rotation;
  24. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.Precision;
  27. import org.orekit.errors.OrekitException;
  28. import org.orekit.errors.OrekitInternalError;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.forces.drag.DragSensitive;
  31. import org.orekit.forces.radiation.RadiationSensitive;
  32. import org.orekit.frames.Frame;
  33. import org.orekit.time.AbsoluteDate;
  34. import org.orekit.utils.PVCoordinatesProvider;
  35. import org.orekit.utils.ParameterDriver;
  36. import org.orekit.utils.ParameterObserver;

  37. /** experimental class representing the features of a classical satellite
  38.  * with a convex body shape and rotating flat solar arrays.
  39.  * <p>
  40.  * As of 5.0, this class is still considered experimental, so use it with care.
  41.  * </p>
  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.  * This model does not take cast shadow between body and solar array into account.
  57.  * </p>
  58.  *
  59.  * @author Luc Maisonobe
  60.  * @author Pascal Parraud
  61.  */
  62. public class BoxAndSolarArraySpacecraft implements RadiationSensitive, DragSensitive {

  63.     /** Parameters scaling factor.
  64.      * <p>
  65.      * We use a power of 2 to avoid numeric noise introduction
  66.      * in the multiplications/divisions sequences.
  67.      * </p>
  68.      */
  69.     private final double SCALE = FastMath.scalb(1.0, -3);

  70.     /** Drivers for drag coefficient parameter. */
  71.     private final ParameterDriver[] dragParametersDrivers;

  72.     /** Drivers for radiation pressure coefficient parameter. */
  73.     private final ParameterDriver[] radiationParametersDrivers;

  74.     /** Surface vectors for body facets. */
  75.     private final List<Facet> facets;

  76.     /** Solar array area (m²). */
  77.     private final double solarArrayArea;

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

  80.     /** Rotation rate for linear rotation angle. */
  81.     private final double rotationRate;

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

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

  86.     /** Solar array rotation axis in spacecraft frame. */
  87.     private final Vector3D saZ;

  88.     /** Drag coefficient. */
  89.     private double dragCoeff;

  90.     /** Absorption coefficient. */
  91.     private double absorptionCoeff;

  92.     /** Specular reflection coefficient. */
  93.     private double specularReflectionCoeff;

  94.     /** Diffuse reflection coefficient. */
  95.     private double diffuseReflectionCoeff;

  96.     /** Sun model. */
  97.     private final PVCoordinatesProvider sun;

  98.     /** Build a spacecraft model with best lighting of solar array.
  99.      * <p>
  100.      * Solar arrays orientation will be such that at each time the Sun direction
  101.      * will always be in the solar array meridian plane defined by solar array
  102.      * rotation axis and solar array normal vector.
  103.      * </p>
  104.      * @param xLength length of the body along its X axis (m)
  105.      * @param yLength length of the body along its Y axis (m)
  106.      * @param zLength length of the body along its Z axis (m)
  107.      * @param sun sun model
  108.      * @param solarArrayArea area of the solar array (m²)
  109.      * @param solarArrayAxis solar array rotation axis in satellite frame
  110.      * @param dragCoeff drag coefficient (used only for drag)
  111.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  112.      * (used only for radiation pressure)
  113.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  114.      * (used only for radiation pressure)
  115.      */
  116.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  117.                                       final double zLength,
  118.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  119.                                       final Vector3D solarArrayAxis,
  120.                                       final double dragCoeff,
  121.                                       final double absorptionCoeff,
  122.                                       final double reflectionCoeff) {
  123.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  124.              dragCoeff, absorptionCoeff, reflectionCoeff);
  125.     }

  126.     /** Build a spacecraft model with best lighting of solar array.
  127.      * <p>
  128.      * The spacecraft body is described by an array of surface vectors. Each facet of
  129.      * the body is describe by a vector normal to the facet (pointing outward of the spacecraft)
  130.      * and whose norm is the surface area in m².
  131.      * </p>
  132.      * <p>
  133.      * Solar arrays orientation will be such that at each time the Sun direction
  134.      * will always be in the solar array meridian plane defined by solar array
  135.      * rotation axis and solar array normal vector.
  136.      * </p>
  137.      * @param facets body facets (only the facets with strictly positive area will be stored)
  138.      * @param sun sun model
  139.      * @param solarArrayArea area of the solar array (m²)
  140.      * @param solarArrayAxis solar array rotation axis in satellite frame
  141.      * @param dragCoeff drag coefficient (used only for drag)
  142.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  143.      * (used only for radiation pressure)
  144.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  145.      * (used only for radiation pressure)
  146.      */
  147.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  148.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  149.                                       final Vector3D solarArrayAxis,
  150.                                       final double dragCoeff,
  151.                                       final double absorptionCoeff,
  152.                                       final double reflectionCoeff) {

  153.         this.dragParametersDrivers      = new ParameterDriver[1];
  154.         this.radiationParametersDrivers = new ParameterDriver[2];
  155.         try {
  156.             dragParametersDrivers[0] = new ParameterDriver(DragSensitive.DRAG_COEFFICIENT,
  157.                                                            dragCoeff, SCALE, 0.0, Double.POSITIVE_INFINITY);
  158.             dragParametersDrivers[0].addObserver(new ParameterObserver() {
  159.                 /** {@inheridDoc} */
  160.                 @Override
  161.                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
  162.                     BoxAndSolarArraySpacecraft.this.dragCoeff = driver.getValue();
  163.                 }
  164.             });
  165.             radiationParametersDrivers[0] = new ParameterDriver(RadiationSensitive.ABSORPTION_COEFFICIENT,
  166.                                                                 absorptionCoeff, SCALE, 0.0, 1.0);
  167.             radiationParametersDrivers[0].addObserver(new ParameterObserver() {
  168.                 /** {@inheridDoc} */
  169.                 @Override
  170.                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
  171.                     BoxAndSolarArraySpacecraft.this.absorptionCoeff = driver.getValue();
  172.                     BoxAndSolarArraySpacecraft.this.diffuseReflectionCoeff =
  173.                                     1 - (driver.getValue() + BoxAndSolarArraySpacecraft.this.specularReflectionCoeff);
  174.                 }
  175.             });
  176.             radiationParametersDrivers[1] = new ParameterDriver(RadiationSensitive.REFLECTION_COEFFICIENT,
  177.                                                                 reflectionCoeff, SCALE, 0.0, 1.0);
  178.             radiationParametersDrivers[1].addObserver(new ParameterObserver() {
  179.                 /** {@inheridDoc} */
  180.                 @Override
  181.                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
  182.                     BoxAndSolarArraySpacecraft.this.specularReflectionCoeff = driver.getValue();
  183.                     BoxAndSolarArraySpacecraft.this.diffuseReflectionCoeff  =
  184.                                     1 - (BoxAndSolarArraySpacecraft.this.absorptionCoeff + driver.getValue());
  185.                 }
  186.             });
  187.         } catch (OrekitException oe) {
  188.             // this should never happen
  189.             throw new OrekitInternalError(oe);
  190.         }

  191.         this.facets = filter(facets);

  192.         this.sun            = sun;
  193.         this.solarArrayArea = solarArrayArea;
  194.         this.referenceDate  = null;
  195.         this.rotationRate   = 0;

  196.         this.saZ = solarArrayAxis.normalize();
  197.         this.saY = null;
  198.         this.saX = null;

  199.         this.dragCoeff               = dragCoeff;
  200.         this.absorptionCoeff         = absorptionCoeff;
  201.         this.specularReflectionCoeff = reflectionCoeff;
  202.         this.diffuseReflectionCoeff  = 1 - (absorptionCoeff + reflectionCoeff);
  203.     }

  204.     /** Build a spacecraft model with linear rotation of solar array.
  205.      * <p>
  206.      * Solar arrays orientation will be a regular rotation from the
  207.      * reference orientation at reference date and using a constant
  208.      * rotation rate.
  209.      * </p>
  210.      * @param xLength length of the body along its X axis (m)
  211.      * @param yLength length of the body along its Y axis (m)
  212.      * @param zLength length of the body along its Z axis (m)
  213.      * @param sun sun model
  214.      * @param solarArrayArea area of the solar array (m²)
  215.      * @param solarArrayAxis solar array rotation axis in satellite frame
  216.      * @param referenceDate reference date for the solar array rotation
  217.      * @param referenceNormal direction of the solar array normal at reference date
  218.      * in spacecraft frame
  219.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  220.      * @param dragCoeff drag coefficient (used only for drag)
  221.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  222.      * (used only for radiation pressure)
  223.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  224.      * (used only for radiation pressure)
  225.      */
  226.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  227.                                       final double zLength,
  228.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  229.                                       final Vector3D solarArrayAxis,
  230.                                       final AbsoluteDate referenceDate,
  231.                                       final Vector3D referenceNormal,
  232.                                       final double rotationRate,
  233.                                       final double dragCoeff,
  234.                                       final double absorptionCoeff,
  235.                                       final double reflectionCoeff) {
  236.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  237.              referenceDate, referenceNormal, rotationRate,
  238.              dragCoeff, absorptionCoeff, reflectionCoeff);
  239.     }

  240.     /** Build a spacecraft model with linear rotation of solar array.
  241.      * <p>
  242.      * The spacecraft body is described by an array of surface vectors. Each facet of
  243.      * the body is describe by a vector normal to the facet (pointing outward of the spacecraft)
  244.      * and whose norm is the surface area in m².
  245.      * </p>
  246.      * <p>
  247.      * Solar arrays orientation will be a regular rotation from the
  248.      * reference orientation at reference date and using a constant
  249.      * rotation rate.
  250.      * </p>
  251.      * @param facets body facets (only the facets with strictly positive area will be stored)
  252.      * @param sun sun model
  253.      * @param solarArrayArea area of the solar array (m²)
  254.      * @param solarArrayAxis solar array rotation axis in satellite frame
  255.      * @param referenceDate reference date for the solar array rotation
  256.      * @param referenceNormal direction of the solar array normal at reference date
  257.      * in spacecraft frame
  258.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  259.      * @param dragCoeff drag coefficient (used only for drag)
  260.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  261.      * (used only for radiation pressure)
  262.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  263.      * (used only for radiation pressure)
  264.      */
  265.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  266.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  267.                                       final Vector3D solarArrayAxis,
  268.                                       final AbsoluteDate referenceDate,
  269.                                       final Vector3D referenceNormal,
  270.                                       final double rotationRate,
  271.                                       final double dragCoeff,
  272.                                       final double absorptionCoeff,
  273.                                       final double reflectionCoeff) {

  274.         this.dragParametersDrivers      = new ParameterDriver[1];
  275.         this.radiationParametersDrivers = new ParameterDriver[2];
  276.         try {
  277.             dragParametersDrivers[0] = new ParameterDriver(DragSensitive.DRAG_COEFFICIENT,
  278.                                                            dragCoeff, SCALE,
  279.                                                            Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
  280.             dragParametersDrivers[0].addObserver(new ParameterObserver() {
  281.                 /** {@inheridDoc} */
  282.                 @Override
  283.                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
  284.                     BoxAndSolarArraySpacecraft.this.dragCoeff = driver.getValue();
  285.                 }
  286.             });
  287.             radiationParametersDrivers[0] = new ParameterDriver(RadiationSensitive.ABSORPTION_COEFFICIENT,
  288.                                                                 absorptionCoeff, SCALE, 0.0, 1.0);
  289.             radiationParametersDrivers[0].addObserver(new ParameterObserver() {
  290.                 /** {@inheridDoc} */
  291.                 @Override
  292.                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
  293.                     BoxAndSolarArraySpacecraft.this.absorptionCoeff = driver.getValue();
  294.                     BoxAndSolarArraySpacecraft.this.diffuseReflectionCoeff =
  295.                                     1 - (driver.getValue() + BoxAndSolarArraySpacecraft.this.specularReflectionCoeff);
  296.                 }
  297.             });
  298.             radiationParametersDrivers[1] = new ParameterDriver(RadiationSensitive.REFLECTION_COEFFICIENT,
  299.                                                                 reflectionCoeff, SCALE, 0.0, 1.0);
  300.             radiationParametersDrivers[1].addObserver(new ParameterObserver() {
  301.                 /** {@inheridDoc} */
  302.                 @Override
  303.                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
  304.                     BoxAndSolarArraySpacecraft.this.specularReflectionCoeff = driver.getValue();
  305.                     BoxAndSolarArraySpacecraft.this.diffuseReflectionCoeff  =
  306.                                     1 - (BoxAndSolarArraySpacecraft.this.absorptionCoeff + driver.getValue());
  307.                 }
  308.             });
  309.         } catch (OrekitException oe) {
  310.             // this should never happen
  311.             throw new OrekitInternalError(oe);
  312.         }

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

  314.         this.sun            = sun;
  315.         this.solarArrayArea = solarArrayArea;
  316.         this.referenceDate  = referenceDate;
  317.         this.rotationRate   = rotationRate;

  318.         this.saZ = solarArrayAxis.normalize();
  319.         this.saY = Vector3D.crossProduct(saZ, referenceNormal).normalize();
  320.         this.saX = Vector3D.crossProduct(saY, saZ);

  321.         this.dragCoeff               = dragCoeff;
  322.         this.absorptionCoeff         = absorptionCoeff;
  323.         this.specularReflectionCoeff = reflectionCoeff;
  324.         this.diffuseReflectionCoeff  = 1 - (absorptionCoeff + reflectionCoeff);

  325.     }

  326.     /** {@inheritDoc} */
  327.     @Override
  328.     public ParameterDriver[] getDragParametersDrivers() {
  329.         return dragParametersDrivers.clone();
  330.     }

  331.     /** {@inheritDoc} */
  332.     @Override
  333.     public ParameterDriver[] getRadiationParametersDrivers() {
  334.         return radiationParametersDrivers.clone();
  335.     }

  336.     /** Get solar array normal in spacecraft frame.
  337.      * @param date current date
  338.      * @param frame inertial reference frame for state (both orbit and attitude)
  339.      * @param position position of spacecraft in reference frame
  340.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  341.      * @return solar array normal in spacecraft frame
  342.      * @exception OrekitException if sun direction cannot be computed in best lighting
  343.      * configuration
  344.      */
  345.     public synchronized Vector3D getNormal(final AbsoluteDate date, final Frame frame,
  346.                                            final Vector3D position, final Rotation rotation)
  347.         throws OrekitException {

  348.         if (referenceDate != null) {
  349.             // use a simple rotation at fixed rate
  350.             final double alpha = rotationRate * date.durationFrom(referenceDate);
  351.             return new Vector3D(FastMath.cos(alpha), saX, FastMath.sin(alpha), saY);
  352.         }

  353.         // compute orientation for best lighting
  354.         final Vector3D sunInert = sun.getPVCoordinates(date, frame).getPosition().subtract(position).normalize();
  355.         final Vector3D sunSpacecraft = rotation.applyTo(sunInert);
  356.         final double d = Vector3D.dotProduct(sunSpacecraft, saZ);
  357.         final double f = 1 - d * d;
  358.         if (f < Precision.EPSILON) {
  359.             // extremely rare case: the sun is along solar array rotation axis
  360.             // (there will not be much output power ...)
  361.             // we set up an arbitrary normal
  362.             return saZ.orthogonal();
  363.         }

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

  366.     }


  367.     /** Get solar array normal in spacecraft frame.
  368.      * @param date current date
  369.      * @param frame inertial reference frame for state (both orbit and attitude)
  370.      * @param position position of spacecraft in reference frame
  371.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  372.      * @return solar array normal in spacecraft frame
  373.      * @exception OrekitException if sun direction cannot be computed in best lighting
  374.      * configuration
  375.      */
  376.     public synchronized FieldVector3D<DerivativeStructure> getNormal(final AbsoluteDate date, final Frame frame,
  377.                                                                      final FieldVector3D<DerivativeStructure> position,
  378.                                                                      final FieldRotation<DerivativeStructure> rotation)
  379.         throws OrekitException {

  380.         final DerivativeStructure one = position.getX().getField().getOne();

  381.         if (referenceDate != null) {
  382.             // use a simple rotation at fixed rate
  383.             final DerivativeStructure alpha = one.multiply(rotationRate * date.durationFrom(referenceDate));
  384.             return new FieldVector3D<DerivativeStructure>(alpha.cos(), saX, alpha.sin(), saY);
  385.         }

  386.         // compute orientation for best lighting
  387.         final FieldVector3D<DerivativeStructure> sunInert =
  388.                 position.subtract(sun.getPVCoordinates(date, frame).getPosition()).negate().normalize();
  389.         final FieldVector3D<DerivativeStructure> sunSpacecraft = rotation.applyTo(sunInert);
  390.         final DerivativeStructure d = FieldVector3D.dotProduct(sunSpacecraft, saZ);
  391.         final DerivativeStructure f = d.multiply(d).subtract(1).negate();
  392.         if (f.getValue() < Precision.EPSILON) {
  393.             // extremely rare case: the sun is along solar array rotation axis
  394.             // (there will not be much output power ...)
  395.             // we set up an arbitrary normal
  396.             return new FieldVector3D<DerivativeStructure>(one, saZ.orthogonal());
  397.         }

  398.         final DerivativeStructure s = f.sqrt().reciprocal();
  399.         return new FieldVector3D<DerivativeStructure>(s, sunSpacecraft).subtract(new FieldVector3D<DerivativeStructure>(s.multiply(d), saZ));

  400.     }


  401.     /** {@inheritDoc} */
  402.     public Vector3D dragAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
  403.                                      final Rotation rotation, final double mass,
  404.                                      final double density, final Vector3D relativeVelocity)
  405.         throws OrekitException {

  406.         // relative velocity in spacecraft frame
  407.         final Vector3D v = rotation.applyTo(relativeVelocity);

  408.         // solar array contribution
  409.         final Vector3D solarArrayFacet = new Vector3D(solarArrayArea, getNormal(date, frame, position, rotation));
  410.         double sv = FastMath.abs(Vector3D.dotProduct(solarArrayFacet, v));

  411.         // body facets contribution
  412.         for (final Facet facet : facets) {
  413.             final double dot = Vector3D.dotProduct(facet.getNormal(), v);
  414.             if (dot < 0) {
  415.                 // the facet intercepts the incoming flux
  416.                 sv -= facet.getArea() * dot;
  417.             }
  418.         }

  419.         return new Vector3D(sv * density * dragCoeff / (2.0 * mass), relativeVelocity);

  420.     }

  421.     /** {@inheritDoc} */
  422.     public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame,
  423.                                                                final FieldVector3D<DerivativeStructure> position,
  424.                                                                final FieldRotation<DerivativeStructure> rotation,
  425.                                                                final DerivativeStructure mass,
  426.                                                                final DerivativeStructure density,
  427.                                                                final FieldVector3D<DerivativeStructure> relativeVelocity)
  428.         throws OrekitException {

  429.         // relative velocity in spacecraft frame
  430.         final FieldVector3D<DerivativeStructure> v = rotation.applyTo(relativeVelocity);

  431.         // solar array contribution
  432.         final FieldVector3D<DerivativeStructure> solarArrayFacet =
  433.                 new FieldVector3D<DerivativeStructure>(solarArrayArea, getNormal(date, frame, position, rotation));
  434.         DerivativeStructure sv = FieldVector3D.dotProduct(v, solarArrayFacet).abs();

  435.         // body facets contribution
  436.         for (final Facet facet : facets) {
  437.             final DerivativeStructure dot = FieldVector3D.dotProduct(v, facet.getNormal());
  438.             if (dot.getValue() < 0) {
  439.                 // the facet intercepts the incoming flux
  440.                 sv = sv.subtract(dot.multiply(facet.getArea()));
  441.             }
  442.         }

  443.         return new FieldVector3D<DerivativeStructure>(sv.multiply(density.multiply(dragCoeff / 2.0)).divide(mass),
  444.                                                       relativeVelocity);

  445.     }

  446.     /** {@inheritDoc} */
  447.     public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame,
  448.                                                                final Vector3D position, final Rotation rotation,
  449.                                                                final double mass, final  double density,
  450.                                                                final Vector3D relativeVelocity,
  451.                                                                final String paramName)
  452.         throws OrekitException {

  453.         if (!DRAG_COEFFICIENT.equals(paramName)) {
  454.             throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName, DRAG_COEFFICIENT);
  455.         }

  456.         final DerivativeStructure dragCoeffDS = new DerivativeStructure(1, 1, 0, dragCoeff);

  457.         // relative velocity in spacecraft frame
  458.         final Vector3D v = rotation.applyTo(relativeVelocity);

  459.         // solar array contribution
  460.         final Vector3D solarArrayFacet = new Vector3D(solarArrayArea, getNormal(date, frame, position, rotation));
  461.         double sv = FastMath.abs(Vector3D.dotProduct(solarArrayFacet, v));

  462.         // body facets contribution
  463.         for (final Facet facet : facets) {
  464.             final double dot = Vector3D.dotProduct(facet.getNormal(), v);
  465.             if (dot < 0) {
  466.                 // the facet intercepts the incoming flux
  467.                 sv -= facet.getArea() * dot;
  468.             }
  469.         }

  470.         return new FieldVector3D<DerivativeStructure>(dragCoeffDS.multiply(sv * density / (2.0 * mass)),
  471.                                                       relativeVelocity);

  472.     }

  473.     /** {@inheritDoc} */
  474.     public Vector3D radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
  475.                                                   final Rotation rotation, final double mass, final Vector3D flux)
  476.         throws OrekitException {

  477.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  478.             // null illumination (we are probably in umbra)
  479.             return Vector3D.ZERO;
  480.         }

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

  483.         // solar array contribution
  484.         Vector3D normal = getNormal(date, frame, position, rotation);
  485.         double dot = Vector3D.dotProduct(normal, fluxSat);
  486.         if (dot > 0) {
  487.             // the solar array is illuminated backward,
  488.             // fix signs to compute contribution correctly
  489.             dot   = -dot;
  490.             normal = normal.negate();
  491.         }
  492.         Vector3D force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot);

  493.         // body facets contribution
  494.         for (final Facet bodyFacet : facets) {
  495.             normal = bodyFacet.getNormal();
  496.             dot = Vector3D.dotProduct(normal, fluxSat);
  497.             if (dot < 0) {
  498.                 // the facet intercepts the incoming flux
  499.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot));
  500.             }
  501.         }

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

  504.     }

  505.     /** {@inheritDoc} */
  506.     public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame,
  507.                                                                             final FieldVector3D<DerivativeStructure> position,
  508.                                                                             final FieldRotation<DerivativeStructure> rotation,
  509.                                                                             final DerivativeStructure mass,
  510.                                                                             final FieldVector3D<DerivativeStructure> flux)
  511.         throws OrekitException {

  512.         if (flux.getNormSq().getValue() < Precision.SAFE_MIN) {
  513.             // null illumination (we are probably in umbra)
  514.             return new FieldVector3D<DerivativeStructure>(0.0, flux);
  515.         }

  516.         // radiation flux in spacecraft frame
  517.         final FieldVector3D<DerivativeStructure> fluxSat = rotation.applyTo(flux);

  518.         // solar array contribution
  519.         FieldVector3D<DerivativeStructure> normal = getNormal(date, frame, position, rotation);
  520.         DerivativeStructure dot = FieldVector3D.dotProduct(normal, fluxSat);
  521.         if (dot.getValue() > 0) {
  522.             // the solar array is illuminated backward,
  523.             // fix signs to compute contribution correctly
  524.             dot    = dot.negate();
  525.             normal = normal.negate();
  526.         }
  527.         FieldVector3D<DerivativeStructure> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot);

  528.         // body facets contribution
  529.         for (final Facet bodyFacet : facets) {
  530.             normal = new FieldVector3D<DerivativeStructure>(mass.getField().getOne(), bodyFacet.getNormal());
  531.             dot = FieldVector3D.dotProduct(normal, fluxSat);
  532.             if (dot.getValue() < 0) {
  533.                 // the facet intercepts the incoming flux
  534.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot));
  535.             }
  536.         }

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

  539.     }

  540.     /** {@inheritDoc} */
  541.     public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame,
  542.                                                                             final Vector3D position, final Rotation rotation,
  543.                                                                             final double mass, final Vector3D flux,
  544.                                                                             final String paramName)
  545.         throws OrekitException {

  546.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  547.             // null illumination (we are probably in umbra)
  548.             final DerivativeStructure zero = new DerivativeStructure(1, 1, 0.0);
  549.             return new FieldVector3D<DerivativeStructure>(zero, zero, zero);
  550.         }

  551.         final DerivativeStructure absorptionCoeffDS;
  552.         final DerivativeStructure specularReflectionCoeffDS;
  553.         if (ABSORPTION_COEFFICIENT.equals(paramName)) {
  554.             absorptionCoeffDS         = new DerivativeStructure(1, 1, 0, absorptionCoeff);
  555.             specularReflectionCoeffDS = new DerivativeStructure(1, 1,    specularReflectionCoeff);
  556.         } else if (REFLECTION_COEFFICIENT.equals(paramName)) {
  557.             absorptionCoeffDS         = new DerivativeStructure(1, 1,    absorptionCoeff);
  558.             specularReflectionCoeffDS = new DerivativeStructure(1, 1, 0, specularReflectionCoeff);
  559.         } else {
  560.             throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName,
  561.                                       ABSORPTION_COEFFICIENT + ", " + REFLECTION_COEFFICIENT);
  562.         }
  563.         final DerivativeStructure diffuseReflectionCoeffDS =
  564.                 absorptionCoeffDS.add(specularReflectionCoeffDS).subtract(1).negate();


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

  567.         // solar array contribution
  568.         Vector3D normal = getNormal(date, frame, position, rotation);
  569.         double dot = Vector3D.dotProduct(normal, fluxSat);
  570.         if (dot > 0) {
  571.             // the solar array is illuminated backward,
  572.             // fix signs to compute contribution correctly
  573.             dot   = -dot;
  574.             normal = normal.negate();
  575.         }
  576.         FieldVector3D<DerivativeStructure> force =
  577.                 facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot,
  578.                                            specularReflectionCoeffDS, diffuseReflectionCoeffDS);

  579.         // body facets contribution
  580.         for (final Facet bodyFacet : facets) {
  581.             normal = bodyFacet.getNormal();
  582.             dot = Vector3D.dotProduct(normal, fluxSat);
  583.             if (dot < 0) {
  584.                 // the facet intercepts the incoming flux
  585.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot,
  586.                                                              specularReflectionCoeffDS, diffuseReflectionCoeffDS));
  587.             }
  588.         }

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

  591.     }

  592.     /** Compute contribution of one facet to force.
  593.      * <p>This method implements equation 8-44 from David A. Vallado's
  594.      * Fundamentals of Astrodynamics and Applications, third edition,
  595.      * 2007, Microcosm Press.</p>
  596.      * @param normal facet normal
  597.      * @param area facet area
  598.      * @param fluxSat radiation pressure flux in spacecraft frame
  599.      * @param dot dot product of facet and fluxSat (must be negative)
  600.      * @return contribution of the facet to force in spacecraft frame
  601.      */
  602.     private Vector3D facetRadiationAcceleration(final Vector3D normal, final double area, final Vector3D fluxSat,
  603.                                                 final double dot) {
  604.         final double psr  = fluxSat.getNorm();

  605.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  606.         // cos (phi) = -dot / (psr * area)
  607.         // n         = facet / area
  608.         // s         = -fluxSat / psr
  609.         final double cN = 2 * area * dot * (diffuseReflectionCoeff / 3 - specularReflectionCoeff * dot / psr);
  610.         final double cS = (area * dot / psr) * (specularReflectionCoeff - 1);
  611.         return new Vector3D(cN, normal, cS, fluxSat);

  612.     }

  613.     /** Compute contribution of one facet to force.
  614.      * <p>This method implements equation 8-44 from David A. Vallado's
  615.      * Fundamentals of Astrodynamics and Applications, third edition,
  616.      * 2007, Microcosm Press.</p>
  617.      * @param normal facet normal
  618.      * @param area facet area
  619.      * @param fluxSat radiation pressure flux in spacecraft frame
  620.      * @param dot dot product of facet and fluxSat (must be negative)
  621.      * @return contribution of the facet to force in spacecraft frame
  622.      */
  623.     private FieldVector3D<DerivativeStructure> facetRadiationAcceleration(final FieldVector3D<DerivativeStructure> normal,
  624.                                                                           final double area,
  625.                                                                           final FieldVector3D<DerivativeStructure> fluxSat,
  626.                                                                           final DerivativeStructure dot) {
  627.         final DerivativeStructure psr  = fluxSat.getNorm();

  628.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  629.         // cos (phi) = -dot / (psr * area)
  630.         // n         = facet / area
  631.         // s         = -fluxSat / psr
  632.         final DerivativeStructure cN = dot.multiply(-2 * area).multiply(dot.divide(psr).multiply(specularReflectionCoeff).subtract(diffuseReflectionCoeff / 3));
  633.         final DerivativeStructure cS = dot.divide(psr).multiply(area * (specularReflectionCoeff - 1));
  634.         return new FieldVector3D<DerivativeStructure>(cN, normal, cS, fluxSat);

  635.     }

  636.     /** Compute contribution of one facet to force.
  637.      * <p>This method implements equation 8-44 from David A. Vallado's
  638.      * Fundamentals of Astrodynamics and Applications, third edition,
  639.      * 2007, Microcosm Press.</p>
  640.      * @param normal facet normal
  641.      * @param area facet area
  642.      * @param fluxSat radiation pressure flux in spacecraft frame
  643.      * @param dot dot product of facet and fluxSat (must be negative)
  644.      * @param specularReflectionCoeffDS specular reflection coefficient
  645.      * @param diffuseReflectionCoeffDS diffuse reflection coefficient
  646.      * @return contribution of the facet to force in spacecraft frame
  647.      */
  648.     private FieldVector3D<DerivativeStructure> facetRadiationAcceleration(final Vector3D normal, final double area,
  649.                                                                           final Vector3D fluxSat, final double dot,
  650.                                                                           final DerivativeStructure specularReflectionCoeffDS,
  651.                                                                           final DerivativeStructure diffuseReflectionCoeffDS) {
  652.         final double psr  = fluxSat.getNorm();

  653.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  654.         // cos (phi) = -dot / (psr * area)
  655.         // n         = facet / area
  656.         // s         = -fluxSat / psr
  657.         final DerivativeStructure cN =
  658.                 diffuseReflectionCoeffDS.divide(3).subtract(specularReflectionCoeffDS.multiply(dot / psr)).multiply(2 * area * dot);
  659.         final DerivativeStructure cS = specularReflectionCoeffDS.subtract(1).multiply(area * dot / psr);

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

  661.     }

  662.     /** Class representing a single facet of a convex spacecraft body.
  663.      * <p>Instance of this class are guaranteed to be immutable.</p>
  664.      * @author Luc Maisonobe
  665.      */
  666.     public static class Facet {

  667.         /** Unit Normal vector, pointing outward. */
  668.         private final Vector3D normal;

  669.         /** Area in m². */
  670.         private final double area;

  671.         /** Simple constructor.
  672.          * @param normal vector normal to the facet, pointing outward (will be normalized)
  673.          * @param area facet area in m²
  674.          */
  675.         public Facet(final Vector3D normal, final double area) {
  676.             this.normal = normal.normalize();
  677.             this.area   = area;
  678.         }

  679.         /** Get unit normal vector.
  680.          * @return unit normal vector
  681.          */
  682.         public Vector3D getNormal() {
  683.             return normal;
  684.         }

  685.         /** Get facet area.
  686.          * @return facet area
  687.          */
  688.         public double getArea() {
  689.             return area;
  690.         }

  691.     }

  692.     /** Build the surface vectors for body facets of a simple parallelepipedic box.
  693.      * @param xLength length of the body along its X axis (m)
  694.      * @param yLength length of the body along its Y axis (m)
  695.      * @param zLength length of the body along its Z axis (m)
  696.      * @return surface vectors array
  697.      */
  698.     private static Facet[] simpleBoxFacets(final double xLength, final double yLength, final double zLength) {
  699.         return new Facet[] {
  700.             new Facet(Vector3D.MINUS_I, yLength * zLength),
  701.             new Facet(Vector3D.PLUS_I,  yLength * zLength),
  702.             new Facet(Vector3D.MINUS_J, xLength * zLength),
  703.             new Facet(Vector3D.PLUS_J,  xLength * zLength),
  704.             new Facet(Vector3D.MINUS_K, xLength * yLength),
  705.             new Facet(Vector3D.PLUS_K,  xLength * yLength)
  706.         };
  707.     }

  708.     /** Filter out zero area facets.
  709.      * @param facets original facets (may include zero area facets)
  710.      * @return filtered array
  711.      */
  712.     private static List<Facet> filter(final Facet[] facets) {
  713.         final List<Facet> filtered = new ArrayList<Facet>(facets.length);
  714.         for (Facet facet : facets) {
  715.             if (facet.getArea() > 0) {
  716.                 filtered.add(facet);
  717.             }
  718.         }
  719.         return filtered;
  720.     }

  721. }