1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF 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  
18  package org.apache.commons.math.ode;
19  
20  /**
21   * This class is used in the junit tests for the ODE integrators.
22  
23   * <p>This specific problem is the following differential equation :
24   * <pre>
25   *    y1'' = -y1/r^3  y1 (0) = 1-e  y1' (0) = 0
26   *    y2'' = -y2/r^3  y2 (0) = 0    y2' (0) =sqrt((1+e)/(1-e))
27   *    r = sqrt (y1^2 + y2^2), e = 0.9
28   * </pre>
29   * This is a two-body problem in the plane which can be solved by
30   * Kepler's equation
31   * <pre>
32   *   y1 (t) = ...
33   * </pre>
34   * </p>
35  
36   */
37  public class TestProblem3
38    extends TestProblemAbstract {
39  
40    /** Serializable version identifier. */
41    private static final long serialVersionUID = 8567328542728919999L;
42  
43    /** Eccentricity */
44    double e;
45  
46    /** theoretical state */
47    private double[] y;
48  
49    /**
50     * Simple constructor.
51     * @param e eccentricity
52     */
53    public TestProblem3(double e) {
54      super();
55      this.e = e;
56      double[] y0 = { 1 - e, 0, 0, Math.sqrt((1+e)/(1-e)) };
57      setInitialConditions(0.0, y0);
58      setFinalConditions(20.0);
59      double[] errorScale = { 1.0, 1.0, 1.0, 1.0 };
60      setErrorScale(errorScale);
61      y = new double[y0.length];
62    }
63   
64    /**
65     * Simple constructor.
66     */
67    public TestProblem3() {
68      this(0.1);
69    }
70   
71    /**
72     * Copy constructor.
73     * @param problem problem to copy
74     */
75    public TestProblem3(TestProblem3 problem) {
76      super(problem);
77      e = problem.e;
78      y = problem.y.clone();
79    }
80  
81    /** {@inheritDoc} */
82    public TestProblem3 copy() {
83      return new TestProblem3(this);
84    }
85  
86    @Override
87    public void doComputeDerivatives(double t, double[] y, double[] yDot) {
88  
89      // current radius
90      double r2 = y[0] * y[0] + y[1] * y[1];
91      double invR3 = 1 / (r2 * Math.sqrt(r2));
92  
93      // compute the derivatives
94      yDot[0] = y[2];
95      yDot[1] = y[3];
96      yDot[2] = -invR3  * y[0];
97      yDot[3] = -invR3  * y[1];
98  
99    }
100 
101   @Override
102   public double[] computeTheoreticalState(double t) {
103 
104     // solve Kepler's equation
105     double E = t;
106     double d = 0;
107     double corr = 999.0;
108     for (int i = 0; (i < 50) && (Math.abs(corr) > 1.0e-12); ++i) {
109       double f2  = e * Math.sin(E);
110       double f0  = d - f2;
111       double f1  = 1 - e * Math.cos(E);
112       double f12 = f1 + f1;
113       corr  = f0 * f12 / (f1 * f12 - f0 * f2);
114       d -= corr;
115       E = t + d;
116     };
117 
118     double cosE = Math.cos(E);
119     double sinE = Math.sin(E);
120 
121     y[0] = cosE - e;
122     y[1] = Math.sqrt(1 - e * e) * sinE;
123     y[2] = -sinE / (1 - e * cosE);
124     y[3] = Math.sqrt(1 - e * e) * cosE / (1 - e * cosE);
125 
126     return y;
127   }
128 
129 }