NewcombOperators.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;
  21. import java.util.Map;
  22. import java.util.SortedMap;
  23. import java.util.TreeMap;

  25. import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
  26. import org.apache.commons.math3.analysis.polynomials.PolynomialsUtils;
  27. import org.apache.commons.math3.util.FastMath;

  29. /** Implementation of the Modified Newcomb Operators.
  30.  *  <p>
  31.  *  From equations 2.7.3 - (12)(13) of the Danielson paper, those operators are defined as:
  32.  *  <pre>
  33.  *  4(&rho; + &sigma;)Y<sub>&rho;,&sigma;</sub><sup>n,s</sup> =
  34.  *     2(2s - n)Y<sub>&rho;-1,&sigma;</sub><sup>n,s+1</sup> + (s - n)Y<sub>&rho;-2,&sigma;</sub><sup>n,s+2</sup>
  35.  *   - 2(2s + n)Y<sub>&rho;,&sigma;-1</sub><sup>n,s-1</sup> - (s+n)Y<sub>&rho;,&sigma;-2</sub><sup>n,s-2</sup>
  36.  *   + 2(2&rho; + 2&sigma; + 2 + 3n)Y<sub>&rho;-1,&sigma;-1</sub><sup>n,s</sup>
  37.  *  </pre>
  38.  *  Initialization is given by : <pre>Y<sub>0,0</sub><sup>n,s</sup> = 1</pre>
  39.  *  </p>
  40.  *
  41.  *  Internally, the Modified Newcomb Operators are stored as an array of {@link PolynomialFunction} :
  42.  *
  43.  *  <pre>
  44.  *  Y<sub>&rho;,&sigma;</sub><sup>n,s</sup> = P<sub>k<sub>0</sub></sub> + P<sub>k<sub>1</sub></sub>n + ... + P<sub>k<sub>j</sub></sub>n<sup>j</sup>
  45.  *  </pre>
  46.  *
  47.  * where the P<sub>k<sub>j</sub></sub> are given by
  48.  *
  49.  * <pre>
  50.  *  P<sub>k<sub>j</sub></sub> = &sum;<sub>j=0;&rho;</sub> a<sub>j</sub>s<sup>j</sup>
  51.  * </pre>
  52.  *
  53.  * @author Romain Di Costanzo
  54.  * @author Pascal Parraud
  55.  */
  56. public class NewcombOperators {

  58.     /** Storage map. */
  59.     private static final Map<NewKey, Double> MAP = new TreeMap<NewKey, Double>();

  61.     /** Private constructor as class is a utility.
  62.      */
  63.     private NewcombOperators() {
  64.     }

  66.     /** Get the Newcomb operator evaluated at n, s, &rho;, &sigma;.
  67.      * <p>
  68.      * This method is guaranteed to be thread-safe
  69.      * </p>
  70.      *  @param rho &rho; index
  71.      *  @param sigma &sigma; index
  72.      *  @param n n index
  73.      *  @param s s index
  74.      *  @return Y<sub>&rho;,&sigma;</sub><sup>n,s</sup>
  75.      */
  76.     public static double getValue(final int rho, final int sigma, final int n, final int s) {

  78.         final NewKey key = new NewKey(n, s, rho, sigma);
  79.         synchronized (MAP) {
  80.             if (MAP.containsKey(key)) {
  81.                 return MAP.get(key);
  82.             }
  83.         }

  85.         // Get the Newcomb polynomials for the given rho and sigma
  86.         final List<PolynomialFunction> polynomials = PolynomialsGenerator.getPolynomials(rho, sigma);

  88.         // Compute the value from the list of polynomials for the given n and s
  89.         double nPower = 1.;
  90.         double value = 0.0;
  91.         for (final PolynomialFunction polynomial : polynomials) {
  92.             value += polynomial.value(s) * nPower;
  93.             nPower = n * nPower;
  94.         }
  95.         synchronized (MAP) {
  96.             MAP.put(key, value);
  97.         }

  99.         return value;

  101.     }

  103.     /** Generator for Newcomb polynomials. */
  104.     private static class PolynomialsGenerator {

  106.         /** Polynomials storage. */
  107.         private static final SortedMap<Couple, List<PolynomialFunction>> POLYNOMIALS =
  108.                 new TreeMap<Couple, List<PolynomialFunction>>();

  110.         /** Private constructor as class is a utility.
  111.          */
  112.         private PolynomialsGenerator() {
  113.         }

  115.         /** Get the list of polynomials representing the Newcomb Operator for the (&rho;,&sigma;) couple.
  116.          * <p>
  117.          * This method is guaranteed to be thread-safe
  118.          * </p>
  119.          *  @param rho &rho; value
  120.          *  @param sigma &sigma; value
  121.          *  @return Polynomials representing the Newcomb Operator for the (&rho;,&sigma;) couple.
  122.          */
  123.         private static List<PolynomialFunction> getPolynomials(final int rho, final int sigma) {

  125.             final Couple couple = new Couple(rho, sigma);

  127.             synchronized (POLYNOMIALS) {

  129.                 if (POLYNOMIALS.isEmpty()) {
  130.                     // Initialize lists
  131.                     final List<PolynomialFunction> l00 = new ArrayList<PolynomialFunction>();
  132.                     final List<PolynomialFunction> l01 = new ArrayList<PolynomialFunction>();
  133.                     final List<PolynomialFunction> l10 = new ArrayList<PolynomialFunction>();
  134.                     final List<PolynomialFunction> l11 = new ArrayList<PolynomialFunction>();

  136.                     // Y(rho = 0, sigma = 0) = 1
  137.                     l00.add(new PolynomialFunction(new double[] {
  138.                         1.
  139.                     }));
  140.                     // Y(rho = 0, sigma = 1) =  -s - n/2
  141.                     l01.add(new PolynomialFunction(new double[] {
  142.                         0, -1.
  143.                     }));
  144.                     l01.add(new PolynomialFunction(new double[] {
  145.                         -0.5
  146.                     }));
  147.                     // Y(rho = 1, sigma = 0) =  s - n/2
  148.                     l10.add(new PolynomialFunction(new double[] {
  149.                         0, 1.
  150.                     }));
  151.                     l10.add(new PolynomialFunction(new double[] {
  152.                         -0.5
  153.                     }));
  154.                     // Y(rho = 1, sigma = 1) = 3/2 - s² + 5n/4 + n²/4
  155.                     l11.add(new PolynomialFunction(new double[] {
  156.                         1.5, 0., -1.
  157.                     }));
  158.                     l11.add(new PolynomialFunction(new double[] {
  159.                         1.25
  160.                     }));
  161.                     l11.add(new PolynomialFunction(new double[] {
  162.                         0.25
  163.                     }));

  165.                     // Initialize polynomials
  166.                     POLYNOMIALS.put(new Couple(0, 0), l00);
  167.                     POLYNOMIALS.put(new Couple(0, 1), l01);
  168.                     POLYNOMIALS.put(new Couple(1, 0), l10);
  169.                     POLYNOMIALS.put(new Couple(1, 1), l11);

  171.                 }

  173.                 // If order hasn't been computed yet, update the Newcomb polynomials
  174.                 if (!POLYNOMIALS.containsKey(couple)) {
  175.                     PolynomialsGenerator.computeFor(rho, sigma);
  176.                 }

  178.                 return POLYNOMIALS.get(couple);

  180.             }
  181.         }

  183.         /** Compute the Modified Newcomb Operators up to a given (&rho;, &sigma;) couple.
  184.          *  <p>
  185.          *  The recursive computation uses equation 2.7.3-(12) of the Danielson paper.
  186.          *  </p>
  187.          *  @param rho &rho; value to reach
  188.          *  @param sigma &sigma; value to reach
  189.          */
  190.         private static void computeFor(final int rho, final int sigma) {

  192.             // Initialize result :
  193.             List<PolynomialFunction> result = new ArrayList<PolynomialFunction>();

  195.             // Get the coefficient from the recurrence relation
  196.             final Map<Integer, List<PolynomialFunction>> map = generateRecurrenceCoefficients(rho, sigma);

  198.             // Compute (s - n) * Y[rho - 2, sigma][n, s + 2]
  199.             if (rho >= 2) {
  200.                 final List<PolynomialFunction> poly = map.get(0);
  201.                 final List<PolynomialFunction> list = getPolynomials(rho - 2, sigma);
  202.                 result = multiplyPolynomialList(poly, shiftList(list, 2));
  203.             }

  205.             // Compute 2(2rho + 2sigma + 2 + 3n) * Y[rho - 1, sigma - 1][n, s]
  206.             if (rho >= 1 && sigma >= 1) {
  207.                 final List<PolynomialFunction> poly = map.get(1);
  208.                 final List<PolynomialFunction> list = getPolynomials(rho - 1, sigma - 1);
  209.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, list));
  210.             }

  212.             // Compute 2(2s - n) * Y[rho - 1, sigma][n, s + 1]
  213.             if (rho >= 1) {
  214.                 final List<PolynomialFunction> poly = map.get(2);
  215.                 final List<PolynomialFunction> list = getPolynomials(rho - 1, sigma);
  216.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, 1)));
  217.             }

  219.             // Compute -(s + n) * Y[rho, sigma - 2][n, s - 2]
  220.             if (sigma >= 2) {
  221.                 final List<PolynomialFunction> poly = map.get(3);
  222.                 final List<PolynomialFunction> list = getPolynomials(rho, sigma - 2);
  223.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, -2)));
  224.             }

  226.             // Compute -2(2s + n) * Y[rho, sigma - 1][n, s - 1]
  227.             if (sigma >= 1) {
  228.                 final List<PolynomialFunction> poly = map.get(4);
  229.                 final List<PolynomialFunction> list = getPolynomials(rho, sigma - 1);
  230.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, -1)));
  231.             }

  233.             // Save polynomials for current (rho, sigma) couple
  234.             final Couple couple = new Couple(rho, sigma);
  235.             POLYNOMIALS.put(couple, result);
  236.         }

  238.         /** Multiply two lists of polynomials defined as the internal representation of the Newcomb Operator.
  239.          *  <p>
  240.          *  Let's call R<sub>s</sub>(n) the result returned by the method :
  241.          *  <pre>
  242.          *  R<sub>s</sub>(n) = (P<sub>s<sub>0</sub></sub> + P<sub>s<sub>1</sub></sub>n + ... + P<sub>s<sub>j</sub></sub>n<sup>j</sup>) *(Q<sub>s<sub>0</sub></sub> + Q<sub>s<sub>1</sub></sub>n + ... + Q<sub>s<sub>k</sub></sub>n<sup>k</sup>)
  243.          *  </pre>
  244.          *
  245.          *  where the P<sub>s<sub>j</sub></sub> and Q<sub>s<sub>k</sub></sub> are polynomials in s,
  246.          *  s being the index of the Y<sub>&rho;,&sigma;</sub><sup>n,s</sup> function
  247.          *
  248.          *  @param poly1 first list of polynomials
  249.          *  @param poly2 second list of polynomials
  250.          *  @return R<sub>s</sub>(n) as a list of {@link PolynomialFunction}
  251.          */
  252.         private static List<PolynomialFunction> multiplyPolynomialList(final List<PolynomialFunction> poly1,
  253.                                                                        final List<PolynomialFunction> poly2) {
  254.             // Initialize the result list of polynomial function
  255.             final List<PolynomialFunction> result = new ArrayList<PolynomialFunction>();
  256.             initializeListOfPolynomials(poly1.size() + poly2.size() - 1, result);

  258.             int i = 0;
  259.             // Iterate over first polynomial list
  260.             for (PolynomialFunction f1 : poly1) {
  261.                 // Iterate over second polynomial list
  262.                 for (int j = i; j < poly2.size() + i; j++) {
  263.                     final PolynomialFunction p2 = poly2.get(j - i);
  264.                     // Get previous polynomial for current 'n' order
  265.                     final PolynomialFunction previousP2 = result.get(j);
  266.                     // Replace the current order by summing the product of both of the polynomials
  267.                     result.set(j, previousP2.add(f1.multiply(p2)));
  268.                 }
  269.                 // shift polynomial order in 'n'
  270.                 i++;
  271.             }
  272.             return result;
  273.         }

  275.         /** Sum two lists of {@link PolynomialFunction}.
  276.          *  @param poly1 first list
  277.          *  @param poly2 second list
  278.          *  @return the summed list
  279.          */
  280.         private static List<PolynomialFunction> sumPolynomialList(final List<PolynomialFunction> poly1,
  281.                                                                   final List<PolynomialFunction> poly2) {
  282.             // identify the lowest degree polynomial
  283.             final int lowLength  = FastMath.min(poly1.size(), poly2.size());
  284.             final int highLength = FastMath.max(poly1.size(), poly2.size());
  285.             // Initialize the result list of polynomial function
  286.             final List<PolynomialFunction> result = new ArrayList<PolynomialFunction>();
  287.             initializeListOfPolynomials(highLength, result);

  289.             for (int i = 0; i < lowLength; i++) {
  290.                 // Add polynomials by increasing order of 'n'
  291.                 result.set(i, poly1.get(i).add(poly2.get(i)));
  292.             }
  293.             // Complete the list if lists are of different size:
  294.             for (int i = lowLength; i < highLength; i++) {
  295.                 if (poly1.size() < poly2.size()) {
  296.                     result.set(i, poly2.get(i));
  297.                 } else {
  298.                     result.set(i, poly1.get(i));
  299.                 }
  300.             }
  301.             return result;
  302.         }

  304.         /** Initialize an empty list of polynomials.
  305.          *  @param i order
  306.          *  @param result list into which polynomials should be added
  307.          */
  308.         private static void initializeListOfPolynomials(final int i,
  309.                                                         final List<PolynomialFunction> result) {
  310.             for (int k = 0; k < i; k++) {
  311.                 result.add(new PolynomialFunction(new double[] {0.}));
  312.             }
  313.         }

  315.         /** Shift a list of {@link PolynomialFunction}, from the
  316.          *  {@link PolynomialsUtils#shift(double[], double)} method.
  317.          *
  318.          *  @param polynomialList list of {@link PolynomialFunction}
  319.          *  @param shift shift value
  320.          *  @return new list of shifted {@link PolynomialFunction}
  321.          */
  322.         private static List<PolynomialFunction> shiftList(final List<PolynomialFunction> polynomialList,
  323.                                                           final int shift) {
  324.             final List<PolynomialFunction> shiftedList = new ArrayList<PolynomialFunction>();
  325.             for (PolynomialFunction function : polynomialList) {
  326.                 shiftedList.add(new PolynomialFunction(PolynomialsUtils.shift(function.getCoefficients(), shift)));
  327.             }
  328.             return shiftedList;
  329.         }

  331.         /** Generate recurrence coefficients.
  332.          *
  333.          *  @param rho &rho; value
  334.          *  @param sigma &sigma; value
  335.          *  @return recurrence coefficients
  336.          */
  337.         private static Map<Integer, List<PolynomialFunction>> generateRecurrenceCoefficients(final int rho, final int sigma) {
  338.             final double den   = 1. / (4. * (rho + sigma));
  339.             final double denx2 = 2. * den;
  340.             final double denx4 = 4. * den;
  341.             // Initialization :
  342.             final Map<Integer, List<PolynomialFunction>> list = new TreeMap<Integer, List<PolynomialFunction>>();
  343.             final List<PolynomialFunction> poly0 = new ArrayList<PolynomialFunction>();
  344.             final List<PolynomialFunction> poly1 = new ArrayList<PolynomialFunction>();
  345.             final List<PolynomialFunction> poly2 = new ArrayList<PolynomialFunction>();
  346.             final List<PolynomialFunction> poly3 = new ArrayList<PolynomialFunction>();
  347.             final List<PolynomialFunction> poly4 = new ArrayList<PolynomialFunction>();
  348.             // (s - n)
  349.             poly0.add(new PolynomialFunction(new double[] {0., den}));
  350.             poly0.add(new PolynomialFunction(new double[] {-den}));
  351.             // 2(2 * rho + 2 * sigma + 2 + 3*n)
  352.             poly1.add(new PolynomialFunction(new double[] {1. + denx4}));
  353.             poly1.add(new PolynomialFunction(new double[] {denx2 + denx4}));
  354.             // 2(2s - n)
  355.             poly2.add(new PolynomialFunction(new double[] {0., denx4}));
  356.             poly2.add(new PolynomialFunction(new double[] {-denx2}));
  357.             // - (s + n)
  358.             poly3.add(new PolynomialFunction(new double[] {0., -den}));
  359.             poly3.add(new PolynomialFunction(new double[] {-den}));
  360.             // - 2(2s + n)
  361.             poly4.add(new PolynomialFunction(new double[] {0., -denx4}));
  362.             poly4.add(new PolynomialFunction(new double[] {-denx2}));

  364.             // Fill the map :
  365.             list.put(0, poly0);
  366.             list.put(1, poly1);
  367.             list.put(2, poly2);
  368.             list.put(3, poly3);
  369.             list.put(4, poly4);
  370.             return list;
  371.         }

  373.     }

  375.     /** Private class to define a couple of value. */
  376.     private static class Couple implements Comparable<Couple> {

  378.         /** first couple value. */
  379.         private final int rho;

  381.         /** second couple value. */
  382.         private final int sigma;

  384.         /** Constructor.
  385.          * @param rho first couple value
  386.          * @param sigma second couple value
  387.          */
  388.         private Couple(final int rho, final int sigma) {
  389.             this.rho = rho;
  390.             this.sigma = sigma;
  391.         }

  393.         /** {@inheritDoc} */
  394.         public int compareTo(final Couple c) {
  395.             int result = 1;
  396.             if (rho == c.rho) {
  397.                 if (sigma < c.sigma) {
  398.                     result = -1;
  399.                 } else if (sigma == c.sigma) {
  400.                     result = 0;
  401.                 }
  402.             } else if (rho < c.rho) {
  403.                 result = -1;
  404.             }
  405.             return result;
  406.         }

  408.         /** {@inheritDoc} */
  409.         public boolean equals(final Object couple) {

  411.             if (couple == this) {
  412.                 // first fast check
  413.                 return true;
  414.             }

  416.             if ((couple != null) && (couple instanceof Couple)) {
  417.                 return (rho == ((Couple) couple).rho) && (sigma == ((Couple) couple).sigma);
  418.             }

  420.             return false;

  422.         }

  424.         /** {@inheritDoc} */
  425.         public int hashCode() {
  426.             return 0x7ab17c0c ^ (rho << 8) ^ sigma;
  427.         }

  429.     }

  431.     /** Newcomb operator's key. */
  432.     private static class NewKey implements Comparable<NewKey> {

  434.         /** n value. */
  435.         private final int n;

  437.         /** s value. */
  438.         private final int s;

  440.         /** &rho; value. */
  441.         private final int rho;

  443.         /** &sigma; value. */
  444.         private final int sigma;

  446.         /** Simpleconstructor.
  447.          * @param n n
  448.          * @param s s
  449.          * @param rho &rho;
  450.          * @param sigma &sigma;
  451.          */
  452.         public NewKey(final int n, final int s, final int rho, final int sigma) {
  453.             this.n = n;
  454.             this.s = s;
  455.             this.rho = rho;
  456.             this.sigma = sigma;
  457.         }

  459.         /** {@inheritDoc} */
  460.         public int compareTo(final NewKey key) {
  461.             int result = 1;
  462.             if (n == key.n) {
  463.                 if (s == key.s) {
  464.                     if (rho == key.rho) {
  465.                         if (sigma < key.sigma) {
  466.                             result = -1;
  467.                         } else if (sigma == key.sigma) {
  468.                             result = 0;
  469.                         } else {
  470.                             result = 1;
  471.                         }
  472.                     } else if (rho < key.rho) {
  473.                         result = -1;
  474.                     } else {
  475.                         result = 1;
  476.                     }
  477.                 } else if (s < key.s) {
  478.                     result = -1;
  479.                 } else {
  480.                     result = 1;
  481.                 }
  482.             } else if (n < key.n) {
  483.                 result = -1;
  484.             }
  485.             return result;
  486.         }

  488.         /** {@inheritDoc} */
  489.         public boolean equals(final Object key) {

  491.             if (key == this) {
  492.                 // first fast check
  493.                 return true;
  494.             }

  496.             if ((key != null) && (key instanceof NewKey)) {
  497.                 return (n     == ((NewKey) key).n) &&
  498.                        (s     == ((NewKey) key).s) &&
  499.                        (rho   == ((NewKey) key).rho) &&
  500.                        (sigma == ((NewKey) key).sigma);
  501.             }

  503.             return false;

  505.         }

  507.         /** {@inheritDoc} */
  508.         public int hashCode() {
  509.             return 0x25baa451 ^ (n << 24) ^ (s << 16) ^ (rho << 8) ^ sigma;
  510.         }

  512.     }

  514. }
Created with JaCoCo 0.6.3.201306030806