GHmsjPolynomials.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.propagation.semianalytical.dsst.utilities;

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

  22. import org.apache.commons.math3.complex.Complex;
  23. import org.apache.commons.math3.util.FastMath;

  25. /** Compute the G<sub>ms</sub><sup>j</sup> and the H<sub>ms</sub><sup>j</sup>
  26.  *  polynomials in the equinoctial elements h, k and the direction cosines &alpha; and &beta;
  27.  *  and their partial derivatives with respect to k, h, &alpha; and &beta;.
  28.  *  <p>
  29.  *  The expressions used are equations 2.7.5-(1)(2) from the Danielson paper.
  30.  *  </p>
  31.  *  @author Romain Di Costanzo
  32.  */
  33. public class GHmsjPolynomials {

  35.     /** C<sub>j</sub>(k, h), S<sub>j</sub>(k, h) coefficient.
  36.      * (k, h) are the (x, y) component of the eccentricity vector in equinoctial elements
  37.      */
  38.     private final CjSjCoefficient cjsjKH;

  40.     /** C<sub>j</sub>(&alpha;, &beta;), S<sub>j</sub>(&alpha;, &beta;) coefficient.
  41.      * (&alpha;, &beta;) are the direction cosines
  42.      */
  43.     private final CjSjCoefficient cjsjAB;

  45.     /** Is the orbit represented as a retrograde orbit.
  46.      *  I = -1 if yes, +1 otherwise.
  47.      */
  48.     private int                   I;

  50.     /** Create a set of G<sub>ms</sub><sup>j</sup> and H<sub>ms</sub><sup>j</sup> polynomials.
  51.      *  @param k X component of the eccentricity vector
  52.      *  @param h Y component of the eccentricity vector
  53.      *  @param alpha direction cosine &alpha;
  54.      *  @param beta direction cosine &beta;
  55.      *  @param retroFactor -1 if the orbit is represented as retrograde, +1 otherwise
  56.      **/
  57.     public GHmsjPolynomials(final double k, final double h,
  58.                             final double alpha, final double beta,
  59.                             final int retroFactor) {
  60.         this.cjsjKH = new CjSjCoefficient(k, h);
  61.         this.cjsjAB = new CjSjCoefficient(alpha, beta);
  62.         this.I = retroFactor;
  63.     }

  65.     /** Get the G<sub>ms</sub><sup>j</sup> coefficient.
  66.      * @param m m subscript
  67.      * @param s s subscript
  68.      * @param j order
  69.      * @return the G<sub>ms</sub><sup>j</sup>
  70.      */
  71.     public double getGmsj(final int m, final int s, final int j) {
  72.         final int sMj = FastMath.abs(s - j);
  73.         double gms = 0d;
  74.         if (FastMath.abs(s) <= m) {
  75.             final int mMis = m - I * s;
  76.             gms = cjsjKH.getCj(sMj) * cjsjAB.getCj(mMis) -
  77.                   I * sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getSj(mMis);
  78.         } else {
  79.             final int sMim = FastMath.abs(s - I * m);
  80.             gms = cjsjKH.getCj(sMj) * cjsjAB.getCj(sMim) +
  81.                   sgn(s - j) * sgn(s - m) * cjsjKH.getSj(sMj) * cjsjAB.getSj(sMim);
  82.         }
  83.         return gms;
  84.     }

  86.     /** Get the H<sub>ms</sub><sup>j</sup> coefficient.
  87.      * @param m m subscript
  88.      * @param s s subscript
  89.      * @param j order
  90.      * @return the H<sub>ms</sub><sup>j</sup>
  91.      */
  92.     public double getHmsj(final int m, final int s, final int j) {
  93.         final int sMj = FastMath.abs(s - j);
  94.         double hms = 0d;
  95.         if (FastMath.abs(s) <= m) {
  96.             final int mMis = m - I * s;
  97.             hms = I * cjsjKH.getCj(sMj) * cjsjAB.getSj(mMis) + sgn(s - j) *
  98.                   cjsjKH.getSj(sMj) * cjsjAB.getCj(mMis);
  99.         } else {
  100.             final int sMim = FastMath.abs(s - I * m);
  101.             hms = -sgn(s - m) * cjsjKH.getCj(sMj) * cjsjAB.getSj(sMim) +
  102.                    sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getCj(sMim);
  103.         }
  104.         return hms;
  105.     }

  107.     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>k</sub> coefficient.
  108.      * @param m m subscript
  109.      * @param s s subscript
  110.      * @param j order
  111.      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>k</sub>
  112.      */
  113.     public double getdGmsdk(final int m, final int s, final int j) {
  114.         final int sMj = FastMath.abs(s - j);
  115.         double dGmsdk = 0d;
  116.         if (FastMath.abs(s) <= m) {
  117.             final int mMis = m - I * s;
  118.             dGmsdk = cjsjKH.getDcjDk(sMj) * cjsjAB.getCj(mMis) -
  119.                    I * sgn(s - j) * cjsjKH.getDsjDk(sMj) * cjsjAB.getSj(mMis);
  120.         } else {
  121.             final int sMim = FastMath.abs(s - I * m);
  122.             dGmsdk = cjsjKH.getDcjDk(sMj) * cjsjAB.getCj(sMim) +
  123.                     sgn(s - j) * sgn(s - m) * cjsjKH.getDsjDk(sMj) * cjsjAB.getSj(sMim);
  124.         }
  125.         return dGmsdk;
  126.     }

  128.     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>h</sub> coefficient.
  129.      * @param m m subscript
  130.      * @param s s subscript
  131.      * @param j order
  132.      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>h</sub>
  133.      */
  134.     public double getdGmsdh(final int m, final int s, final int j) {
  135.         final int sMj = FastMath.abs(s - j);
  136.         double dGmsdh = 0d;
  137.         if (FastMath.abs(s) <= m) {
  138.             final int mMis = m - I * s;
  139.             dGmsdh = cjsjKH.getDcjDh(sMj) * cjsjAB.getCj(mMis) -
  140.                     I * sgn(s - j) * cjsjKH.getDsjDh(sMj) * cjsjAB.getSj(mMis);
  141.         } else {
  142.             final int sMim = FastMath.abs(s - I * m);
  143.             dGmsdh = cjsjKH.getDcjDh(sMj) * cjsjAB.getCj(sMim) +
  144.                     sgn(s - j) * sgn(s - m) * cjsjKH.getDsjDh(sMj) * cjsjAB.getSj(sMim);
  145.         }
  146.         return dGmsdh;
  147.     }

  149.     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub> coefficient.
  150.      * @param m m subscript
  151.      * @param s s subscript
  152.      * @param j order
  153.      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub>
  154.      */
  155.     public double getdGmsdAlpha(final int m, final int s, final int j) {
  156.         final int sMj  = FastMath.abs(s - j);
  157.         double dGmsdAl = 0d;
  158.         if (FastMath.abs(s) <= m) {
  159.             final int mMis = m - I * s;
  160.             dGmsdAl = cjsjKH.getCj(sMj) * cjsjAB.getDcjDk(mMis) -
  161.                    I * sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDk(mMis);
  162.         } else {
  163.             final int sMim = FastMath.abs(s - I * m);
  164.             dGmsdAl = cjsjKH.getCj(sMj) * cjsjAB.getDcjDk(sMim) +
  165.                     sgn(s - j) * sgn(s - m) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDk(sMim);
  166.         }
  167.         return dGmsdAl;
  168.     }

  170.     /** Get the dG<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub> coefficient.
  171.      * @param m m subscript
  172.      * @param s s subscript
  173.      * @param j order
  174.      * @return dG<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub>
  175.      */
  176.     public double getdGmsdBeta(final int m, final int s, final int j) {
  177.         final int sMj = FastMath.abs(s - j);
  178.         double dGmsdBe = 0d;
  179.         if (FastMath.abs(s) <= m) {
  180.             final int mMis = m - I * s;
  181.             dGmsdBe = cjsjKH.getCj(sMj) * cjsjAB.getDcjDh(mMis) -
  182.                     I * sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDh(mMis);
  183.         } else {
  184.             final int sMim = FastMath.abs(s - I * m);
  185.             dGmsdBe = cjsjKH.getCj(sMj) * cjsjAB.getDcjDh(sMim) +
  186.                     sgn(s - j) * sgn(s - m) * cjsjKH.getSj(sMj) * cjsjAB.getDsjDh(sMim);
  187.         }
  188.         return dGmsdBe;
  189.     }

  191.     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>k</sub> coefficient.
  192.      * @param m m subscript
  193.      * @param s s subscript
  194.      * @param j order
  195.      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>k</sub>
  196.      */
  197.     public double getdHmsdk(final int m, final int s, final int j) {
  198.         final int sMj = FastMath.abs(s - j);
  199.         double dHmsdk = 0d;
  200.         if (FastMath.abs(s) <= m) {
  201.             final int mMis = m - I * s;
  202.             dHmsdk = I * cjsjKH.getDcjDk(sMj) * cjsjAB.getSj(mMis) +
  203.                     sgn(s - j) * cjsjKH.getDsjDk(sMj) * cjsjAB.getCj(mMis);
  204.         } else {
  205.             final int sMim = FastMath.abs(s - I * m);
  206.             dHmsdk = -sgn(s - m) * cjsjKH.getDcjDk(sMj) * cjsjAB.getSj(sMim) +
  207.                     sgn(s - j) * cjsjKH.getDsjDk(sMj) * cjsjAB.getCj(sMim);
  208.         }
  209.         return dHmsdk;
  210.     }

  212.     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>h</sub> coefficient.
  213.      * @param m m subscript
  214.      * @param s s subscript
  215.      * @param j order
  216.      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>h</sub>
  217.      */
  218.     public double getdHmsdh(final int m,  final int s, final int j) {
  219.         final int sMj = FastMath.abs(s - j);
  220.         double dHmsdh = 0d;
  221.         if (FastMath.abs(s) <= m) {
  222.             final int mMis = m - I * s;
  223.             dHmsdh = I * cjsjKH.getDcjDh(sMj) * cjsjAB.getSj(mMis) +
  224.                     sgn(s - j) * cjsjKH.getDsjDh(sMj) * cjsjAB.getCj(mMis);
  225.         } else {
  226.             final int sMim = FastMath.abs(s - I * m);
  227.             dHmsdh = -sgn(s - m) * cjsjKH.getDcjDh(sMj) * cjsjAB.getSj(sMim) +
  228.                     sgn(s - j) * cjsjKH.getDsjDh(sMj) * cjsjAB.getCj(sMim);
  229.         }
  230.         return dHmsdh;
  231.     }

  233.     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub> coefficient.
  234.      * @param m m subscript
  235.      * @param s s subscript
  236.      * @param j order
  237.      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>&alpha;</sub>
  238.      */
  239.     public double getdHmsdAlpha(final int m, final int s, final int j) {
  240.         final int sMj  = FastMath.abs(s - j);
  241.         double dHmsdAl = 0d;
  242.         if (FastMath.abs(s) <= m) {
  243.             final int mMis = m - I * s;
  244.             dHmsdAl = I * cjsjKH.getCj(sMj) * cjsjAB.getDsjDk(mMis) +
  245.                     sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDk(mMis);
  246.         } else {
  247.             final int sMim = FastMath.abs(s - I * m);
  248.             dHmsdAl = -sgn(s - m) * cjsjKH.getCj(sMj) * cjsjAB.getDsjDk(sMim) +
  249.                     sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDk(sMim);
  250.         }
  251.         return dHmsdAl;
  252.     }

  254.     /** Get the dH<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub> coefficient.
  255.      * @param m m subscript
  256.      * @param s s subscript
  257.      * @param j order
  258.      * @return dH<sub>ms</sub><sup>j</sup> / d<sub>&beta;</sub>
  259.      */
  260.     public double getdHmsdBeta(final int m, final int s, final int j) {
  261.         final int sMj = FastMath.abs(s - j);
  262.         double dHmsdBe = 0d;
  263.         if (FastMath.abs(s) <= m) {
  264.             final int mMis = m - I * s;
  265.             dHmsdBe = I * cjsjKH.getCj(sMj) * cjsjAB.getDsjDh(mMis) +
  266.                    sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDh(mMis);
  267.         } else {
  268.             final int sMim = FastMath.abs(s - I * m);
  269.             dHmsdBe = -sgn(s - m) * cjsjKH.getCj(sMj) * cjsjAB.getDsjDh(sMim) +
  270.                     sgn(s - j) * cjsjKH.getSj(sMj) * cjsjAB.getDcjDh(sMim);
  271.         }
  272.         return dHmsdBe;
  273.     }

  275.     /** Get the sign of an integer.
  276.      *  @param i number on which evaluation is done
  277.      *  @return -1 or +1 depending on sign of i
  278.      */
  279.     private int sgn(final int i) {
  280.         return (i < 0) ? -1 : 1;
  281.     }

  283.     /** Compute the S<sub>j</sub>(k, h) and the C<sub>j</sub>(k, h) series
  284.      *  and their partial derivatives with respect to k and h.
  285.      *  <p>
  286.      *  Those series are given in Danielson paper by expression 2.5.3-(5):
  287.      *  <pre>C<sub>j</sub>(k, h) + i S<sub>j</sub>(k, h) = (k+ih)<sup>j</sup> </pre>
  288.      *  </p>
  289.      *  The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) elements are store as an
  290.      *  {@link ArrayList} of {@link Complex} number, the C<sub>j</sub>(k, h) being
  291.      *  represented by the real and the S<sub>j</sub>(k, h) by the imaginary part.
  292.      */
  293.     private static class CjSjCoefficient {

  295.         /** Last computed order j. */
  296.         private int jLast;

  298.         /** Complex base (k + ih) of the C<sub>j</sub>, S<sub>j</sub> series. */
  299.         private final Complex kih;

  301.         /** List of computed elements. */
  302.         private final List<Complex> cjsj;

  304.         /** C<sub>j</sub>(k, h) and S<sub>j</sub>(k, h) constructor.
  305.          * @param k k value
  306.          * @param h h value
  307.          */
  308.         public CjSjCoefficient(final double k, final double h) {
  309.             kih  = new Complex(k, h);
  310.             cjsj = new ArrayList<Complex>();
  311.             cjsj.add(new Complex(1, 0));
  312.             cjsj.add(kih);
  313.             jLast = 1;
  314.         }

  316.         /** Get the C<sub>j</sub> coefficient.
  317.          * @param j order
  318.          * @return C<sub>j</sub>
  319.          */
  320.         public double getCj(final int j) {
  321.             if (j > jLast) {
  322.                 // Update to order j
  323.                 updateCjSj(j);
  324.             }
  325.             return cjsj.get(j).getReal();
  326.         }

  328.         /** Get the S<sub>j</sub> coefficient.
  329.          * @param j order
  330.          * @return S<sub>j</sub>
  331.          */
  332.         public double getSj(final int j) {
  333.             if (j > jLast) {
  334.                 // Update to order j
  335.                 updateCjSj(j);
  336.             }
  337.             return cjsj.get(j).getImaginary();
  338.         }

  340.         /** Get the dC<sub>j</sub> / dk coefficient.
  341.          * @param j order
  342.          * @return dC<sub>j</sub> / d<sub>k</sub>
  343.          */
  344.         public double getDcjDk(final int j) {
  345.             return j == 0 ? 0 : j * getCj(j - 1);
  346.         }

  348.         /** Get the dS<sub>j</sub> / dk coefficient.
  349.          * @param j order
  350.          * @return dS<sub>j</sub> / d<sub>k</sub>
  351.          */
  352.         public double getDsjDk(final int j) {
  353.             return j == 0 ? 0 : j * getSj(j - 1);
  354.         }

  356.         /** Get the dC<sub>j</sub> / dh coefficient.
  357.          * @param j order
  358.          * @return dC<sub>i</sub> / d<sub>k</sub>
  359.          */
  360.         public double getDcjDh(final int j) {
  361.             return j == 0 ? 0 : -j * getSj(j - 1);
  362.         }

  364.         /** Get the dS<sub>j</sub> / dh coefficient.
  365.          * @param j order
  366.          * @return dS<sub>j</sub> / d<sub>h</sub>
  367.          */
  368.         public double getDsjDh(final int j) {
  369.             return j == 0 ? 0 : j * getCj(j - 1);
  370.         }

  372.         /** Update the cjsj up to order j.
  373.          * @param j order
  374.          */
  375.         private void updateCjSj(final int j) {
  376.             Complex last = cjsj.get(cjsj.size() - 1);
  377.             for (int i = jLast; i < j; i++) {
  378.                 final Complex next = last.multiply(kih);
  379.                 cjsj.add(next);
  380.                 last = next;
  381.             }
  382.             jLast = j;
  383.         }

  385.     }

  387. }
Created with JaCoCo 0.6.3.201306030806