defmodule Astro.Lunar do
@moduledoc """
Calulates lunar phases.
Each of the phases of the Moon is defined by the
angle between the Moon and Sun in the sky. When the Moon
is in between the Earth and the Sun, so that there is nearly a
zero degree separation, we see a New Moon.
Because the orbit of the Moon is tilted in relation to the
Earth’s orbit around the Sun, a New Moon can still be as much
as 5.2 degrees away from the Sun, thus why there isn't a
solar eclipse every month.
A crescent moon is 45 degrees from the Sun, a quarter moon
is 90 degrees from the Sun, a gibbous moon is 135 degrees
from the Sun, and the Full Moon is 180 degrees away from
the Sun.
"""
alias Astro.{Math, Time, Solar}
import Astro.Math, only: [
deg: 1,
sin: 1,
cos: 1,
mt: 1,
asin: 1,
sigma: 2,
mod: 2,
degrees: 1,
poly: 2,
invert_angular: 4
]
import Astro.Time, only: [
j2000: 0,
julian_centuries_from_moment: 1,
mean_synodic_month: 0
]
import Astro.Earth, only: [
nutation: 1
]
@months_epoch_to_j2000 24_724
@average_distance_earth_to_moon 385_000_560.0
@doc """
Returns the date time of the new
moon before a given moment.
## Arguments
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Returns
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Example
iex> Astro.Lunar.date_time_new_moon_before 738390
738375.5757777032
"""
@doc since: "0.5.0"
@spec date_time_new_moon_before(Time.moment()) :: Time.moment()
def date_time_new_moon_before(t) when is_number(t) do
t0 = nth_new_moon(0)
phi = lunar_phase_at(t)
n = round((t - t0) / mean_synodic_month() - phi / deg(360)) |> trunc()
nth_new_moon(Math.final(n - 1, &(nth_new_moon(&1) < t)))
end
@doc """
Returns the date time of the new
moon at or after a given date or
date time.
## Arguments
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Returns
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Example
iex> Astro.Lunar.date_time_new_moon_at_or_after 738390
738405.0352292997
"""
@doc since: "0.5.0"
@spec date_time_new_moon_at_or_after(Time.moment()) :: Time.moment()
def date_time_new_moon_at_or_after(t) when is_number(t) do
t0 = nth_new_moon(0)
phi = lunar_phase_at(t)
n = round((t - t0) / mean_synodic_month() - phi / deg(360.0))
nth_new_moon(Math.next(n, &(nth_new_moon(&1) >= t)))
end
@doc """
Returns the lunar phase as a
float number of degrees at a given
moment.
## Arguments
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Returns
* the lunar phase as a float number of
degrees.
## Example
iex> Astro.Lunar.lunar_phase_at 738389.5007195644
180.00001498208536
iex> Astro.Lunar.lunar_phase_at 738346.0544609067
0.021567106773019873
"""
@doc since: "0.5.0"
@spec lunar_phase_at(Time.moment()) :: Time.moment()
def lunar_phase_at(t) when is_number(t) do
phi = mod(lunar_longitude(t) - solar_longitude(t), 360)
t0 = nth_new_moon(0)
n = round((t - t0) / mean_synodic_month())
phi_prime = deg(360) * mod((t - nth_new_moon(n)) / mean_synodic_month(), 1)
if abs(phi - phi_prime) > deg(180.0) do
phi_prime
else
phi
end
end
@doc """
Returns the date time of a given
lunar phase at or before a given
moment.
## Arguments
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
* `phase` is the required lunar phase expressed
as a float number of degrees between `0.0` and
`360.0`
## Returns
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Example
iex> Astro.Lunar.date_time_lunar_phase_at_or_before(738368, Astro.Lunar.new_moon())
738346.0524695957
"""
@doc since: "0.5.0"
@spec date_time_lunar_phase_at_or_before(Time.moment(), Astro.phase()) :: Time.moment()
def date_time_lunar_phase_at_or_before(t, phase) do
tau = t - mean_synodic_month() * (1.0 / deg(360.0)) * mod(lunar_phase_at(t) - phase, 360.0)
a = tau - 2
b = min(t, tau + 2)
invert_angular(&lunar_phase_at/1, phase, a, b)
end
@doc """
Returns the date time of a given
lunar phase at or after a given
date time or date.
## Arguments
* a `moment` which is a float number of days
since `0000-01-01`
* `phase` is the required lunar phase expressed
as a float number of degrees between `0` and
`3660`
## Returns
* a `t:Time.moment()` which is a float number of days
since `0000-01-01`
## Example
iex> Astro.Lunar.date_time_lunar_phase_at_or_after(738368, Astro.Lunar.full_moon())
738389.5007195254
"""
@doc since: "0.5.0"
@spec date_time_lunar_phase_at_or_after(Time.moment(), Astro.phase()) :: Time.moment()
def date_time_lunar_phase_at_or_after(t, phase) do
tau = t + mean_synodic_month() * (1 / deg(360.0)) * mod(phase - lunar_phase_at(t), 360.0)
a = max(t, tau - 2)
b = tau + 2
invert_angular(&lunar_phase_at/1, phase, a, b)
end
@doc since: "0.6.0"
@spec lunar_position(Time.moment()) :: {Astro.angle(), Astro.angle(), Astro.meters()}
def lunar_position(t) do
lambda = lunar_longitude(t)
beta = lunar_latitude(t)
distance = lunar_distance(t)
{Astro.right_ascension(t, beta, lambda), Astro.declination(t, beta, lambda), distance}
end
@doc """
Returns the fractional illumination of the moon
at a given time as a fraction between 0.0 and 1.0.
"""
@doc since: "0.6.0"
@spec illuminated_fraction_of_moon(Time.time()) :: float()
def illuminated_fraction_of_moon(t) do
{a0, d0, r0} = lunar_position(t)
{a, d, r} = Solar.solar_position(t)
r = Math.au_to_m(r)
phi = :math.acos(sin(d0) * sin(d) + cos(d0) * cos(d) * cos(a0 - a))
i = Math.atan_r(r * :math.sin(phi), r0 - r * :math.cos(phi))
0.5 * (1 + :math.cos(i))
end
@doc """
Returns the new moon lunar
phase expressed as a float number
of degrees.
"""
@doc since: "0.5.0"
@spec new_moon() :: Astro.phase()
def new_moon() do
deg(0.0)
end
@doc """
Returns the full moon lunar
phase expressed as a float number
of degrees.
"""
@doc since: "0.5.0"
@spec full_moon() :: Astro.phase()
def full_moon() do
deg(180.0)
end
@doc """
Returns the first quarter lunar
phase expressed as a float number
of degrees.
"""
@doc since: "0.5.0"
@spec first_quarter() :: Astro.phase()
def first_quarter() do
deg(90.0)
end
@doc """
Returns the last quarter lunar
phase expressed as a float number
of degrees.
"""
@doc since: "0.5.0"
@spec last_quarter() :: Astro.phase()
def last_quarter() do
deg(270.0)
end
@doc false
@doc since: "0.5.0"
@spec lunar_longitude(Time.moment()) :: Astro.phase()
def lunar_longitude(t) do
c = julian_centuries_from_moment(t)
l = mean_lunar_longitude(c)
d = lunar_elongation(c)
m = solar_anomaly(c)
m_prime = lunar_anomaly(c)
f = moon_node(c)
e = poly(c, [1.0, -0.002516, -0.0000074])
args_lunar_elong = [
0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0, 4, 0, 4, 2, 2, 1,
1, 2, 2, 4, 2, 0, 2, 2, 1, 2, 0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2,
2, 2, 4, 0, 4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2
]
args_solar_anom = [
0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1,
0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2, 1, 2, -2, 0, 0, -1, 0, 0, 1,
-1, 2, 2, 1, -1, 0, 0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0
]
args_lunar_anom = [
1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 3, -2,
-1, 0, -1, 0, 1, 2, 0, -3, -2, -1, -2, 1, 0, 2, 0, -1, 1, 0,
-1, 2, -1, 1, -2, -1, -1, -2, 0, 1, 4, 0, -2, 0, 2, 1, -2, -3,
2, 1, -1, 3
]
args_moon_node = [
0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0,
0, 0, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0
]
sine_coeff = [
6288774, 1274027, 658314, 213618, -185116, -114332,
58793, 57066, 53322, 45758, -40923, -34720, -30383,
15327, -12528, 10980, 10675, 10034, 8548, -7888,
-6766, -5163, 4987, 4036, 3994, 3861, 3665, -2689,
-2602, 2390, -2348, 2236, -2120, -2069, 2048, -1773,
-1595, 1215, -1110, -892, -810, 759, -713, -700, 691,
596, 549, 537, 520, -487, -399, -381, 351, -340, 330,
327, -323, 299, 294
]
correction = deg(1 / 1000000) * sigma(
[sine_coeff, args_lunar_elong, args_solar_anom, args_lunar_anom, args_moon_node],
fn [v,w,x,y,z] ->
v * :math.pow(e, abs(x)) * sin(w * d + x * m + y * m_prime + z * f)
end
)
venus =
deg(3958 / 1000000) *
sin(deg(119.75) + c * deg(131.849))
jupiter =
deg(318 / 1000000) *
sin(deg(53.09) + c * deg(479_264.29))
flat_earth =
deg(1962 / 1000000) *
sin(l - f)
mod(l + correction + venus + jupiter + flat_earth + nutation(c), 360)
end
@doc since: "0.6.0"
@spec lunar_latitude(Time.moment()) :: Astro.angle()
def lunar_latitude(t) do
c = julian_centuries_from_moment(t)
l = mean_lunar_longitude(c)
d = lunar_elongation(c)
m = solar_anomaly(c)
m_prime = lunar_anomaly(c)
f = moon_node(c)
e = poly(c, [1.0, -0.002516, -0.0000074])
lunar_elongation = [
0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,
1,0,0,0,1,0,4,4,0,4,2,2,2,2,0,2,2,2,2,4,2,2,
0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2
]
solar_anomaly = [
0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,
0,1,0,1,1,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,1,1,
0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2]
lunar_anomaly = [
0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,
0,0,-1,0,1,1,0,0,3,0,-1,1,-2,0,2,1,-2,3,2,-3,
-1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-2,
0,0
]
moon_node = [
1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,
-1,1,3,1,1,1,-1,-1,-1,1,-1,1,-3,1,-3,-1,-1,1,
-1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,
-1,-1,-1,-1,-1,1
]
sine_coeff = [
5128122, 280602, 277693, 173237, 55413, 46271, 32573,
17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211,
2065, -1870, 1828, -1794, -1749, -1565, -1491, -1475,
-1410, -1344, -1335, 1107, 1021, 833, 777, 671, 607,
596, 491, -451, 439, 422, 421, -366, -351, 331, 315,
302, -283, -229, 223, 223, -220, -220, -185, 181,
-177, 176, 166, -164, 132, -119, 115, 107
]
beta = deg(1.0 / 1000000.0) * sigma(
[sine_coeff, lunar_elongation, solar_anomaly, lunar_anomaly, moon_node],
fn [v, w, x, y, z] ->
v * :math.pow(e, abs(x)) * sin(w*d + x*m + y*m_prime + z*f)
end
)
venus =
deg(175.0 / 1000000.0) *
sin(deg(119.75) + c * deg(131.849) + f) *
sin(deg(119.75) + c * deg(131.849) - f)
flat_earth =
deg(-2235.0 / 1000000.0) * sin(l) +
deg(127.0 / 1000000.0) * sin(l - m_prime) +
deg(-115.0 / 1000000.0) * sin(l + m_prime)
extra = deg(382.0 / 1000000.0) * sin(deg(313.45) + (c * deg(481266.484)))
beta + venus + flat_earth + extra
end
@doc since: "0.4.0"
@spec lunar_altitude(Time.moment(), Geo.PointZ.t()) :: Astro.angle()
def lunar_altitude(t, %Geo.PointZ{coordinates: {psi, phi, _alt}}) do
lambda = lunar_longitude(t)
beta = lunar_latitude(t)
alpha = Astro.right_ascension(t, beta, lambda)
delta = Astro.declination(t, beta, lambda)
theta = Time.sidereal_from_moment(t)
h = mod(theta + psi - alpha, 360.0)
altitude = asin(sin(phi) * sin(delta) + cos(phi) * cos(delta) * cos(h))
mod(altitude + deg(180.0), 360.0) - deg(180.0)
end
@doc since: "0.6.0"
@spec lunar_distance(Time.moment()) :: Astro.meters()
def lunar_distance(t) do
c = Time.julian_centuries_from_moment(t)
d = lunar_elongation(c)
m = solar_anomaly(c)
m_prime = lunar_anomaly(c)
f = moon_node(c)
e = poly(c, [1.0, -0.002516, -0.0000074])
lunar_elongation = [
0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,
1,2,2,4,2,0,2,2,1,2,0,0,2,2,2,4,0,3,2,4,0,2,
2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2
]
solar_anomaly = [
0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,
0,1,-1,0,0,0,1,0,-1,0,-2,1,2,-2,0,0,-1,0,0,1,
-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0
]
lunar_anomaly = [
1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,
-1,0,-1,0,1,2,0,-3,-2,-1,-2,1,0,2,0,-1,1,0,
-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,
2,1,-1,3,-1
]
moon_node = [
0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,
0,0,0,0,0,0,0,2,0,0,0,0,0,0,-2,2,0,2,0,0,0,0,
0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2
]
cos_coeff = [
-20905355,-3699111,-2955968,-569925,48888,-3149,
246158,-152138,-170733,-204586,-129620,108743,
104755,10321,0,79661,-34782,-23210,-21636,24208,
30824,-8379,-16675,-12831,-10445,-11650,14403,
-7003,0,10056,6322,-9884,5751,0,-4950,4130,0,
-3958,0,3258,2616,-1897,-2117,2354,0,0,-1423,
-1117,-1571,-1739,0,-4421,0,0,0,0,1165,0,0,
8752
]
correction = sigma(
[cos_coeff, lunar_elongation, solar_anomaly, lunar_anomaly, moon_node],
fn [v, w, x, y, z] ->
v * :math.pow(e, abs(x)) * cos((w * d) + (x * m) + (y * m_prime) + (z * f))
end
)
mt(@average_distance_earth_to_moon) + correction
end
@doc false
@doc since: "0.4.0"
@spec nth_new_moon(number()) :: Time.moment()
def nth_new_moon(n) do
k = n - @months_epoch_to_j2000
c = k / 1_236.85
approx =
j2000() + poly(c, [
5.09766, mean_synodic_month() * 1236.85, 0.0001437, -0.000000150, 0.00000000073
])
e = poly(c, [
1, -0.002516, -0.0000074
])
solar_anomaly = poly(c, [
2.5534, 1236.85 * 29.10535669, -0.0000014, -0.00000011
])
lunar_anomaly = poly(c, [
201.5643, 385.81693528 * 1236.85,
0.0107582, 0.00001238, -0.000000058
])
moon_argument = poly(c, [
160.7108, 390.67050284 * 1236.85,
-0.0016118, -0.00000227, 0.000000011
])
omega = poly(c, [
124.7746, -1.56375588 * 1236.85,
0.0020672, 0.00000215
])
e_factor = [
0, 1, 0, 0, 1, 1, 2, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0
]
solar_coeff = [
0, 1, 0, 0, -1, 1, 2, 0, 0, 1, 0, 1, 1, -1, 2,
0, 3, 1, 0, 1, -1, -1, 1, 0
]
lunar_coeff = [
1, 0, 2, 0, 1, 1, 0, 1, 1, 2, 3, 0, 0, 2, 1, 2,
0, 1, 2, 1, 1, 1, 3, 4
]
moon_coeff = [
0, 0, 0, 2, 0, 0, 0, -2, 2, 0, 0, 2, -2, 0, 0,
-2, 0, -2, 2, 2, 2, -2, 0, 0
]
sine_coeff = [
-0.40720, 0.17241, 0.01608,
0.01039, 0.00739, -0.00514,
0.00208, -0.00111, -0.00057,
0.00056, -0.00042, 0.00042,
0.00038, -0.00024, -0.00007,
0.00004, 0.00004, 0.00003,
0.00003, -0.00003, 0.00003,
-0.00002, -0.00002, 0.00002
]
correction =
deg(-0.00017) * sin(omega) +
sigma(
[sine_coeff, e_factor, solar_coeff, lunar_coeff, moon_coeff],
fn [v,w,x,y,z] ->
v * :math.pow(e, w) *
sin((x * solar_anomaly) + (y * lunar_anomaly) + (z * moon_argument))
end
)
extra =
deg(0.000325) *
sin(poly(c, [299.77, 132.8475848, -0.009173]))
add_const = [
251.88, 251.83, 349.42, 84.66, 141.74, 207.14, 154.84,
34.52, 207.19, 291.34, 161.72, 239.56, 331.55
]
add_coeff = [
0.016321, 26.651886, 36.412478, 18.206239, 53.303771,
2.453732, 7.306860, 27.261239, 0.121824, 1.844379,
24.198154, 25.513099, 3.592518
]
add_factor = [
0.000165, 0.000164, 0.000126, 0.000110, 0.000062, 0.000060,
0.000056, 0.000047, 0.000042, 0.000040, 0.000037, 0.000035,
0.000023
]
additional = sigma([add_const, add_coeff, add_factor], fn [i,j,l] -> l * sin(i + j * k) end)
Time.universal_from_dynamical(approx + correction + extra + additional)
end
def lunar_parallax(t, location) do
geo = lunar_altitude(t, location)
delta = lunar_distance(t)
alt = mt(6_378_140) / delta
arg = alt * cos(geo)
asin(arg)
end
def topocentric_lunar_altitude(t, location) do
lunar_altitude(t, location) - lunar_parallax(t, location)
end
@doc false
def mean_lunar_longitude(c) do
c
|> poly([218.3164477, 481267.88123421, -0.0015786, 1 / 538841.0, -1 /65194000.0])
|> degrees()
end
@doc false
def lunar_elongation(c) do
c
|> poly([297.8501921, 445267.1114034, -0.0018819, 1/545868, -1 / 113065000.0])
|> degrees()
end
@doc false
def solar_anomaly(c) do
c
|> poly([357.5291092, 35999.0502909, -0.0001536, 1 / 24490000.0])
|> degrees()
end
@doc false
def lunar_anomaly(c) do
c
|> poly([134.9633964, 477198.8675055, 0.0087414, 1 / 69699.0, -1 / 14712000.0])
|> degrees()
end
defp solar_longitude(t) do
c = julian_centuries_from_moment(t)
Astro.Solar.sun_apparent_longitude_alt(c)
end
@doc false
def lunar_node(t) do
c = julian_centuries_from_moment(t)
moon_node(c + deg(90.0))
|> mod(180.0)
|> Kernel.-(90.0)
end
@doc false
def moon_node(c) do
c
|> poly([93.2720950, 483202.0175233, -0.0036539, -1 / 3526000.0, 1 / 863310000.0])
|> degrees()
end
end
