BoxAndSolarArraySpacecraft.java

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

  34. /** experimental class representing the features of a classical satellite
  35.  * with a convex body shape and rotating flat solar arrays.
  36.  * <p>
  37.  * As of 5.0, this class is still considered experimental, so use it with care.
  38.  * </p>
  39.  * <p>
  40.  * The body can be either a simple parallelepipedic box aligned with
  41.  * spacecraft axes or a set of facets defined by their area and normal vector.
  42.  * This should handle accurately most spacecraft shapes.
  43.  * </p>
  44.  * <p>
  45.  * The solar array rotation with respect to satellite body can be either
  46.  * the best lightning orientation (i.e. Sun exactly in solar array meridian
  47.  * plane defined by solar array rotation axis and solar array normal vector)
  48.  * or a rotation evolving linearly according to a start position and an
  49.  * angular rate (which can be set to 0 for non-rotating panels, which may
  50.  * occur in special modes or during contingencies).
  51.  * </p>
  52.  * <p>
  53.  * This model does not take cast shadow between body and solar array into account.
  54.  * </p>
  55.  * <p>
  56.  * Instances of this class are guaranteed to be immutable.
  57.  * </p>
  58.  *
  59.  * @see SphericalSpacecraft
  60.  * @author Luc Maisonobe
  61.  * @author Pascal Parraud
  62.  */
  63. public class BoxAndSolarArraySpacecraft implements RadiationSensitive, DragSensitive {

  64.     /** Surface vectors for body facets. */
  65.     private final List<Facet> facets;

  66.     /** Solar array area (m<sup>2</sup>). */
  67.     private final double solarArrayArea;

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

  70.     /** Rotation rate for linear rotation angle. */
  71.     private final double rotationRate;

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

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

  76.     /** Solar array rotation axis in spacecraft frame. */
  77.     private final Vector3D saZ;

  78.     /** Drag coefficient. */
  79.     private double dragCoeff;

  80.     /** Absorption coefficient. */
  81.     private double absorptionCoeff;

  82.     /** Specular reflection coefficient. */
  83.     private double specularReflectionCoeff;

  84.     /** Diffuse reflection coefficient. */
  85.     private double diffuseReflectionCoeff;

  86.     /** Sun model. */
  87.     private final PVCoordinatesProvider sun;

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

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

  143.         this.facets = filter(facets);

  144.         this.sun            = sun;
  145.         this.solarArrayArea = solarArrayArea;
  146.         this.referenceDate  = null;
  147.         this.rotationRate   = 0;

  148.         this.saZ = solarArrayAxis.normalize();
  149.         this.saY = null;
  150.         this.saX = null;

  151.         this.dragCoeff               = dragCoeff;
  152.         this.absorptionCoeff         = absorptionCoeff;
  153.         this.specularReflectionCoeff = reflectionCoeff;
  154.         this.diffuseReflectionCoeff  = 1 - (absorptionCoeff + reflectionCoeff);
  155.     }

  156.     /** Build a spacecraft model with linear rotation of solar array.
  157.      * <p>
  158.      * Solar arrays orientation will be a regular rotation from the
  159.      * reference orientation at reference date and using a constant
  160.      * rotation rate.
  161.      * </p>
  162.      * @param xLength length of the body along its X axis (m)
  163.      * @param yLength length of the body along its Y axis (m)
  164.      * @param zLength length of the body along its Z axis (m)
  165.      * @param sun sun model
  166.      * @param solarArrayArea area of the solar array (m<sup>2</sup>)
  167.      * @param solarArrayAxis solar array rotation axis in satellite frame
  168.      * @param referenceDate reference date for the solar array rotation
  169.      * @param referenceNormal direction of the solar array normal at reference date
  170.      * in spacecraft frame
  171.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  172.      * @param dragCoeff drag coefficient (used only for drag)
  173.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  174.      * (used only for radiation pressure)
  175.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  176.      * (used only for radiation pressure)
  177.      */
  178.     public BoxAndSolarArraySpacecraft(final double xLength, final double yLength,
  179.                                       final double zLength,
  180.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  181.                                       final Vector3D solarArrayAxis,
  182.                                       final AbsoluteDate referenceDate,
  183.                                       final Vector3D referenceNormal,
  184.                                       final double rotationRate,
  185.                                       final double dragCoeff,
  186.                                       final double absorptionCoeff,
  187.                                       final double reflectionCoeff) {
  188.         this(simpleBoxFacets(xLength, yLength, zLength), sun, solarArrayArea, solarArrayAxis,
  189.              referenceDate, referenceNormal, rotationRate,
  190.              dragCoeff, absorptionCoeff, reflectionCoeff);
  191.     }

  192.     /** Build a spacecraft model with linear rotation of solar array.
  193.      * <p>
  194.      * The spacecraft body is described by an array of surface vectors. Each facet of
  195.      * the body is describe by a vector normal to the facet (pointing outward of the spacecraft)
  196.      * and whose norm is the surface area in m<sup>2</sup>.
  197.      * </p>
  198.      * <p>
  199.      * Solar arrays orientation will be a regular rotation from the
  200.      * reference orientation at reference date and using a constant
  201.      * rotation rate.
  202.      * </p>
  203.      * @param facets body facets (only the facets with strictly positive area will be stored)
  204.      * @param sun sun model
  205.      * @param solarArrayArea area of the solar array (m<sup>2</sup>)
  206.      * @param solarArrayAxis solar array rotation axis in satellite frame
  207.      * @param referenceDate reference date for the solar array rotation
  208.      * @param referenceNormal direction of the solar array normal at reference date
  209.      * in spacecraft frame
  210.      * @param rotationRate rotation rate of the solar array, may be 0 (rad/s)
  211.      * @param dragCoeff drag coefficient (used only for drag)
  212.      * @param absorptionCoeff absorption coefficient between 0.0 an 1.0
  213.      * (used only for radiation pressure)
  214.      * @param reflectionCoeff specular reflection coefficient between 0.0 an 1.0
  215.      * (used only for radiation pressure)
  216.      */
  217.     public BoxAndSolarArraySpacecraft(final Facet[] facets,
  218.                                       final PVCoordinatesProvider sun, final double solarArrayArea,
  219.                                       final Vector3D solarArrayAxis,
  220.                                       final AbsoluteDate referenceDate,
  221.                                       final Vector3D referenceNormal,
  222.                                       final double rotationRate,
  223.                                       final double dragCoeff,
  224.                                       final double absorptionCoeff,
  225.                                       final double reflectionCoeff) {

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

  227.         this.sun            = sun;
  228.         this.solarArrayArea = solarArrayArea;
  229.         this.referenceDate  = referenceDate;
  230.         this.rotationRate   = rotationRate;

  231.         this.saZ = solarArrayAxis.normalize();
  232.         this.saY = Vector3D.crossProduct(saZ, referenceNormal).normalize();
  233.         this.saX = Vector3D.crossProduct(saY, saZ);

  234.         this.dragCoeff               = dragCoeff;
  235.         this.absorptionCoeff         = absorptionCoeff;
  236.         this.specularReflectionCoeff = reflectionCoeff;
  237.         this.diffuseReflectionCoeff  = 1 - (absorptionCoeff + reflectionCoeff);

  238.     }

  239.     /** Get solar array normal in spacecraft frame.
  240.      * @param date current date
  241.      * @param frame inertial reference frame for state (both orbit and attitude)
  242.      * @param position position of spacecraft in reference frame
  243.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  244.      * @return solar array normal in spacecraft frame
  245.      * @exception OrekitException if sun direction cannot be computed in best lightning
  246.      * configuration
  247.      */
  248.     public synchronized Vector3D getNormal(final AbsoluteDate date, final Frame frame,
  249.                                            final Vector3D position, final Rotation rotation)
  250.         throws OrekitException {

  251.         if (referenceDate != null) {
  252.             // use a simple rotation at fixed rate
  253.             final double alpha = rotationRate * date.durationFrom(referenceDate);
  254.             return new Vector3D(FastMath.cos(alpha), saX, FastMath.sin(alpha), saY);
  255.         }

  256.         // compute orientation for best lightning
  257.         final Vector3D sunInert = sun.getPVCoordinates(date, frame).getPosition().subtract(position).normalize();
  258.         final Vector3D sunSpacecraft = rotation.applyTo(sunInert);
  259.         final double d = Vector3D.dotProduct(sunSpacecraft, saZ);
  260.         final double f = 1 - d * d;
  261.         if (f < Precision.EPSILON) {
  262.             // extremely rare case: the sun is along solar array rotation axis
  263.             // (there will not be much output power ...)
  264.             // we set up an arbitrary normal
  265.             return saZ.orthogonal();
  266.         }

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

  269.     }


  270.     /** Get solar array normal in spacecraft frame.
  271.      * @param date current date
  272.      * @param frame inertial reference frame for state (both orbit and attitude)
  273.      * @param position position of spacecraft in reference frame
  274.      * @param rotation orientation (attitude) of the spacecraft with respect to reference frame
  275.      * @return solar array normal in spacecraft frame
  276.      * @exception OrekitException if sun direction cannot be computed in best lightning
  277.      * configuration
  278.      */
  279.     public synchronized FieldVector3D<DerivativeStructure> getNormal(final AbsoluteDate date, final Frame frame,
  280.                                                                      final FieldVector3D<DerivativeStructure> position,
  281.                                                                      final FieldRotation<DerivativeStructure> rotation)
  282.         throws OrekitException {

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

  284.         if (referenceDate != null) {
  285.             // use a simple rotation at fixed rate
  286.             final DerivativeStructure alpha = one.multiply(rotationRate * date.durationFrom(referenceDate));
  287.             return new FieldVector3D<DerivativeStructure>(alpha.cos(), saX, alpha.sin(), saY);
  288.         }

  289.         // compute orientation for best lightning
  290.         final FieldVector3D<DerivativeStructure> sunInert =
  291.                 position.subtract(sun.getPVCoordinates(date, frame).getPosition()).negate().normalize();
  292.         final FieldVector3D<DerivativeStructure> sunSpacecraft = rotation.applyTo(sunInert);
  293.         final DerivativeStructure d = FieldVector3D.dotProduct(sunSpacecraft, saZ);
  294.         final DerivativeStructure f = d.multiply(d).subtract(1).negate();
  295.         if (f.getValue() < Precision.EPSILON) {
  296.             // extremely rare case: the sun is along solar array rotation axis
  297.             // (there will not be much output power ...)
  298.             // we set up an arbitrary normal
  299.             return new FieldVector3D<DerivativeStructure>(one, saZ.orthogonal());
  300.         }

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

  303.     }


  304.     /** {@inheritDoc} */
  305.     public Vector3D dragAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
  306.                                      final Rotation rotation, final double mass,
  307.                                      final double density, final Vector3D relativeVelocity)
  308.         throws OrekitException {

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

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

  314.         // body facets contribution
  315.         for (final Facet facet : facets) {
  316.             final double dot = Vector3D.dotProduct(facet.getNormal(), v);
  317.             if (dot < 0) {
  318.                 // the facet intercepts the incoming flux
  319.                 sv -= facet.getArea() * dot;
  320.             }
  321.         }

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

  323.     }

  324.     /** {@inheritDoc} */
  325.     public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame,
  326.                                                                final FieldVector3D<DerivativeStructure> position,
  327.                                                                final FieldRotation<DerivativeStructure> rotation,
  328.                                                                final DerivativeStructure mass,
  329.                                                                final double density,
  330.                                                                final FieldVector3D<DerivativeStructure> relativeVelocity)
  331.         throws OrekitException {

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

  334.         // solar array contribution
  335.         final FieldVector3D<DerivativeStructure> solarArrayFacet =
  336.                 new FieldVector3D<DerivativeStructure>(solarArrayArea, getNormal(date, frame, position, rotation));
  337.         DerivativeStructure sv = FieldVector3D.dotProduct(v, solarArrayFacet).abs();

  338.         // body facets contribution
  339.         for (final Facet facet : facets) {
  340.             final DerivativeStructure dot = FieldVector3D.dotProduct(v, facet.getNormal());
  341.             if (dot.getValue() < 0) {
  342.                 // the facet intercepts the incoming flux
  343.                 sv = sv.subtract(dot.multiply(facet.getArea()));
  344.             }
  345.         }

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

  348.     }

  349.     /** {@inheritDoc} */
  350.     public FieldVector3D<DerivativeStructure> dragAcceleration(final AbsoluteDate date, final Frame frame,
  351.                                                                final Vector3D position, final Rotation rotation,
  352.                                                                final double mass, final  double density,
  353.                                                                final Vector3D relativeVelocity,
  354.                                                                final String paramName)
  355.         throws OrekitException {

  356.         if (!DRAG_COEFFICIENT.equals(paramName)) {
  357.             throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName, DRAG_COEFFICIENT);
  358.         }

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

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

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

  365.         // body facets contribution
  366.         for (final Facet facet : facets) {
  367.             final double dot = Vector3D.dotProduct(facet.getNormal(), v);
  368.             if (dot < 0) {
  369.                 // the facet intercepts the incoming flux
  370.                 sv -= facet.getArea() * dot;
  371.             }
  372.         }

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

  375.     }

  376.     /** {@inheritDoc} */
  377.     public Vector3D radiationPressureAcceleration(final AbsoluteDate date, final Frame frame, final Vector3D position,
  378.                                                   final Rotation rotation, final double mass, final Vector3D flux)
  379.         throws OrekitException {

  380.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  381.             // null illumination (we are probably in umbra)
  382.             return Vector3D.ZERO;
  383.         }

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

  386.         // solar array contribution
  387.         Vector3D normal = getNormal(date, frame, position, rotation);
  388.         double dot = Vector3D.dotProduct(normal, fluxSat);
  389.         if (dot > 0) {
  390.             // the solar array is illuminated backward,
  391.             // fix signs to compute contribution correctly
  392.             dot   = -dot;
  393.             normal = normal.negate();
  394.         }
  395.         Vector3D force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot);

  396.         // body facets contribution
  397.         for (final Facet bodyFacet : facets) {
  398.             normal = bodyFacet.getNormal();
  399.             dot = Vector3D.dotProduct(normal, fluxSat);
  400.             if (dot < 0) {
  401.                 // the facet intercepts the incoming flux
  402.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot));
  403.             }
  404.         }

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

  407.     }

  408.     /** {@inheritDoc} */
  409.     public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame,
  410.                                                                             final FieldVector3D<DerivativeStructure> position,
  411.                                                                             final FieldRotation<DerivativeStructure> rotation,
  412.                                                                             final DerivativeStructure mass,
  413.                                                                             final FieldVector3D<DerivativeStructure> flux)
  414.         throws OrekitException {

  415.         if (flux.getNormSq().getValue() < Precision.SAFE_MIN) {
  416.             // null illumination (we are probably in umbra)
  417.             return new FieldVector3D<DerivativeStructure>(0.0, flux);
  418.         }

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

  421.         // solar array contribution
  422.         FieldVector3D<DerivativeStructure> normal = getNormal(date, frame, position, rotation);
  423.         DerivativeStructure dot = FieldVector3D.dotProduct(normal, fluxSat);
  424.         if (dot.getValue() > 0) {
  425.             // the solar array is illuminated backward,
  426.             // fix signs to compute contribution correctly
  427.             dot    = dot.negate();
  428.             normal = normal.negate();
  429.         }
  430.         FieldVector3D<DerivativeStructure> force = facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot);

  431.         // body facets contribution
  432.         for (final Facet bodyFacet : facets) {
  433.             normal = new FieldVector3D<DerivativeStructure>(mass.getField().getOne(), bodyFacet.getNormal());
  434.             dot = FieldVector3D.dotProduct(normal, fluxSat);
  435.             if (dot.getValue() < 0) {
  436.                 // the facet intercepts the incoming flux
  437.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot));
  438.             }
  439.         }

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

  442.     }

  443.     /** {@inheritDoc} */
  444.     public FieldVector3D<DerivativeStructure> radiationPressureAcceleration(final AbsoluteDate date, final Frame frame,
  445.                                                                             final Vector3D position, final Rotation rotation,
  446.                                                                             final double mass, final Vector3D flux,
  447.                                                                             final String paramName)
  448.         throws OrekitException {

  449.         if (flux.getNormSq() < Precision.SAFE_MIN) {
  450.             // null illumination (we are probably in umbra)
  451.             final DerivativeStructure zero = new DerivativeStructure(1, 1, 0.0);
  452.             return new FieldVector3D<DerivativeStructure>(zero, zero, zero);
  453.         }

  454.         final DerivativeStructure absorptionCoeffDS;
  455.         final DerivativeStructure specularReflectionCoeffDS;
  456.         if (ABSORPTION_COEFFICIENT.equals(paramName)) {
  457.             absorptionCoeffDS         = new DerivativeStructure(1, 1, 0, absorptionCoeff);
  458.             specularReflectionCoeffDS = new DerivativeStructure(1, 1,    specularReflectionCoeff);
  459.         } else if (REFLECTION_COEFFICIENT.equals(paramName)) {
  460.             absorptionCoeffDS         = new DerivativeStructure(1, 1,    absorptionCoeff);
  461.             specularReflectionCoeffDS = new DerivativeStructure(1, 1, 0, specularReflectionCoeff);
  462.         } else {
  463.             throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, paramName,
  464.                                       ABSORPTION_COEFFICIENT + ", " + REFLECTION_COEFFICIENT);
  465.         }
  466.         final DerivativeStructure diffuseReflectionCoeffDS =
  467.                 absorptionCoeffDS.add(specularReflectionCoeffDS).subtract(1).negate();


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

  470.         // solar array contribution
  471.         Vector3D normal = getNormal(date, frame, position, rotation);
  472.         double dot = Vector3D.dotProduct(normal, fluxSat);
  473.         if (dot > 0) {
  474.             // the solar array is illuminated backward,
  475.             // fix signs to compute contribution correctly
  476.             dot   = -dot;
  477.             normal = normal.negate();
  478.         }
  479.         FieldVector3D<DerivativeStructure> force =
  480.                 facetRadiationAcceleration(normal, solarArrayArea, fluxSat, dot,
  481.                                            specularReflectionCoeffDS, diffuseReflectionCoeffDS);

  482.         // body facets contribution
  483.         for (final Facet bodyFacet : facets) {
  484.             normal = bodyFacet.getNormal();
  485.             dot = Vector3D.dotProduct(normal, fluxSat);
  486.             if (dot < 0) {
  487.                 // the facet intercepts the incoming flux
  488.                 force = force.add(facetRadiationAcceleration(normal, bodyFacet.getArea(), fluxSat, dot,
  489.                                                              specularReflectionCoeffDS, diffuseReflectionCoeffDS));
  490.             }
  491.         }

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

  494.     }

  495.     /** Compute contribution of one facet to force.
  496.      * <p>This method implements equation 8-44 from David A. Vallado's
  497.      * Fundamentals of Astrodynamics and Applications, third edition,
  498.      * 2007, Microcosm Press.</p>
  499.      * @param normal facet normal
  500.      * @param area facet area
  501.      * @param fluxSat radiation pressure flux in spacecraft frame
  502.      * @param dot dot product of facet and fluxSat (must be negative)
  503.      * @return contribution of the facet to force in spacecraft frame
  504.      */
  505.     private Vector3D facetRadiationAcceleration(final Vector3D normal, final double area, final Vector3D fluxSat,
  506.                                                 final double dot) {
  507.         final double psr  = fluxSat.getNorm();

  508.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  509.         // cos (phi) = -dot / (psr * area)
  510.         // n         = facet / area
  511.         // s         = -fluxSat / psr
  512.         final double cN = 2 * area * dot * (diffuseReflectionCoeff / 3 - specularReflectionCoeff * dot / psr);
  513.         final double cS = (area * dot / psr) * (specularReflectionCoeff - 1);
  514.         return new Vector3D(cN, normal, cS, fluxSat);

  515.     }

  516.     /** Compute contribution of one facet to force.
  517.      * <p>This method implements equation 8-44 from David A. Vallado's
  518.      * Fundamentals of Astrodynamics and Applications, third edition,
  519.      * 2007, Microcosm Press.</p>
  520.      * @param normal facet normal
  521.      * @param area facet area
  522.      * @param fluxSat radiation pressure flux in spacecraft frame
  523.      * @param dot dot product of facet and fluxSat (must be negative)
  524.      * @return contribution of the facet to force in spacecraft frame
  525.      */
  526.     private FieldVector3D<DerivativeStructure> facetRadiationAcceleration(final FieldVector3D<DerivativeStructure> normal,
  527.                                                                           final double area,
  528.                                                                           final FieldVector3D<DerivativeStructure> fluxSat,
  529.                                                                           final DerivativeStructure dot) {
  530.         final DerivativeStructure psr  = fluxSat.getNorm();

  531.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  532.         // cos (phi) = -dot / (psr * area)
  533.         // n         = facet / area
  534.         // s         = -fluxSat / psr
  535.         final DerivativeStructure cN = dot.multiply(-2 * area).multiply(dot.divide(psr).multiply(specularReflectionCoeff).subtract(diffuseReflectionCoeff / 3));
  536.         final DerivativeStructure cS = dot.divide(psr).multiply(area * (specularReflectionCoeff - 1));
  537.         return new FieldVector3D<DerivativeStructure>(cN, normal, cS, fluxSat);

  538.     }

  539.     /** Compute contribution of one facet to force.
  540.      * <p>This method implements equation 8-44 from David A. Vallado's
  541.      * Fundamentals of Astrodynamics and Applications, third edition,
  542.      * 2007, Microcosm Press.</p>
  543.      * @param normal facet normal
  544.      * @param area facet area
  545.      * @param fluxSat radiation pressure flux in spacecraft frame
  546.      * @param dot dot product of facet and fluxSat (must be negative)
  547.      * @param specularReflectionCoeffDS specular reflection coefficient
  548.      * @param diffuseReflectionCoeffDS diffuse reflection coefficient
  549.      * @return contribution of the facet to force in spacecraft frame
  550.      */
  551.     private FieldVector3D<DerivativeStructure> facetRadiationAcceleration(final Vector3D normal, final double area,
  552.                                                                           final Vector3D fluxSat, final double dot,
  553.                                                                           final DerivativeStructure specularReflectionCoeffDS,
  554.                                                                           final DerivativeStructure diffuseReflectionCoeffDS) {
  555.         final double psr  = fluxSat.getNorm();

  556.         // Vallado's equation 8-44 uses different parameters which are related to our parameters as:
  557.         // cos (phi) = -dot / (psr * area)
  558.         // n         = facet / area
  559.         // s         = -fluxSat / psr
  560.         final DerivativeStructure cN =
  561.                 diffuseReflectionCoeffDS.divide(3).subtract(specularReflectionCoeffDS.multiply(dot / psr)).multiply(2 * area * dot);
  562.         final DerivativeStructure cS = specularReflectionCoeffDS.subtract(1).multiply(area * dot / psr);

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

  564.     }

  565.     /** Class representing a single facet of a convex spacecraft body.
  566.      * <p>Instance of this class are guaranteed to be immutable.</p>
  567.      * @author Luc Maisonobe
  568.      */
  569.     public static class Facet {

  570.         /** Unit Normal vector, pointing outward. */
  571.         private final Vector3D normal;

  572.         /** Area in m<sup>2</sup>. */
  573.         private final double area;

  574.         /** Simple constructor.
  575.          * @param normal vector normal to the facet, pointing outward (will be normalized)
  576.          * @param area facet area in m<sup>2</sup>
  577.          */
  578.         public Facet(final Vector3D normal, final double area) {
  579.             this.normal = normal.normalize();
  580.             this.area   = area;
  581.         }

  582.         /** Get unit normal vector.
  583.          * @return unit normal vector
  584.          */
  585.         public Vector3D getNormal() {
  586.             return normal;
  587.         }

  588.         /** Get facet area.
  589.          * @return facet area
  590.          */
  591.         public double getArea() {
  592.             return area;
  593.         }

  594.     }

  595.     /** Build the surface vectors for body facets of a simple parallelepipedic box.
  596.      * @param xLength length of the body along its X axis (m)
  597.      * @param yLength length of the body along its Y axis (m)
  598.      * @param zLength length of the body along its Z axis (m)
  599.      * @return surface vectors array
  600.      */
  601.     private static Facet[] simpleBoxFacets(final double xLength, final double yLength, final double zLength) {
  602.         return new Facet[] {
  603.             new Facet(Vector3D.MINUS_I, yLength * zLength),
  604.             new Facet(Vector3D.PLUS_I,  yLength * zLength),
  605.             new Facet(Vector3D.MINUS_J, xLength * zLength),
  606.             new Facet(Vector3D.PLUS_J,  xLength * zLength),
  607.             new Facet(Vector3D.MINUS_K, xLength * yLength),
  608.             new Facet(Vector3D.PLUS_K,  xLength * yLength)
  609.         };
  610.     }

  611.     /** Filter out zero area facets.
  612.      * @param facets original facets (may include zero area facets)
  613.      * @return filtered array
  614.      */
  615.     private static List<Facet> filter(final Facet[] facets) {
  616.         final List<Facet> filtered = new ArrayList<Facet>(facets.length);
  617.         for (Facet facet : facets) {
  618.             if (facet.getArea() > 0) {
  619.                 filtered.add(facet);
  620.             }
  621.         }
  622.         return filtered;
  623.     }

  624.     /** {@inheritDoc} */
  625.     public void setAbsorptionCoefficient(final double value) {
  626.         absorptionCoeff = value;
  627.         diffuseReflectionCoeff = 1 - (absorptionCoeff + specularReflectionCoeff);
  628.     }

  629.     /** {@inheritDoc} */
  630.     public double getAbsorptionCoefficient() {
  631.         return absorptionCoeff;
  632.     }

  633.     /** {@inheritDoc} */
  634.     public void setReflectionCoefficient(final double value) {
  635.         specularReflectionCoeff = value;
  636.         diffuseReflectionCoeff  = 1 - (absorptionCoeff + specularReflectionCoeff);
  637.     }

  638.     /** {@inheritDoc} */
  639.     public double getReflectionCoefficient() {
  640.         return specularReflectionCoeff;
  641.     }

  642.     /** {@inheritDoc} */
  643.     public void setDragCoefficient(final double value) {
  644.         dragCoeff = value;
  645.     }

  646.     /** {@inheritDoc} */
  647.     public double getDragCoefficient() {
  648.         return dragCoeff;
  649.     }

  650. }