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
018/**
019 * Finds the root of a function by using the Pegasus method.
020 *
021 * @see <a href="https://en.wikipedia.org/wiki/False_position_method">regula falsi</a>
022 * @since 2.3
023 */
024public class Pegasus {
025
026    private static final int MAX_ITERATIONS = 30;
027
028    /**
029     * Find the root of the given function within the boundaries.
030     *
031     * @param lower
032     *            Lower boundary
033     * @param upper
034     *            Upper boundary
035     * @param accuracy
036     *            Desired accuracy
037     * @param f
038     *            Function to be used for calculation
039     * @return root that was found
040     * @throws ArithmeticException
041     *             if the root could not be found in the given accuracy within
042     *             {@value #MAX_ITERATIONS} iterations.
043     */
044    public static Double calculate(double lower, double upper, double accuracy, Function f) {
045        double x1 = lower;
046        double x2 = upper;
047
048        double f1 = f.apply(x1);
049        double f2 = f.apply(x2);
050
051        if (f1 * f2 >= 0.0) {
052            throw new ArithmeticException("No root within the given boundaries");
053        }
054
055        int i = MAX_ITERATIONS;
056
057        while (i-- > 0) {
058            double x3 = x2 - f2 / ((f2 - f1) / (x2 - x1));
059            double f3 = f.apply(x3);
060
061            if (f3 * f2 <= 0.0) {
062                x1 = x2;
063                f1 = f2;
064                x2 = x3;
065                f2 = f3;
066            } else {
067                f1 = f1 * f2 / (f2 + f3);
068                x2 = x3;
069                f2 = f3;
070            }
071
072            if (abs(x2 - x1) <= accuracy) {
073                return abs(f1) < abs(f2) ? x1 : x2;
074            }
075        }
076
077        throw new ArithmeticException("Maximum number of iterations exceeded");
078    }
079
080    /**
081     * The function that is to be solved.
082     *
083     * @since 2.3
084     */
085    public interface Function {
086
087        /**
088         * Calculate the function result for x.
089         *
090         * @param x
091         *            x
092         * @return f(x)
093         */
094        double apply(double x);
095
096    }
097
098}