001/*
002 * Shredzone Commons - suncalc
003 *
004 * Copyright (C) 2017 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.*;
017
018import java.util.Comparator;
019import java.util.function.Function;
020
021/**
022 * Contains constants and mathematical operations that are not available in {@link Math}.
023 */
024public final class ExtendedMath {
025
026    /**
027     * PI * 2
028     */
029    public static final double PI2 = PI * 2.0;
030
031    /**
032     * Arc-Seconds per Radian.
033     */
034    public static final double ARCS = toDegrees(3600.0);
035
036    /**
037     * Mean radius of the earth, in kilometers.
038     */
039    public static final double EARTH_MEAN_RADIUS = 6371.0;
040
041    /**
042     * Refraction at the horizon, in radians.
043     */
044    public static final double REFRACTION_AT_HORIZON = PI / (tan(toRadians(7.31 / 4.4)) * 10800.0);
045
046    private ExtendedMath() {
047        // utility class without constructor
048    }
049
050    /**
051     * Returns the decimal part of a value.
052     *
053     * @param a
054     *            Value
055     * @return Fraction of that value. It has the same sign as the input value.
056     */
057    public static double frac(double a) {
058        return a % 1.0;
059    }
060
061    /**
062     * Performs a safe check if the given double is actually zero (0.0).
063     * <p>
064     * Note that "almost zero" returns {@code false}, so this method should not be used
065     * for comparing calculation results to zero.
066     *
067     * @param d
068     *            double to check for zero.
069     * @return {@code true} if the value was zero, or negative zero.
070     */
071    public static boolean isZero(double d) {
072        // This should keep squid:S1244 silent...
073        return !Double.isNaN(d) && round(signum(d)) == 0L;
074    }
075
076    /**
077     * Converts equatorial coordinates to horizontal coordinates.
078     *
079     * @param tau
080     *            Hour angle (radians)
081     * @param dec
082     *            Declination (radians)
083     * @param dist
084     *            Distance of the object
085     * @param lat
086     *            Latitude of the observer (radians)
087     * @return {@link Vector} containing the horizontal coordinates
088     */
089    public static Vector equatorialToHorizontal(double tau, double dec, double dist, double lat) {
090        return Matrix.rotateY(PI / 2.0 - lat).multiply(Vector.ofPolar(tau, dec, dist));
091    }
092
093    /**
094     * Creates a rotational {@link Matrix} for converting equatorial to ecliptical
095     * coordinates.
096     *
097     * @param t
098     *            {@link JulianDate} to use
099     * @return {@link Matrix} for converting equatorial to ecliptical coordinates
100     */
101    public static Matrix equatorialToEcliptical(JulianDate t) {
102        double jc = t.getJulianCentury();
103        double eps = toRadians(23.43929111 - (46.8150 + (0.00059 - 0.001813 * jc) * jc) * jc / 3600.0);
104        return Matrix.rotateX(eps);
105    }
106
107    /**
108     * Returns the parallax for objects at the horizon.
109     *
110     * @param elevation
111     *            Observer's elevation, in meters above sea level. Must not be negative.
112     * @param distance
113     *            Distance of the object, in kilometers.
114     * @return parallax, in radians
115     */
116    public static double parallax(double elevation, double distance) {
117        return asin(EARTH_MEAN_RADIUS / distance)
118             - acos(EARTH_MEAN_RADIUS / (EARTH_MEAN_RADIUS + (elevation / 1000.0)));
119    }
120
121    /**
122     * Calculates the atmospheric refraction of an object at the given apparent altitude.
123     * <p>
124     * The result is only valid for positive altitude angles. If negative, 0.0 is
125     * returned.
126     * <p>
127     * Assumes an atmospheric pressure of 1010 hPa and a temperature of 10 °C.
128     *
129     * @param ha
130     *            Apparent altitude, in radians.
131     * @return Refraction at this altitude
132     * @see <a href="https://en.wikipedia.org/wiki/Atmospheric_refraction">Wikipedia:
133     *      Atmospheric Refraction</a>
134     */
135    public static double apparentRefraction(double ha) {
136        if (ha < 0.0) {
137            return 0.0;
138        }
139
140        if (isZero(ha)) {
141            return REFRACTION_AT_HORIZON;
142        }
143
144        return PI / (tan(toRadians(ha + (7.31 / (ha + 4.4)))) * 10800.0);
145    }
146
147    /**
148     * Calculates the atmospheric refraction of an object at the given altitude.
149     * <p>
150     * The result is only valid for positive altitude angles. If negative, 0.0 is
151     * returned.
152     * <p>
153     * Assumes an atmospheric pressure of 1010 hPa and a temperature of 10 °C.
154     *
155     * @param h
156     *            True altitude, in radians.
157     * @return Refraction at this altitude
158     * @see <a href="https://en.wikipedia.org/wiki/Atmospheric_refraction">Wikipedia:
159     *      Atmospheric Refraction</a>
160     */
161    public static double refraction(double h) {
162        if (h < 0.0) {
163            return 0.0;
164        }
165
166        // refraction formula, converted to radians
167        return 0.000296706 / tan(h + 0.00312537 / (h + 0.0890118));
168    }
169
170    /**
171     * Converts dms to double.
172     *
173     * @param d
174     *            Degrees. Sign is used for result.
175     * @param m
176     *            Minutes. Sign is ignored.
177     * @param s
178     *            Seconds and fractions. Sign is ignored.
179     * @return angle, in degrees
180     */
181    public static double dms(int d, int m, double s) {
182        double sig = d < 0 ? -1.0 : 1.0;
183        return sig * ((abs(s) / 60.0 + abs(m)) / 60.0 + abs(d));
184    }
185
186    /**
187     * Locates the true maximum within the given time frame.
188     *
189     * @param time
190     *         Base time
191     * @param frame
192     *         Time frame, which is added to and subtracted from the base time for the
193     *         interval
194     * @param depth
195     *         Maximum recursion depth. For each recursion, the function is invoked once.
196     * @param f
197     *         Function to be used for calculation
198     * @return time of the true maximum
199     */
200    public static double readjustMax(double time, double frame, int depth, Function<Double, Double> f) {
201        double left = time - frame;
202        double right = time + frame;
203        double leftY = f.apply(left);
204        double rightY = f.apply(right);
205
206        return readjustInterval(left, right, leftY, rightY, depth, f, Double::compare);
207    }
208
209    /**
210     * Locates the true minimum within the given time frame.
211     *
212     * @param time
213     *         Base time
214     * @param frame
215     *         Time frame, which is added to and subtracted from the base time for the
216     *         interval
217     * @param depth
218     *         Maximum recursion depth. For each recursion, the function is invoked once.
219     * @param f
220     *         Function to be used for calculation
221     * @return time of the true minimum
222     */
223    public static double readjustMin(double time, double frame, int depth, Function<Double, Double> f) {
224        double left = time - frame;
225        double right = time + frame;
226        double leftY = f.apply(left);
227        double rightY = f.apply(right);
228
229        return readjustInterval(left, right, leftY, rightY, depth, f, (yl, yr) -> Double.compare(yr, yl));
230    }
231
232    /**
233     * Recursively find the true maximum/minimum within the given time frame.
234     *
235     * @param left
236     *         Left interval border
237     * @param right
238     *         Right interval border
239     * @param yl
240     *         Function result at the left interval
241     * @param yr
242     *         Function result at the right interval
243     * @param depth
244     *         Maximum recursion depth. For each recursion, the function is invoked once.
245     * @param f
246     *         Function to invoke
247     * @param cmp
248     *         Comparator to decide whether the left or right side of the interval half is
249     *         to be used
250     * @return Position of the approximated minimum/maximum
251     */
252    private static double readjustInterval(double left, double right, double yl, double yr, int depth,
253                                           Function<Double, Double> f, Comparator<Double> cmp) {
254        if (depth <= 0) {
255            return (cmp.compare(yl, yr) < 0) ? right : left;
256        }
257
258        double middle = (left + right) / 2.0;
259        double ym = f.apply(middle);
260        if (cmp.compare(yl, yr) < 0) {
261            return readjustInterval(middle, right, ym, yr, depth - 1, f, cmp);
262        } else {
263            return readjustInterval(left, middle, yl, ym, depth - 1, f, cmp);
264        }
265    }
266
267}