Solve W=C tan W using Newton-Raphson - for https://stackoverflow.com/q/46289473

SO46289473.java

/**
 * A Newton-Raphson implementation of https://stackoverflow.com/q/46289473
 */
public class SO46289473 
{
    public static double nrSolveW(double c, double epsilon) {
        double guessValue;
        if (c < 0 && c != Double.NEGATIVE_INFINITY) {
            // Then 1st intersection is between π/2 and π (towards -infinite C)
            guessValue = Math.PI * 0.75;
        } else if (c > 0 && c < 1) {
            // Than 1st intersection is between  π/2 and 0
            guessValue = Math.PI * 0.25;
        } else if (c > 1 && c != Double.POSITIVE_INFINITY) {
            // Than 1st intersection is between NOT 3π/2, but 4.493 and π/2 (towards infinite C)
            guessValue = 3.032;
        } else {
            // C = Infinity, so it's not possible to fit a sinusoidal equation between these two points :/
            return Double.NaN;
        }

        // Newton-Raphson method to solve f(w) = w - c tan w = 0
        // where f'(w) = df/dw = -c * (cos w)^-2
        // and next w = w - (f(w)/f'(w))

        double w = guessValue;
        double lastAbsF = 1000000.0;
        for (int iteration = 0; iteration < 1000; ++iteration) {
            double f = w - c * Math.tan(w);

            double absF = Math.abs(f);
            if (absF < epsilon) {
                // Within margin of accuracy required
                return w;
            } else if (absF > lastAbsF) {
                // We're diverging away from the solution
                throw new RuntimeException("Newton-Raphson diverging: c = " + c);
            }
            lastAbsF = absF;

            double cosW = Math.cos(w);
            double correction = f * cosW * cosW / c;
            w += correction;
        }

        throw new RuntimeException("Newton-Raphson did not converge: c = " + c);
    }

    public static void main( String[] args )
    {
        double c = 0.4;
        double epsilon = Math.pow(2.0, -40);
        System.out.println(nrSolveW(c, epsilon));
    }
}