private double evolveDistance(final KinematicModel km, final double dist, final double deltaH, final double fAir) {
        // elevation force
        final double fh = deltaH * km.totalweight * 9.81 / dist;

        final double effectiveSpeedLimit = km.getEffectiveSpeedLimit();
        final double emax = 0.5 * km.totalweight * effectiveSpeedLimit * effectiveSpeedLimit;
        if (emax <= 0.) {
            return -1.;
        }
        final double vb = km.getBreakingSpeed(effectiveSpeedLimit);
        final double elow = 0.5 * km.totalweight * vb * vb;

        double elapsedTime = 0.;
        double dissipatedEnergy = 0.;

        double v = Math.sqrt(2. * ekin / km.totalweight);
        double d = dist;
        while (d > 0.) {
            final boolean slow = ekin < elow;
            final boolean fast = ekin >= emax;
            final double etarget = slow ? elow : emax;
            double f = km.fRoll + fAir * v * v + fh;
            final double fRecup = Math.max(0., fast ? -f : (slow ? km.fRecup : 0) - fh); // additional recup for slow part
            f += fRecup;

            double deltaEkin;
            final double timeStep;
            double x;
            if (fast) {
                x = d;
                deltaEkin = x * f;
                timeStep = x / v;
                ekin = etarget;
            } else {
                deltaEkin = etarget - ekin;
                final double b = 2. * fAir / km.totalweight;
                final double x0 = deltaEkin / f;
                final double x0b = x0 * b;
                x = x0 * (1. - x0b * (0.5 + x0b * (0.333333333 - x0b * 0.25))); // = ln( deltaEkin*b/f + 1.) / b;
                final double maxstep = Math.min(50., d);
                if (x >= maxstep) {
                    x = maxstep;
                    final double xb = x * b;
                    deltaEkin = x * f * (1. + xb * (0.5 + xb * (0.166666667 + xb * 0.0416666667))); // = f/b* exp(xb-1)
                    ekin += deltaEkin;
                } else {
                    ekin = etarget;
                }
                final double v2 = Math.sqrt(2. * ekin / km.totalweight);
                final double a = f / km.totalweight; // TODO: average force?
                timeStep = (v2 - v) / a;
                v = v2;
            }
            d -= x;
            elapsedTime += timeStep;

            // dissipated energy does not contain elevation and efficient recup
            dissipatedEnergy += deltaEkin - x * (fh + fRecup * km.recupEfficiency);

            // correction: inefficient recup going into heating is half efficient
            final double ieRecup = x * fRecup * (1. - km.recupEfficiency);
            final double eaux = timeStep * km.pStandby;
            dissipatedEnergy -= Math.max(ieRecup, eaux) * 0.5;
        }

        dissipatedEnergy += elapsedTime * km.pStandby;

        totalTime += elapsedTime;
        totalEnergy += dissipatedEnergy + dist * fh;

        return (km.pw * elapsedTime + dissipatedEnergy) / km.cost0; // =cost
    }
The path-instance of the kinematic model

@author ab
main/src/main/java/cgeo/geocaching/brouter/core/KinematicPath.java