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.*;
017import static org.shredzone.commons.suncalc.util.ExtendedMath.*;
018
019/**
020 * Calculations and constants for the Moon.
021 *
022 * @see "Astronomy on the Personal Computer, 4th edition
023 *      (Oliver Montenbruck, Thomas Pfleger) -
024 *      ISBN 978-3-540-67221-0"
025 */
026public final class Moon {
027
028    private static final double MOON_MEAN_RADIUS = 1737.1;
029
030    private Moon() {
031        // Utility class without constructor
032    }
033
034    /**
035     * Calculates the equatorial position of the moon.
036     *
037     * @param date
038     *            {@link JulianDate} to be used
039     * @return {@link Vector} of equatorial moon position
040     */
041    public static Vector positionEquatorial(JulianDate date) {
042        double T  = date.getJulianCentury();
043        double L0 =       frac(0.606433 + 1336.855225 * T);
044        double l  = PI2 * frac(0.374897 + 1325.552410 * T);
045        double ls = PI2 * frac(0.993133 +   99.997361 * T);
046        double D  = PI2 * frac(0.827361 + 1236.853086 * T);
047        double F  = PI2 * frac(0.259086 + 1342.227825 * T);
048        double D2 = 2.0 * D;
049        double l2 = 2.0 * l;
050        double F2 = 2.0 * F;
051
052        double dL = 22640.0 * sin(l)
053                  -  4586.0 * sin(l - D2)
054                  +  2370.0 * sin(D2)
055                  +   769.0 * sin(l2)
056                  -   668.0 * sin(ls)
057                  -   412.0 * sin(F2)
058                  -   212.0 * sin(l2 - D2)
059                  -   206.0 * sin(l + ls - D2)
060                  +   192.0 * sin(l + D2)
061                  -   165.0 * sin(ls - D2)
062                  -   125.0 * sin(D)
063                  -   110.0 * sin(l + ls)
064                  +   148.0 * sin(l - ls)
065                  -    55.0 * sin(F2 - D2);
066
067        double S  = F + (dL + 412.0 * sin(F2) + 541.0 * sin(ls)) / ARCS;
068        double h  = F - D2;
069        double N  =  -526.0 * sin(h)
070                  +    44.0 * sin(l + h)
071                  -    31.0 * sin(-l + h)
072                  -    23.0 * sin(ls + h)
073                  +    11.0 * sin(-ls + h)
074                  -    25.0 * sin(-l2 + F)
075                  +    21.0 * sin(-l + F);
076
077        double l_Moon = PI2 * frac(L0 + dL / 1296.0e3);
078        double b_Moon = (18520.0 * sin(S) + N) / ARCS;
079
080        double dt = 385000.5584
081                  -  20905.3550 * cos(l)
082                  -   3699.1109 * cos(D2 - l)
083                  -   2955.9676 * cos(D2)
084                  -    569.9251 * cos(l2);
085
086        return Vector.ofPolar(l_Moon, b_Moon, dt);
087    }
088
089    /**
090     * Calculates the geocentric position of the moon.
091     *
092     * @param date
093     *            {@link JulianDate} to be used
094     * @return {@link Vector} of geocentric moon position
095     */
096    public static Vector position(JulianDate date) {
097        Matrix rotateMatrix = equatorialToEcliptical(date).transpose();
098        return rotateMatrix.multiply(positionEquatorial(date));
099    }
100
101    /**
102     * Calculates the horizontal position of the moon.
103     *
104     * @param date
105     *            {@link JulianDate} to be used
106     * @param lat
107     *            Latitude, in radians
108     * @param lng
109     *            Longitute, in radians
110     * @return {@link Vector} of horizontal moon position
111     */
112    public static Vector positionHorizontal(JulianDate date, double lat, double lng) {
113        Vector mc = position(date);
114        double h = date.getGreenwichMeanSiderealTime() + lng - mc.getPhi();
115        return equatorialToHorizontal(h, mc.getTheta(), mc.getR(), lat);
116    }
117
118    /**
119     * Calculates the topocentric position of the moon.
120     * <p>
121     * Atmospheric refraction is <em>not</em> taken into account.
122     *
123     * @param date
124     *            {@link JulianDate} to be used
125     * @param lat
126     *            Latitude, in radians
127     * @param lng
128     *            Longitute, in radians
129     * @param elev
130     *            Elevation, in meters
131     * @return {@link Vector} of topocentric moon position
132     * @since 3.9
133     */
134    public static Vector positionTopocentric(JulianDate date, double lat, double lng, double elev) {
135        Vector pos = positionHorizontal(date, lat, lng);
136        return Vector.ofPolar(
137                pos.getPhi(),
138                pos.getTheta() - parallax(elev, pos.getR()),
139                pos.getR()
140        );
141    }
142
143    /**
144     * Returns the angular radius of the moon.
145     *
146     * @param distance
147     *            Distance of the moon, in kilometers.
148     * @return Angular radius of the moon, in radians.
149     * @see <a href="https://en.wikipedia.org/wiki/Angular_diameter">Wikipedia: Angular
150     *      Diameter</a>
151     */
152    public static double angularRadius(double distance) {
153        return asin(MOON_MEAN_RADIUS / distance);
154    }
155
156}