NewcombOperators.java

  1. /* Copyright 2002-2020 CS Group
  2.  * Licensed to CS Group (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;

  18. import java.util.ArrayList;
  19. import java.util.List;
  20. import java.util.Map;
  21. import java.util.SortedMap;
  22. import java.util.TreeMap;

  23. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  24. import org.hipparchus.util.FastMath;

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

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

  55.     /** Private constructor as class is a utility.
  56.      */
  57.     private NewcombOperators() {
  58.     }

  59.     /** Get the Newcomb operator evaluated at n, s, ρ, σ.
  60.      * <p>
  61.      * This method is guaranteed to be thread-safe
  62.      * </p>
  63.      *  @param rho ρ index
  64.      *  @param sigma σ index
  65.      *  @param n n index
  66.      *  @param s s index
  67.      *  @return Y<sub>ρ,σ</sub><sup>n,s</sup>
  68.      */
  69.     public static double getValue(final int rho, final int sigma, final int n, final int s) {

  70.         final NewKey key = new NewKey(n, s, rho, sigma);
  71.         synchronized (MAP) {
  72.             if (MAP.containsKey(key)) {
  73.                 return MAP.get(key);
  74.             }
  75.         }

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

  78.         // Compute the value from the list of polynomials for the given n and s
  79.         double nPower = 1.;
  80.         double value = 0.0;
  81.         for (final PolynomialFunction polynomial : polynomials) {
  82.             value += polynomial.value(s) * nPower;
  83.             nPower = n * nPower;
  84.         }
  85.         synchronized (MAP) {
  86.             MAP.put(key, value);
  87.         }

  88.         return value;

  89.     }

  90.     /** Generator for Newcomb polynomials. */
  91.     private static class PolynomialsGenerator {

  92.         /** Polynomials storage. */
  93.         private static final SortedMap<Couple, List<PolynomialFunction>> POLYNOMIALS =
  94.                 new TreeMap<Couple, List<PolynomialFunction>>();

  95.         /** Private constructor as class is a utility.
  96.          */
  97.         private PolynomialsGenerator() {
  98.         }

  99.         /** Get the list of polynomials representing the Newcomb Operator for the (ρ,σ) couple.
  100.          * <p>
  101.          * This method is guaranteed to be thread-safe
  102.          * </p>
  103.          *  @param rho ρ value
  104.          *  @param sigma σ value
  105.          *  @return Polynomials representing the Newcomb Operator for the (ρ,σ) couple.
  106.          */
  107.         private static List<PolynomialFunction> getPolynomials(final int rho, final int sigma) {

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

  109.             synchronized (POLYNOMIALS) {

  110.                 if (POLYNOMIALS.isEmpty()) {
  111.                     // Initialize lists
  112.                     final List<PolynomialFunction> l00 = new ArrayList<PolynomialFunction>();
  113.                     final List<PolynomialFunction> l01 = new ArrayList<PolynomialFunction>();
  114.                     final List<PolynomialFunction> l10 = new ArrayList<PolynomialFunction>();
  115.                     final List<PolynomialFunction> l11 = new ArrayList<PolynomialFunction>();

  116.                     // Y(rho = 0, sigma = 0) = 1
  117.                     l00.add(new PolynomialFunction(new double[] {
  118.                         1.
  119.                     }));
  120.                     // Y(rho = 0, sigma = 1) =  -s - n/2
  121.                     l01.add(new PolynomialFunction(new double[] {
  122.                         0, -1.
  123.                     }));
  124.                     l01.add(new PolynomialFunction(new double[] {
  125.                         -0.5
  126.                     }));
  127.                     // Y(rho = 1, sigma = 0) =  s - n/2
  128.                     l10.add(new PolynomialFunction(new double[] {
  129.                         0, 1.
  130.                     }));
  131.                     l10.add(new PolynomialFunction(new double[] {
  132.                         -0.5
  133.                     }));
  134.                     // Y(rho = 1, sigma = 1) = 3/2 - s² + 5n/4 + n²/4
  135.                     l11.add(new PolynomialFunction(new double[] {
  136.                         1.5, 0., -1.
  137.                     }));
  138.                     l11.add(new PolynomialFunction(new double[] {
  139.                         1.25
  140.                     }));
  141.                     l11.add(new PolynomialFunction(new double[] {
  142.                         0.25
  143.                     }));

  144.                     // Initialize polynomials
  145.                     POLYNOMIALS.put(new Couple(0, 0), l00);
  146.                     POLYNOMIALS.put(new Couple(0, 1), l01);
  147.                     POLYNOMIALS.put(new Couple(1, 0), l10);
  148.                     POLYNOMIALS.put(new Couple(1, 1), l11);

  149.                 }

  150.                 // If order hasn't been computed yet, update the Newcomb polynomials
  151.                 if (!POLYNOMIALS.containsKey(couple)) {
  152.                     PolynomialsGenerator.computeFor(rho, sigma);
  153.                 }

  154.                 return POLYNOMIALS.get(couple);

  155.             }
  156.         }

  157.         /** Compute the Modified Newcomb Operators up to a given (ρ, σ) couple.
  158.          *  <p>
  159.          *  The recursive computation uses equation 2.7.3-(12) of the Danielson paper.
  160.          *  </p>
  161.          *  @param rho ρ value to reach
  162.          *  @param sigma σ value to reach
  163.          */
  164.         private static void computeFor(final int rho, final int sigma) {

  165.             // Initialize result :
  166.             List<PolynomialFunction> result = new ArrayList<PolynomialFunction>();

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

  169.             // Compute (s - n) * Y[rho - 2, sigma][n, s + 2]
  170.             if (rho >= 2) {
  171.                 final List<PolynomialFunction> poly = map.get(0);
  172.                 final List<PolynomialFunction> list = getPolynomials(rho - 2, sigma);
  173.                 result = multiplyPolynomialList(poly, shiftList(list, 2));
  174.             }

  175.             // Compute 2(2rho + 2sigma + 2 + 3n) * Y[rho - 1, sigma - 1][n, s]
  176.             if (rho >= 1 && sigma >= 1) {
  177.                 final List<PolynomialFunction> poly = map.get(1);
  178.                 final List<PolynomialFunction> list = getPolynomials(rho - 1, sigma - 1);
  179.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, list));
  180.             }

  181.             // Compute 2(2s - n) * Y[rho - 1, sigma][n, s + 1]
  182.             if (rho >= 1) {
  183.                 final List<PolynomialFunction> poly = map.get(2);
  184.                 final List<PolynomialFunction> list = getPolynomials(rho - 1, sigma);
  185.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, 1)));
  186.             }

  187.             // Compute -(s + n) * Y[rho, sigma - 2][n, s - 2]
  188.             if (sigma >= 2) {
  189.                 final List<PolynomialFunction> poly = map.get(3);
  190.                 final List<PolynomialFunction> list = getPolynomials(rho, sigma - 2);
  191.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, -2)));
  192.             }

  193.             // Compute -2(2s + n) * Y[rho, sigma - 1][n, s - 1]
  194.             if (sigma >= 1) {
  195.                 final List<PolynomialFunction> poly = map.get(4);
  196.                 final List<PolynomialFunction> list = getPolynomials(rho, sigma - 1);
  197.                 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, -1)));
  198.             }

  199.             // Save polynomials for current (rho, sigma) couple
  200.             final Couple couple = new Couple(rho, sigma);
  201.             POLYNOMIALS.put(couple, result);
  202.         }

  203.         /** Multiply two lists of polynomials defined as the internal representation of the Newcomb Operator.
  204.          *  <p>
  205.          *  Let's call R<sub>s</sub>(n) the result returned by the method :
  206.          *  <pre>
  207.          *  R<sub>s</sub>(n) = (P<sub>s₀</sub> + P<sub>s₁</sub>n + ... + P<sub>s<sub>j</sub></sub>n<sup>j</sup>) *(Q<sub>s₀</sub> + Q<sub>s₁</sub>n + ... + Q<sub>s<sub>k</sub></sub>n<sup>k</sup>)
  208.          *  </pre>
  209.          *
  210.          *  where the P<sub>s<sub>j</sub></sub> and Q<sub>s<sub>k</sub></sub> are polynomials in s,
  211.          *  s being the index of the Y<sub>ρ,σ</sub><sup>n,s</sup> function
  212.          *
  213.          *  @param poly1 first list of polynomials
  214.          *  @param poly2 second list of polynomials
  215.          *  @return R<sub>s</sub>(n) as a list of {@link PolynomialFunction}
  216.          */
  217.         private static List<PolynomialFunction> multiplyPolynomialList(final List<PolynomialFunction> poly1,
  218.                                                                        final List<PolynomialFunction> poly2) {
  219.             // Initialize the result list of polynomial function
  220.             final List<PolynomialFunction> result = new ArrayList<PolynomialFunction>();
  221.             initializeListOfPolynomials(poly1.size() + poly2.size() - 1, result);

  222.             int i = 0;
  223.             // Iterate over first polynomial list
  224.             for (PolynomialFunction f1 : poly1) {
  225.                 // Iterate over second polynomial list
  226.                 for (int j = i; j < poly2.size() + i; j++) {
  227.                     final PolynomialFunction p2 = poly2.get(j - i);
  228.                     // Get previous polynomial for current 'n' order
  229.                     final PolynomialFunction previousP2 = result.get(j);
  230.                     // Replace the current order by summing the product of both of the polynomials
  231.                     result.set(j, previousP2.add(f1.multiply(p2)));
  232.                 }
  233.                 // shift polynomial order in 'n'
  234.                 i++;
  235.             }
  236.             return result;
  237.         }

  238.         /** Sum two lists of {@link PolynomialFunction}.
  239.          *  @param poly1 first list
  240.          *  @param poly2 second list
  241.          *  @return the summed list
  242.          */
  243.         private static List<PolynomialFunction> sumPolynomialList(final List<PolynomialFunction> poly1,
  244.                                                                   final List<PolynomialFunction> poly2) {
  245.             // identify the lowest degree polynomial
  246.             final int lowLength  = FastMath.min(poly1.size(), poly2.size());
  247.             final int highLength = FastMath.max(poly1.size(), poly2.size());
  248.             // Initialize the result list of polynomial function
  249.             final List<PolynomialFunction> result = new ArrayList<PolynomialFunction>();
  250.             initializeListOfPolynomials(highLength, result);

  251.             for (int i = 0; i < lowLength; i++) {
  252.                 // Add polynomials by increasing order of 'n'
  253.                 result.set(i, poly1.get(i).add(poly2.get(i)));
  254.             }
  255.             // Complete the list if lists are of different size:
  256.             for (int i = lowLength; i < highLength; i++) {
  257.                 if (poly1.size() < poly2.size()) {
  258.                     result.set(i, poly2.get(i));
  259.                 } else {
  260.                     result.set(i, poly1.get(i));
  261.                 }
  262.             }
  263.             return result;
  264.         }

  265.         /** Initialize an empty list of polynomials.
  266.          *  @param i order
  267.          *  @param result list into which polynomials should be added
  268.          */
  269.         private static void initializeListOfPolynomials(final int i,
  270.                                                         final List<PolynomialFunction> result) {
  271.             for (int k = 0; k < i; k++) {
  272.                 result.add(new PolynomialFunction(new double[] {0.}));
  273.             }
  274.         }

  275.         /** Shift a list of {@link PolynomialFunction}.
  276.          *
  277.          *  @param polynomialList list of {@link PolynomialFunction}
  278.          *  @param shift shift value
  279.          *  @return new list of shifted {@link PolynomialFunction}
  280.          */
  281.         private static List<PolynomialFunction> shiftList(final List<PolynomialFunction> polynomialList,
  282.                                                           final int shift) {
  283.             final List<PolynomialFunction> shiftedList = new ArrayList<PolynomialFunction>();
  284.             for (PolynomialFunction function : polynomialList) {
  285.                 shiftedList.add(new PolynomialFunction(shift(function.getCoefficients(), shift)));
  286.             }
  287.             return shiftedList;
  288.         }

  289.         /**
  290.          * Compute the coefficients of the polynomial \(P_s(x)\)
  291.          * whose values at point {@code x} will be the same as the those from the
  292.          * original polynomial \(P(x)\) when computed at {@code x + shift}.
  293.          * <p>
  294.          * More precisely, let \(\Delta = \) {@code shift} and let
  295.          * \(P_s(x) = P(x + \Delta)\).  The returned array
  296.          * consists of the coefficients of \(P_s\).  So if \(a_0, ..., a_{n-1}\)
  297.          * are the coefficients of \(P\), then the returned array
  298.          * \(b_0, ..., b_{n-1}\) satisfies the identity
  299.          * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
  300.          * </p>
  301.          * <p>
  302.          * This method is a modified version of the method with the same name
  303.          * in Hipparchus {@code PolynomialsUtils} class, simply changing
  304.          * computation of binomial coefficients so degrees higher than 66 can be used.
  305.          * </p>
  306.          *
  307.          * @param coefficients Coefficients of the original polynomial.
  308.          * @param shift Shift value.
  309.          * @return the coefficients \(b_i\) of the shifted
  310.          * polynomial.
  311.          */
  312.         public static double[] shift(final double[] coefficients,
  313.                                      final double shift) {
  314.             final int dp1 = coefficients.length;
  315.             final double[] newCoefficients = new double[dp1];

  316.             // Pascal triangle.
  317.             final double[][] coeff = new double[dp1][dp1];
  318.             coeff[0][0] = 1;
  319.             for (int i = 1; i < dp1; i++) {
  320.                 coeff[i][0] = 1;
  321.                 for (int j = 1; j < i; j++) {
  322.                     coeff[i][j] = coeff[i - 1][j - 1] + coeff[i - 1][j];
  323.                 }
  324.                 coeff[i][i] = 1;
  325.             }

  326.             // First polynomial coefficient.
  327.             double shiftI = 1;
  328.             for (int i = 0; i < dp1; i++) {
  329.                 newCoefficients[0] += coefficients[i] * shiftI;
  330.                 shiftI *= shift;
  331.             }

  332.             // Superior order.
  333.             final int d = dp1 - 1;
  334.             for (int i = 0; i < d; i++) {
  335.                 double shiftJmI = 1;
  336.                 for (int j = i; j < d; j++) {
  337.                     newCoefficients[i + 1] += coeff[j + 1][j - i] * coefficients[j + 1] * shiftJmI;
  338.                     shiftJmI *= shift;
  339.                 }
  340.             }

  341.             return newCoefficients;
  342.         }

  343.         /** Generate recurrence coefficients.
  344.          *
  345.          *  @param rho ρ value
  346.          *  @param sigma σ value
  347.          *  @return recurrence coefficients
  348.          */
  349.         private static Map<Integer, List<PolynomialFunction>> generateRecurrenceCoefficients(final int rho, final int sigma) {
  350.             final double den   = 1. / (4. * (rho + sigma));
  351.             final double denx2 = 2. * den;
  352.             final double denx4 = 4. * den;
  353.             // Initialization :
  354.             final Map<Integer, List<PolynomialFunction>> list = new TreeMap<Integer, List<PolynomialFunction>>();
  355.             final List<PolynomialFunction> poly0 = new ArrayList<PolynomialFunction>();
  356.             final List<PolynomialFunction> poly1 = new ArrayList<PolynomialFunction>();
  357.             final List<PolynomialFunction> poly2 = new ArrayList<PolynomialFunction>();
  358.             final List<PolynomialFunction> poly3 = new ArrayList<PolynomialFunction>();
  359.             final List<PolynomialFunction> poly4 = new ArrayList<PolynomialFunction>();
  360.             // (s - n)
  361.             poly0.add(new PolynomialFunction(new double[] {0., den}));
  362.             poly0.add(new PolynomialFunction(new double[] {-den}));
  363.             // 2(2 * rho + 2 * sigma + 2 + 3*n)
  364.             poly1.add(new PolynomialFunction(new double[] {1. + denx4}));
  365.             poly1.add(new PolynomialFunction(new double[] {denx2 + denx4}));
  366.             // 2(2s - n)
  367.             poly2.add(new PolynomialFunction(new double[] {0., denx4}));
  368.             poly2.add(new PolynomialFunction(new double[] {-denx2}));
  369.             // - (s + n)
  370.             poly3.add(new PolynomialFunction(new double[] {0., -den}));
  371.             poly3.add(new PolynomialFunction(new double[] {-den}));
  372.             // - 2(2s + n)
  373.             poly4.add(new PolynomialFunction(new double[] {0., -denx4}));
  374.             poly4.add(new PolynomialFunction(new double[] {-denx2}));

  375.             // Fill the map :
  376.             list.put(0, poly0);
  377.             list.put(1, poly1);
  378.             list.put(2, poly2);
  379.             list.put(3, poly3);
  380.             list.put(4, poly4);
  381.             return list;
  382.         }

  383.     }

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

  386.         /** first couple value. */
  387.         private final int rho;

  388.         /** second couple value. */
  389.         private final int sigma;

  390.         /** Constructor.
  391.          * @param rho first couple value
  392.          * @param sigma second couple value
  393.          */
  394.         private Couple(final int rho, final int sigma) {
  395.             this.rho = rho;
  396.             this.sigma = sigma;
  397.         }

  398.         /** {@inheritDoc} */
  399.         public int compareTo(final Couple c) {
  400.             int result = 1;
  401.             if (rho == c.rho) {
  402.                 if (sigma < c.sigma) {
  403.                     result = -1;
  404.                 } else if (sigma == c.sigma) {
  405.                     result = 0;
  406.                 }
  407.             } else if (rho < c.rho) {
  408.                 result = -1;
  409.             }
  410.             return result;
  411.         }

  412.         /** {@inheritDoc} */
  413.         public boolean equals(final Object couple) {

  414.             if (couple == this) {
  415.                 // first fast check
  416.                 return true;
  417.             }

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

  421.             return false;

  422.         }

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

  427.     }

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

  430.         /** n value. */
  431.         private final int n;

  432.         /** s value. */
  433.         private final int s;

  434.         /** ρ value. */
  435.         private final int rho;

  436.         /** σ value. */
  437.         private final int sigma;

  438.         /** Simpleconstructor.
  439.          * @param n n
  440.          * @param s s
  441.          * @param rho ρ
  442.          * @param sigma σ
  443.          */
  444.         NewKey(final int n, final int s, final int rho, final int sigma) {
  445.             this.n = n;
  446.             this.s = s;
  447.             this.rho = rho;
  448.             this.sigma = sigma;
  449.         }

  450.         /** {@inheritDoc} */
  451.         public int compareTo(final NewKey key) {
  452.             int result = 1;
  453.             if (n == key.n) {
  454.                 if (s == key.s) {
  455.                     if (rho == key.rho) {
  456.                         if (sigma < key.sigma) {
  457.                             result = -1;
  458.                         } else if (sigma == key.sigma) {
  459.                             result = 0;
  460.                         } else {
  461.                             result = 1;
  462.                         }
  463.                     } else if (rho < key.rho) {
  464.                         result = -1;
  465.                     } else {
  466.                         result = 1;
  467.                     }
  468.                 } else if (s < key.s) {
  469.                     result = -1;
  470.                 } else {
  471.                     result = 1;
  472.                 }
  473.             } else if (n < key.n) {
  474.                 result = -1;
  475.             }
  476.             return result;
  477.         }

  478.         /** {@inheritDoc} */
  479.         public boolean equals(final Object key) {

  480.             if (key == this) {
  481.                 // first fast check
  482.                 return true;
  483.             }

  484.             if ((key != null) && (key instanceof NewKey)) {
  485.                 return (n     == ((NewKey) key).n) &&
  486.                        (s     == ((NewKey) key).s) &&
  487.                        (rho   == ((NewKey) key).rho) &&
  488.                        (sigma == ((NewKey) key).sigma);
  489.             }

  490.             return false;

  491.         }

  492.         /** {@inheritDoc} */
  493.         public int hashCode() {
  494.             return 0x25baa451 ^ (n << 24) ^ (s << 16) ^ (rho << 8) ^ sigma;
  495.         }

  496.     }

  497. }