ZonalContribution.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.forces;

  19. import java.util.TreeMap;

  21. import org.apache.commons.math3.util.FastMath;
  22. import org.orekit.errors.OrekitException;
  23. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  24. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  25. import org.orekit.propagation.SpacecraftState;
  26. import org.orekit.propagation.events.EventDetector;
  27. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  28. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  29. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
  30. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
  31. import org.orekit.time.AbsoluteDate;

  33. /** Zonal contribution to the {@link DSSTCentralBody central body gravitational perturbation}.
  34.  *
  35.  *   @author Romain Di Costanzo
  36.  *   @author Pascal Parraud
  37.  */
  38. class ZonalContribution implements DSSTForceModel {

  40.     /** Truncation tolerance. */
  41.     private static final double TRUNCATION_TOLERANCE = 1e-4;

  43.     /** Provider for spherical harmonics. */
  44.     private final UnnormalizedSphericalHarmonicsProvider provider;

  46.     /** Maximal degree to consider for harmonics potential. */
  47.     private final int maxDegree;

  49.     /** Maximal degree to consider for harmonics potential. */
  50.     private final int maxOrder;

  52.     /** Factorial. */
  53.     private final double[] fact;

  55.     /** Coefficient used to define the mean disturbing function V<sub>ns</sub> coefficient. */
  56.     private final TreeMap<NSKey, Double> Vns;

  58.     /** Highest power of the eccentricity to be used in series expansion. */
  59.     private int maxEccPow;

  61.     // Equinoctial elements (according to DSST notation)
  62.     /** a. */
  63.     private double a;
  64.     /** ex. */
  65.     private double k;
  66.     /** ey. */
  67.     private double h;
  68.     /** hx. */
  69.     private double q;
  70.     /** hy. */
  71.     private double p;
  72.     /** Retrograde factor. */
  73.     private int    I;

  75.     /** Eccentricity. */
  76.     private double ecc;

  78.     /** Direction cosine &alpha. */
  79.     private double alpha;
  80.     /** Direction cosine &beta. */
  81.     private double beta;
  82.     /** Direction cosine &gamma. */
  83.     private double gamma;

  85.     // Common factors for potential computation
  86.     /** &Chi; = 1 / sqrt(1 - e<sup>2</sup>) = 1 / B. */
  87.     private double X;
  88.     /** &Chi;<sup>2</sup>. */
  89.     private double XX;
  90.     /** &Chi;<sup>3</sup>. */
  91.     private double XXX;
  92.     /** 1 / (A * B) .*/
  93.     private double ooAB;
  94.     /** B / A .*/
  95.     private double BoA;
  96.     /** B / A(1 + B) .*/
  97.     private double BoABpo;
  98.     /** -C / (2 * A * B) .*/
  99.     private double mCo2AB;
  100.     /** -2 * a / A .*/
  101.     private double m2aoA;
  102.     /** &mu; / a .*/
  103.     private double muoa;

  105.     /** Simple constructor.
  106.      * @param provider provider for spherical harmonics
  107.      */
  108.     public ZonalContribution(final UnnormalizedSphericalHarmonicsProvider provider) {

  110.         this.provider  = provider;
  111.         this.maxDegree = provider.getMaxDegree();
  112.         this.maxOrder  = provider.getMaxOrder();

  114.         // Vns coefficients
  115.         this.Vns = CoefficientsFactory.computeVns(maxDegree + 1);

  117.         // Factorials computation
  118.         final int maxFact = 2 * maxDegree + 1;
  119.         this.fact = new double[maxFact];
  120.         fact[0] = 1.;
  121.         for (int i = 1; i < maxFact; i++) {
  122.             fact[i] = i * fact[i - 1];
  123.         }

  125.         // Initialize default values
  126.         this.maxEccPow = (maxDegree == 2) ? 0 : Integer.MIN_VALUE;
  127.     }

  129.     /** Get the spherical harmonics provider.
  130.      *  @return the spherical harmonics provider
  131.      */
  132.     public UnnormalizedSphericalHarmonicsProvider getProvider() {
  133.         return provider;
  134.     }

  136.     /** Computes the highest power of the eccentricity to appear in the truncated
  137.      *  analytical power series expansion.
  138.      *  <p>
  139.      *  This method computes the upper value for the central body potential and
  140.      *  determines the maximal power for the eccentricity producing potential
  141.      *  terms bigger than a defined tolerance.
  142.      *  </p>
  143.      *  @param aux auxiliary elements related to the current orbit
  144.      *  @throws OrekitException if some specific error occurs
  145.      */
  146.     public void initialize(final AuxiliaryElements aux)
  147.         throws OrekitException {

  149.         if (maxDegree == 2) {
  150.             maxEccPow = 0;
  151.         } else {
  152.             // Initializes specific parameters.
  153.             initializeStep(aux);
  154.             final UnnormalizedSphericalHarmonics harmonics = provider.onDate(aux.getDate());

  156.             // Utilities for truncation
  157.             final double ax2or = 2. * a / provider.getAe();
  158.             double xmuran = provider.getMu() / a;
  159.             // Set a lower bound for eccentricity
  160.             final double eo2  = FastMath.max(0.0025, ecc / 2.);
  161.             final double x2o2 = XX / 2.;
  162.             final double[] eccPwr = new double[maxDegree + 1];
  163.             final double[] chiPwr = new double[maxDegree + 1];
  164.             final double[] hafPwr = new double[maxDegree + 1];
  165.             eccPwr[0] = 1.;
  166.             chiPwr[0] = X;
  167.             hafPwr[0] = 1.;
  168.             for (int i = 1; i <= maxDegree; i++) {
  169.                 eccPwr[i] = eccPwr[i - 1] * eo2;
  170.                 chiPwr[i] = chiPwr[i - 1] * x2o2;
  171.                 hafPwr[i] = hafPwr[i - 1] * 0.5;
  172.                 xmuran  /= ax2or;
  173.             }

  175.             // Set highest power of e and degree of current spherical harmonic.
  176.             maxEccPow = 0;
  177.             int n = maxDegree;
  178.             // Loop over n
  179.             do {
  180.                 // Set order of current spherical harmonic.
  181.                 int m = 0;
  182.                 // Loop over m
  183.                 do {
  184.                     // Compute magnitude of current spherical harmonic coefficient.
  185.                     final double cnm = harmonics.getUnnormalizedCnm(n, m);
  186.                     final double snm = harmonics.getUnnormalizedSnm(n, m);
  187.                     final double csnm = FastMath.hypot(cnm, snm);
  188.                     if (csnm == 0.) break;
  189.                     // Set magnitude of last spherical harmonic term.
  190.                     double lastTerm = 0.;
  191.                     // Set current power of e and related indices.
  192.                     int nsld2 = (n - maxEccPow - 1) / 2;
  193.                     int l = n - 2 * nsld2;
  194.                     // Loop over l
  195.                     double term = 0.;
  196.                     do {
  197.                         // Compute magnitude of current spherical harmonic term.
  198.                         if (m < l) {
  199.                             term = csnm * xmuran *
  200.                                     (fact[n - l] / (fact[n - m])) *
  201.                                     (fact[n + l] / (fact[nsld2] * fact[nsld2 + l])) *
  202.                                     eccPwr[l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) *
  203.                                     (UpperBounds.getRnml(gamma, n, l, m, 1, I) + UpperBounds.getRnml(gamma, n, l, m, -1, I));
  204.                         } else {
  205.                             term = csnm * xmuran *
  206.                                     (fact[n + m] / (fact[nsld2] * fact[nsld2 + l])) *
  207.                                     eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) *
  208.                                     (UpperBounds.getRnml(gamma, n, m, l, 1, I) + UpperBounds.getRnml(gamma, n, m, l, -1, I));
  209.                         }
  210.                         // Is the current spherical harmonic term bigger than the truncation tolerance ?
  211.                         if (term >= TRUNCATION_TOLERANCE) {
  212.                             maxEccPow = l;
  213.                         } else {
  214.                             // Is the current term smaller than the last term ?
  215.                             if (term < lastTerm) {
  216.                                 break;
  217.                             }
  218.                         }
  219.                         // Proceed to next power of e.
  220.                         lastTerm = term;
  221.                         l += 2;
  222.                         nsld2--;
  223.                     } while (l < n);
  224.                     // Is the current spherical harmonic term bigger than the truncation tolerance ?
  225.                     if (term >= TRUNCATION_TOLERANCE) {
  226.                         maxEccPow = FastMath.min(maxDegree - 2, maxEccPow);
  227.                         return;
  228.                     }
  229.                     // Proceed to next order.
  230.                     m++;
  231.                 } while (m <= FastMath.min(n, maxOrder));
  232.                 // Procced to next degree.
  233.                 xmuran *= ax2or;
  234.                 n--;
  235.             } while (n > maxEccPow + 2);

  237.             maxEccPow = FastMath.min(maxDegree - 2, maxEccPow);
  238.         }
  239.     }

  241.     /** {@inheritDoc} */
  242.     public void initializeStep(final AuxiliaryElements aux) throws OrekitException {

  244.         // Equinoctial elements
  245.         a = aux.getSma();
  246.         k = aux.getK();
  247.         h = aux.getH();
  248.         q = aux.getQ();
  249.         p = aux.getP();

  251.         // Retrograde factor
  252.         I = aux.getRetrogradeFactor();

  254.         // Eccentricity
  255.         ecc = aux.getEcc();

  257.         // Direction cosines
  258.         alpha = aux.getAlpha();
  259.         beta  = aux.getBeta();
  260.         gamma = aux.getGamma();

  262.         // Equinoctial coefficients
  263.         final double AA = aux.getA();
  264.         final double BB = aux.getB();
  265.         final double CC = aux.getC();

  267.         // &Chi; = 1 / B
  268.         X = 1. / BB;
  269.         XX = X * X;
  270.         XXX = X * XX;
  271.         // 1 / AB
  272.         ooAB = 1. / (AA * BB);
  273.         // B / A
  274.         BoA = BB / AA;
  275.         // -C / 2AB
  276.         mCo2AB = -CC * ooAB / 2.;
  277.         // B / A(1 + B)
  278.         BoABpo = BoA / (1. + BB);
  279.         // -2 * a / A
  280.         m2aoA = -2 * a / AA;
  281.         // &mu; / a
  282.         muoa = provider.getMu() / a;

  284.     }

  286.     /** {@inheritDoc} */
  287.     public double[] getMeanElementRate(final SpacecraftState spacecraftState) throws OrekitException {

  289.         // Compute potential derivative
  290.         final double[] dU  = computeUDerivatives(spacecraftState.getDate());
  291.         final double dUda  = dU[0];
  292.         final double dUdk  = dU[1];
  293.         final double dUdh  = dU[2];
  294.         final double dUdAl = dU[3];
  295.         final double dUdBe = dU[4];
  296.         final double dUdGa = dU[5];

  298.         // Compute cross derivatives [Eq. 2.2-(8)]
  299.         // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
  300.         final double UAlphaGamma   = alpha * dUdGa - gamma * dUdAl;
  301.         // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
  302.         final double UBetaGamma    =  beta * dUdGa - gamma * dUdBe;
  303.         // Common factor
  304.         final double pUAGmIqUBGoAB = (p * UAlphaGamma - I * q * UBetaGamma) * ooAB;

  306.         // Compute mean elements rates [Eq. 3.1-(1)]
  307.         final double da =  0.;
  308.         final double dh =  BoA * dUdk + k * pUAGmIqUBGoAB;
  309.         final double dk = -BoA * dUdh - h * pUAGmIqUBGoAB;
  310.         final double dp =  mCo2AB * UBetaGamma;
  311.         final double dq =  mCo2AB * UAlphaGamma * I;
  312.         final double dM =  m2aoA * dUda + BoABpo * (h * dUdh + k * dUdk) + pUAGmIqUBGoAB;

  314.         return new double[] {da, dk, dh, dq, dp, dM};
  315.     }

  317.     /** {@inheritDoc} */
  318.     public double[] getShortPeriodicVariations(final AbsoluteDate date,
  319.                                                final double[] meanElements)
  320.         throws OrekitException {
  321.         // TODO: not implemented yet, Short Periodic Variations are set to null
  322.         return new double[] {0., 0., 0., 0., 0., 0.};
  323.     }

  325.     /** {@inheritDoc} */
  326.     public EventDetector[] getEventsDetectors() {
  327.         return null;
  328.     }

  330.     /** Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
  331.      *  <p>
  332.      *  The result is the array
  333.      *  [dU/da, dU/dk, dU/dh, dU/d&alpha;, dU/d&beta;, dU/d&gamma;]
  334.      *  </p>
  335.      *  @param date current date
  336.      *  @return potential derivatives
  337.      *  @throws OrekitException if an error occurs in hansen computation
  338.      */
  339.     private double[] computeUDerivatives(final AbsoluteDate date) throws OrekitException {

  341.         final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);

  343.         // Hansen coefficients
  344.         final HansenZonal hansen = new HansenZonal();
  345.         // Gs and Hs coefficients
  346.         final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPow);
  347.         // Qns coefficients
  348.         final double[][] Qns  = CoefficientsFactory.computeQns(gamma, maxDegree, maxEccPow);
  349.         // r / a up to power degree
  350.         final double roa = provider.getAe() / a;

  352.         final double[] roaPow = new double[maxDegree + 1];
  353.         roaPow[0] = 1.;
  354.         for (int i = 1; i <= maxDegree; i++) {
  355.             roaPow[i] = roa * roaPow[i - 1];
  356.         }

  358.         // Potential derivatives
  359.         double dUda  = 0.;
  360.         double dUdk  = 0.;
  361.         double dUdh  = 0.;
  362.         double dUdAl = 0.;
  363.         double dUdBe = 0.;
  364.         double dUdGa = 0.;

  366.         for (int s = 0; s <= maxEccPow; s++) {
  367.             // Get the current Gs coefficient
  368.             final double gs = GsHs[0][s];

  370.             // Compute Gs partial derivatives from 3.1-(9)
  371.             double dGsdh  = 0.;
  372.             double dGsdk  = 0.;
  373.             double dGsdAl = 0.;
  374.             double dGsdBe = 0.;
  375.             if (s > 0) {
  376.                 // First get the G(s-1) and the H(s-1) coefficients
  377.                 final double sxgsm1 = s * GsHs[0][s - 1];
  378.                 final double sxhsm1 = s * GsHs[1][s - 1];
  379.                 // Then compute derivatives
  380.                 dGsdh  = beta  * sxgsm1 - alpha * sxhsm1;
  381.                 dGsdk  = alpha * sxgsm1 + beta  * sxhsm1;
  382.                 dGsdAl = k * sxgsm1 - h * sxhsm1;
  383.                 dGsdBe = h * sxgsm1 + k * sxhsm1;
  384.             }

  386.             // Kronecker symbol (2 - delta(0,s))
  387.             final double d0s = (s == 0) ? 1 : 2;

  389.             for (int n = s + 2; n <= maxDegree; n++) {
  390.                 // (n - s) must be even
  391.                 if ((n - s) % 2 == 0) {
  392.                     // Extract data from previous computation :
  393.                     final double kns   = hansen.getValue(-n - 1, s);
  394.                     final double dkns  = hansen.getDerivative(-n - 1, s);
  395.                     final double vns   = Vns.get(new NSKey(n, s));
  396.                     final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
  397.                     final double coef1 = coef0 * Qns[n][s];
  398.                     final double coef2 = coef1 * kns;
  399.                     // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
  400.                     final double dqns  = Qns[n][s + 1];

  402.                     // Compute dU / da :
  403.                     dUda += coef2 * (n + 1) * gs;
  404.                     // Compute dU / dEx
  405.                     dUdk += coef1 * (kns * dGsdk + k * XXX * gs * dkns);
  406.                     // Compute dU / dEy
  407.                     dUdh += coef1 * (kns * dGsdh + h * XXX * gs * dkns);
  408.                     // Compute dU / dAlpha
  409.                     dUdAl += coef2 * dGsdAl;
  410.                     // Compute dU / dBeta
  411.                     dUdBe += coef2 * dGsdBe;
  412.                     // Compute dU / dGamma
  413.                     dUdGa += coef0 * kns * dqns * gs;
  414.                 }
  415.             }
  416.         }

  418.         return new double[] {
  419.             dUda  *  muoa / a,
  420.             dUdk  * -muoa,
  421.             dUdh  * -muoa,
  422.             dUdAl * -muoa,
  423.             dUdBe * -muoa,
  424.             dUdGa * -muoa
  425.         };
  426.     }

  428.     /** Hansen coefficients for zonal contribution to central body force model.
  429.      *  <p>
  430.      *  Hansen coefficients are functions of the eccentricity.
  431.      *  For a given eccentricity, the computed elements are stored in a map.
  432.      *  </p>
  433.      */
  434.     private class HansenZonal {

  436.         /** Map to store every Hansen coefficient computed. */
  437.         private TreeMap<NSKey, Double> coefficients;

  439.         /** Map to store every Hansen coefficient derivative computed. */
  440.         private TreeMap<NSKey, Double> derivatives;

  442.         /** Simple constructor. */
  443.         public HansenZonal() {
  444.             coefficients = new TreeMap<CoefficientsFactory.NSKey, Double>();
  445.             derivatives  = new TreeMap<CoefficientsFactory.NSKey, Double>();
  446.             initialize();
  447.         }

  449.         /** Get the K<sub>0</sub><sup>-n-1,s</sup> coefficient value.
  450.          *  @param mnm1 -n-1 value
  451.          *  @param s s value
  452.          *  @return K<sub>0</sub><sup>-n-1,s</sup>
  453.          */
  454.         public double getValue(final int mnm1, final int s) {
  455.             if (coefficients.containsKey(new NSKey(mnm1, s))) {
  456.                 return coefficients.get(new NSKey(mnm1, s));
  457.             } else {
  458.                 // Compute the K0(-n-1,s) coefficients
  459.                 return computeValue(-mnm1 - 1, s);
  460.             }
  461.         }

  463.         /** Get the dK<sub>0</sub><sup>-n-1,s</sup> / d&chi; coefficient derivative.
  464.          *  @param mnm1 -n-1 value
  465.          *  @param s s value
  466.          *  @return dK<sub>0</sub><sup>-n-1,s</sup> / d&chi;
  467.          */
  468.         public double getDerivative(final int mnm1, final int s) {
  469.             if (derivatives.containsKey(new NSKey(mnm1, s))) {
  470.                 return derivatives.get(new NSKey(mnm1, s));
  471.             } else {
  472.                 // Compute the dK0(-n-1,s) / dX derivative
  473.                 return computeDerivative(-mnm1 - 1, s);
  474.             }
  475.         }

  477.         /** Kernels initialization. */
  478.         private void initialize() {
  479.             // coefficients
  480. //            coefficients.put(new NSKey(0, 0), 1.);
  481.             coefficients.put(new NSKey(-1, 0), 0.);
  482.             coefficients.put(new NSKey(-1, 1), 0.);
  483.             coefficients.put(new NSKey(-2, 0), X);
  484.             coefficients.put(new NSKey(-2, 1), 0.);
  485.             coefficients.put(new NSKey(-3, 0), XXX);
  486.             coefficients.put(new NSKey(-3, 1), 0.5 * XXX);
  487.             // derivatives
  488.             derivatives.put(new NSKey(-1, 0), 0.);
  489.             derivatives.put(new NSKey(-2, 0), 1.);
  490.             derivatives.put(new NSKey(-2, 1), 0.);
  491.             derivatives.put(new NSKey(-3, 1), 1.5 * XX);
  492.         }

  494.         /** Compute the K<sub>0</sub><sup>-n-1,s</sup> coefficient from equation 2.7.3-(6).
  495.          * @param n n  positive value. For a given 'n', the K<sub>0</sub><sup>-n-1,s</sup> is returned.
  496.          * @param s s value
  497.          * @return K<sub>0</sub><sup>-n-1,s</sup>
  498.          */
  499.         private double computeValue(final int n, final int s) {
  500.             // Initialize return value
  501.             double kns = 0.;

  503.             final NSKey key = new NSKey(-n - 1, s);

  505.             if (coefficients.containsKey(key)) {
  506.                 kns = coefficients.get(key);
  507.             } else {
  508.                 if (n == s) {
  509.                     kns = 0.;
  510.                 } else if (n == (s + 1)) {
  511.                     kns = FastMath.pow(X, 1 + 2 * s) / FastMath.pow(2, s);
  512.                 } else {
  513.                     final NSKey keymNS = new NSKey(-n, s);
  514.                     double kmNS = 0.;
  515.                     if (coefficients.containsKey(keymNS)) {
  516.                         kmNS = coefficients.get(keymNS);
  517.                     } else {
  518.                         kmNS = computeValue(n - 1, s);
  519.                     }

  521.                     final NSKey keymNp1S = new NSKey(-(n - 1), s);
  522.                     double kmNp1S = 0.;
  523.                     if (coefficients.containsKey(keymNp1S)) {
  524.                         kmNp1S = coefficients.get(keymNp1S);
  525.                     } else {
  526.                         kmNp1S = computeValue(n - 2, s);
  527.                     }

  529.                     kns = (n - 1.) * XX * ((2. * n - 3.) * kmNS - (n - 2.) * kmNp1S);
  530.                     kns /= (n + s - 1.) * (n - s - 1.);
  531.                 }
  532.                 // Add K(-n-1, s)
  533.                 coefficients.put(key, kns);
  534.             }
  535.             return kns;
  536.         }

  538.         /** Compute dK<sub>0</sub><sup>-n-1,s</sup> / d&chi; from Equation 3.1-(7).
  539.          *  @param n n positive value. For a given 'n', the dK<sub>0</sub><sup>-n-1,s</sup> / d&chi; is returned.
  540.          *  @param s s value
  541.          *  @return dK<sub>0</sub><sup>-n-1,s</sup> / d&chi;
  542.          */
  543.         private double computeDerivative(final int n, final int s) {
  544.             // Initialize return value
  545.             double dkdxns = 0.;

  547.             final NSKey key = new NSKey(-n - 1, s);
  548.             if (n == s) {
  549.                 derivatives.put(key, 0.);
  550.             } else if (n == s + 1) {
  551.                 dkdxns = (1. + 2. * s) * FastMath.pow(X, 2 * s) / FastMath.pow(2, s);
  552.             } else {
  553.                 final NSKey keymN = new NSKey(-n, s);
  554.                 double dkmN = 0.;
  555.                 if (derivatives.containsKey(keymN)) {
  556.                     dkmN = derivatives.get(keymN);
  557.                 } else {
  558.                     dkmN = computeDerivative(n - 1, s);
  559.                 }
  560.                 final NSKey keymNp1 = new NSKey(-n + 1, s);
  561.                 double dkmNp1 = 0.;
  562.                 if (derivatives.containsKey(keymNp1)) {
  563.                     dkmNp1 = derivatives.get(keymNp1);
  564.                 } else {
  565.                     dkmNp1 = computeDerivative(n - 2, s);
  566.                 }
  567.                 final double kns = getValue(-n - 1, s);
  568.                 dkdxns = (n - 1) * XX * ((2. * n - 3.) * dkmN - (n - 2.) * dkmNp1) / ((n + s - 1.) * (n - s - 1.));
  569.                 dkdxns += 2. * kns / X;
  570.             }

  572.             derivatives.put(key, dkdxns);
  573.             return dkdxns;
  574.         }

  576.     }

  578. }
Created with JaCoCo 0.6.3.201306030806