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  class TestProblem3
38    extends TestProblemAbstract {
39  
40    /** Eccentricity */
41    double e;
42  
43    /** theoretical state */
44    private double[] y;
45  
46    /**
47     * Simple constructor.
48     * @param e eccentricity
49     */
50    public TestProblem3(double e) {
51      super();
52      this.e = e;
53      double[] y0 = { 1 - e, 0, 0, Math.sqrt((1+e)/(1-e)) };
54      setInitialConditions(0.0, y0);
55      setFinalConditions(20.0);
56      double[] errorScale = { 1.0, 1.0, 1.0, 1.0 };
57      setErrorScale(errorScale);
58      y = new double[y0.length];
59    }
60   
61    /**
62     * Simple constructor.
63     */
64    public TestProblem3() {
65      this(0.1);
66    }
67   
68    /**
69     * Copy constructor.
70     * @param problem problem to copy
71     */
72    public TestProblem3(TestProblem3 problem) {
73      super(problem);
74      e = problem.e;
75      y = (double[]) problem.y.clone();
76    }
77  
78    /**
79     * Clone operation.
80     * @return a copy of the instance
81     */
82    public Object clone() {
83      return new TestProblem3(this);
84    }
85  
86    public void doComputeDerivatives(double t, double[] y, double[] yDot) {
87  
88      // current radius
89      double r2 = y[0] * y[0] + y[1] * y[1];
90      double invR3 = 1 / (r2 * Math.sqrt(r2));
91  
92      // compute the derivatives
93      yDot[0] = y[2];
94      yDot[1] = y[3];
95      yDot[2] = -invR3  * y[0];
96      yDot[3] = -invR3  * y[1];
97  
98    }
99  
100   public double[] computeTheoreticalState(double t) {
101 
102     // solve Kepler's equation
103     double E = t;
104     double d = 0;
105     double corr = 0;
106     do {
107       double f2  = e * Math.sin(E);
108       double f0  = d - f2;
109       double f1  = 1 - e * Math.cos(E);
110       double f12 = f1 + f1;
111       corr  = f0 * f12 / (f1 * f12 - f0 * f2);
112       d -= corr;
113       E = t + d;
114     } while (Math.abs(corr) > 1.0e-12);
115 
116     double cosE = Math.cos(E);
117     double sinE = Math.sin(E);
118 
119     y[0] = cosE - e;
120     y[1] = Math.sqrt(1 - e * e) * sinE;
121     y[2] = -sinE / (1 - e * cosE);
122     y[3] = Math.sqrt(1 - e * e) * cosE / (1 - e * cosE);
123 
124     return y;
125   }
126 
127 }