gistfile1.txt

do
            {
                steps++;
                activeStep = steps;
                activeDt = dt;
                repeatWithSmallerStep = false;
                Vector3d error;
                Vector3d zv;
                // Perform the RK4 calculation
                {
                    Vector3d dv1 = dt * TotalAccel(x, v, true);
                    Vector3d dx1 = dt * v;

                    Vector3d dv2 = dt * TotalAccel(x + 0.5 * dx1, v + 0.5 * dv1);
                    Vector3d dx2 = dt * (v + 0.5 * dv1);

                    Vector3d dv3 = dt * TotalAccel(x + 0.75 * dx2, v + 0.75 * dv2);
                    Vector3d dx3 = dt * (v + 0.75 * dv2);

                    Vector3d dv4 = dt * TotalAccel(x + 2d/9 * dx1 + 1d/3 * dx2 + 4d/9 * dx3, v + 2d / 9 * dv1 + 1d / 3 * dv2 + 4d / 9 * dv3);
                    //Vector3d dx4 = dt * (v + 2d / 9 * dv1 + 1d / 3 * dv2 + 4d / 9 * dv3);

                    dx = (2 * dx1 + 3 * dx2 + 4 * dx3 ) / 9.0;
                    dv = (2 * dv1 + 3 * dv2 + 4 * dv3 ) / 9.0;

                    //Vector3d zx = (7 * dx1 + 6 * dx2 + 8 * dx3 + 3 * dx4) / 24.0;
                    zv = (7 * dv1 + 6 * dv2 + 8 * dv3 + 3 * dv4) / 24.0;
                    error = zv - dv;
                }

                dt = Math.Max(dt * Math.Pow(tol / error.magnitude, 1 / (zv.magnitude + 1)), min_dt);

                // If the change in velocity is more than half the current velocity, then we need to try again with a smaller delta-t
                // or if dt is already small enough then continue anyway.
                if (error.magnitude > tol && dt > min_dt && dt < 100)
                {
                    repeatWithSmallerStep = true;
                }
                else
                {
                    // Consider opening the parachutes. If we do open them, and the dt is not as small as it could me, make it smaller and repeat,
                    Vector3 xForChuteSim = x + dx;
                    double altASL = xForChuteSim.magnitude - bodyRadius;
                    double altAGL = altASL - probableLandingSiteASL;
                    double pressure = Pressure(xForChuteSim);

                    bool willChutesOpen = vessel.WillChutesDeploy(altAGL, altASL, probableLandingSiteASL, pressure, t, parachuteSemiDeployMultiplier);
                    maxDragGees = Math.Max(maxDragGees, lastRecordedDrag.magnitude / 9.81f);

                    // If parachutes are about to open and we are running with a dt larger than the physics frame then drop dt to the physics frame rate and start again
                    if (willChutesOpen && dt > min_dt) // TODO test to see if we are better off just setting a minimum dt of the physics frame rate.
                    {
                        dt = Math.Max(dt / 2, min_dt);
                        repeatWithSmallerStep = true;
                    }
                    else
                    {
                        parachutesDeploying = vessel.Simulate(altAGL, altASL, probableLandingSiteASL, pressure, t, parachuteSemiDeployMultiplier);
                    }
                }
            }
            while (repeatWithSmallerStep);