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 Sun.
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 Sun {
027
028    private static final double SUN_DISTANCE = 149598000.0;
029    private static final double SUN_MEAN_RADIUS = 695700.0;
030
031    private Sun() {
032        // Utility class without constructor
033    }
034
035    /**
036     * Calculates the equatorial position of the sun.
037     *
038     * @param date
039     *            {@link JulianDate} to be used
040     * @return {@link Vector} containing the sun position
041     */
042    public static Vector positionEquatorial(JulianDate date) {
043        double T = date.getJulianCentury();
044        double M = PI2 * frac(0.993133 + 99.997361 * T);
045        double L = PI2 * frac(0.7859453 + M / PI2
046            + (6893.0 * sin(M) + 72.0 * sin(2.0 * M) + 6191.2 * T) / 1296.0e3);
047
048        double d = SUN_DISTANCE
049            * (1 - 0.016718 * cos(date.getTrueAnomaly()));
050
051        return Vector.ofPolar(L, 0.0, d);
052    }
053
054    /**
055     * Calculates the geocentric position of the sun.
056     *
057     * @param date
058     *            {@link JulianDate} to be used
059     * @return {@link Vector} containing the sun position
060     */
061    public static Vector position(JulianDate date) {
062        Matrix rotateMatrix = equatorialToEcliptical(date).transpose();
063        return rotateMatrix.multiply(positionEquatorial(date));
064    }
065
066    /**
067     * Calculates the horizontal position of the sun.
068     *
069     * @param date
070     *            {@link JulianDate} to be used
071     * @param lat
072     *            Latitude, in radians
073     * @param lng
074     *            Longitute, in radians
075     * @return {@link Vector} of horizontal sun position
076     */
077    public static Vector positionHorizontal(JulianDate date, double lat, double lng) {
078        Vector mc = position(date);
079        double h = date.getGreenwichMeanSiderealTime() + lng - mc.getPhi();
080        return equatorialToHorizontal(h, mc.getTheta(), mc.getR(), lat);
081    }
082
083    /**
084     * Calculates the topocentric position of the sun.
085     * <p>
086     * Atmospheric refraction is <em>not</em> taken into account.
087     *
088     * @param date
089     *            {@link JulianDate} to be used
090     * @param lat
091     *            Latitude, in radians
092     * @param lng
093     *            Longitute, in radians
094     * @param elev
095     *            Elevation, in meters
096     * @return {@link Vector} of topocentric sun position
097     * @since 3.9
098     */
099    public static Vector positionTopocentric(JulianDate date, double lat, double lng, double elev) {
100        Vector pos = positionHorizontal(date, lat, lng);
101        return Vector.ofPolar(
102                pos.getPhi(),
103                pos.getTheta() - parallax(elev, pos.getR()),
104                pos.getR()
105        );
106    }
107
108    /**
109     * Returns the angular radius of the sun.
110     *
111     * @param distance
112     *            Distance of the sun, in kilometers.
113     * @return Angular radius of the sun, in radians.
114     * @see <a href="https://en.wikipedia.org/wiki/Angular_diameter">Wikipedia: Angular
115     *      Diameter</a>
116     */
117    public static double angularRadius(double distance) {
118        return asin(SUN_MEAN_RADIUS / distance);
119    }
120
121}