001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.ode;
019    
020    /**
021     * This class is used in the junit tests for the ODE integrators.
022    
023     * <p>This specific problem is the following differential equation :
024     * <pre>
025     *    y1'' = -y1/r^3  y1 (0) = 1-e  y1' (0) = 0
026     *    y2'' = -y2/r^3  y2 (0) = 0    y2' (0) =sqrt((1+e)/(1-e))
027     *    r = sqrt (y1^2 + y2^2), e = 0.9
028     * </pre>
029     * This is a two-body problem in the plane which can be solved by
030     * Kepler's equation
031     * <pre>
032     *   y1 (t) = ...
033     * </pre>
034     * </p>
035    
036     */
037    public class TestProblem3
038      extends TestProblemAbstract {
039    
040      /** Serializable version identifier. */
041      private static final long serialVersionUID = 8567328542728919999L;
042    
043      /** Eccentricity */
044      double e;
045    
046      /** theoretical state */
047      private double[] y;
048    
049      /**
050       * Simple constructor.
051       * @param e eccentricity
052       */
053      public TestProblem3(double e) {
054        super();
055        this.e = e;
056        double[] y0 = { 1 - e, 0, 0, Math.sqrt((1+e)/(1-e)) };
057        setInitialConditions(0.0, y0);
058        setFinalConditions(20.0);
059        double[] errorScale = { 1.0, 1.0, 1.0, 1.0 };
060        setErrorScale(errorScale);
061        y = new double[y0.length];
062      }
063     
064      /**
065       * Simple constructor.
066       */
067      public TestProblem3() {
068        this(0.1);
069      }
070     
071      /**
072       * Copy constructor.
073       * @param problem problem to copy
074       */
075      public TestProblem3(TestProblem3 problem) {
076        super(problem);
077        e = problem.e;
078        y = problem.y.clone();
079      }
080    
081      /** {@inheritDoc} */
082      public TestProblem3 copy() {
083        return new TestProblem3(this);
084      }
085    
086      @Override
087      public void doComputeDerivatives(double t, double[] y, double[] yDot) {
088    
089        // current radius
090        double r2 = y[0] * y[0] + y[1] * y[1];
091        double invR3 = 1 / (r2 * Math.sqrt(r2));
092    
093        // compute the derivatives
094        yDot[0] = y[2];
095        yDot[1] = y[3];
096        yDot[2] = -invR3  * y[0];
097        yDot[3] = -invR3  * y[1];
098    
099      }
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    }