001/*
002 * Shredzone Commons - suncalc
003 *
004 * Copyright (C) 2018 Richard "Shred" Körber
005 *   http://commons.shredzone.org
006 *
007 * Licensed under the Apache License, Version 2.0 (the "License");
008 * you may not use this file except in compliance with the License.
009 *
010 * This program is distributed in the hope that it will be useful,
011 * but WITHOUT ANY WARRANTY; without even the implied warranty of
012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
013 */
014package org.shredzone.commons.suncalc.util;
015
016import static java.lang.Math.abs;
017
018import java.util.function.Function;
019
020/**
021 * Finds the root of a function by using the Pegasus method.
022 *
023 * @see <a href="https://en.wikipedia.org/wiki/False_position_method">regula falsi</a>
024 */
025public class Pegasus {
026
027    private static final int MAX_ITERATIONS = 30;
028
029    /**
030     * Find the root of the given function within the boundaries.
031     *
032     * @param lower
033     *            Lower boundary
034     * @param upper
035     *            Upper boundary
036     * @param accuracy
037     *            Desired accuracy
038     * @param f
039     *            Function to be used for calculation
040     * @return root that was found
041     * @throws ArithmeticException
042     *             if the root could not be found in the given accuracy within
043     *             {@value #MAX_ITERATIONS} iterations.
044     */
045    public static Double calculate(double lower, double upper, double accuracy, Function<Double, Double> f) {
046        double x1 = lower;
047        double x2 = upper;
048
049        double f1 = f.apply(x1);
050        double f2 = f.apply(x2);
051
052        if (f1 * f2 >= 0.0) {
053            throw new ArithmeticException("No root within the given boundaries");
054        }
055
056        int i = MAX_ITERATIONS;
057
058        while (i-- > 0) {
059            double x3 = x2 - f2 / ((f2 - f1) / (x2 - x1));
060            double f3 = f.apply(x3);
061
062            if (f3 * f2 <= 0.0) {
063                x1 = x2;
064                f1 = f2;
065                x2 = x3;
066                f2 = f3;
067            } else {
068                f1 = f1 * f2 / (f2 + f3);
069                x2 = x3;
070                f2 = f3;
071            }
072
073            if (abs(x2 - x1) <= accuracy) {
074                return abs(f1) < abs(f2) ? x1 : x2;
075            }
076        }
077
078        throw new ArithmeticException("Maximum number of iterations exceeded");
079    }
080
081}