/*
Copyright (C) 2024 Samuel Jimenez <therealsamueljimenez@gmail.com>
JavaScript version
Copyright 2011 Glad Deschrijver <glad.deschrijver@gmail.com>
with the same license as below.
*/
/*
This file is part of the kmoon application with explicit permission by the author
Copyright 1996 Christopher Osburn <chris@speakeasy.org>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
.pragma library
// var LUNATION_OFFSET = 953//lunation of LUNATION_EPOCH in Brown Lunation Number
const LUNATION_OFFSET = 0//lunation of LUNATION_EPOCH using Meeus's Lunation Number
// LUNATION_EPOCH
//date of first new moon of 2000 (2000-01-06 18:14:00)
const LUNATION_EPOCH_MS = 947182440000// number of milliseconds between 1970-01-01 00:00:00 and epoch
const LUNATION_EPOCH_JD = 2451550.259722222//epoch in JD
// milliseconds in one day
const MS_PER_DAY = 86400000
// average lunation duration
const LUNATION_CYCLE_DURATION = 29.530588853
const LUNATION_CYCLE_DURATION_MS = 2551442877
/* DatefromJD(int) => Date
*
* convert a Julian Date to a JavaScript Date object
*/
function DatefromJD(jd)
{
// convert to epoch time, then instantiate Date
return new Date((jd - 2440587.5) * MS_PER_DAY)
}
/* radian(x) : x in Reals -> [-2 PI, 2 PI]
*
* convert x to radians from degrees
*/
function radian(x)
{
return ((x % 360.0) * Math.PI/180)
}
/* getLunation(Date time) => int lunation
*
* Given time, returns number of lunation time falls in
*
*/
function getLunation(time)
{
// find lunation
var lunation = Math.floor((time.getTime() - LUNATION_EPOCH_MS) / LUNATION_CYCLE_DURATION_MS) + LUNATION_OFFSET
//verify result
var currNew = DatefromJD(moonphasebylunation(lunation, 0))
var nextNew
while (currNew > time) {
// math is off: log result
console.log("Lunation number", lunation, "too high for", time.toISOString())
lunation--
currNew = DatefromJD(moonphasebylunation(lunation, 0))
}
nextNew = DatefromJD(moonphasebylunation(lunation + 1, 0))
while (nextNew < time) {
// math is off: log result
console.log("Lunation number", lunation, "too low for", time.toISOString())
lunation++
nextNew = DatefromJD(moonphasebylunation(lunation + 1, 0))
}
return lunation;
}
/*
** moonphase.c
** 1996/02/11
**
** Copyright 1996, Christopher Osburn, Lunar Outreach Services,
** Non-commercial usage license granted to all.
**
** calculate phase of the moon per Meeus Ch. 47 (ISBN:9780943396613)
**
** Parameters:
** int lun: phase parameter. This is the number of lunations
** since the New Moon of 2000 January 6.
**
** int phi: another phase parameter, selecting the phase of the
** moon. 0 = New, 1 = First Qtr, 2 = Full, 3 = Last Qtr
**
** Return: Apparent JD of the needed phase
*/
function moonphasebylunation(lun, phi)
{
var k = lun - LUNATION_OFFSET + phi / 4.0
// Julian Ephemeris Day of phase event
var JDE
// time parameter, Julian Centuries since LUNATION_EPOCH (47.3)
var T = k / 1236.85
// Eccentricity anomaly (45.6)
var E = 1.0 + T * (-0.002516 + -0.0000074 * T)
// Sun's mean anomaly (47.4)
var M = radian(2.5534 + 29.10535669 * k
+ T * T * (-0.0000218 + -0.00000011 * T))
// Moon's mean anomaly (47.5)
var M1 = radian(201.5643 + 385.81693528 * k
+ T * T * (0.0107438 + T * (0.00001239 + -0.000000058 * T)))
// Moon's argument of latitude (47.6)
var F = radian(160.7108 + 390.67050274 * k
+ T * T * (-0.0016341 * T * (-0.00000227 + 0.000000011 * T)))
// Moon's longitude of ascending node (47.7)
var O = radian(124.7746 - 1.56375580 * k
+ T * T * (0.0020691 + 0.00000215 * T))
// planetary arguments
var planetaryArguments =
0.000325 * Math.sin(radian(299.77 + 0.107408 * k - 0.009173 * T * T))
+ 0.000165 * Math.sin(radian(251.88 + 0.016321 * k))
+ 0.000164 * Math.sin(radian(251.83 + 26.651886 * k))
+ 0.000126 * Math.sin(radian(349.42 + 36.412478 * k))
+ 0.000110 * Math.sin(radian(84.66 + 18.206239 * k))
+ 0.000062 * Math.sin(radian(141.74 + 53.303771 * k))
+ 0.000060 * Math.sin(radian(207.14 + 2.453732 * k))
+ 0.000056 * Math.sin(radian(154.84 + 7.306860 * k))
+ 0.000047 * Math.sin(radian(34.52 + 27.261239 * k))
+ 0.000042 * Math.sin(radian(207.19 + 0.121824 * k))
+ 0.000040 * Math.sin(radian(291.34 + 1.844379 * k))
+ 0.000037 * Math.sin(radian(161.72 + 24.198154 * k))
+ 0.000035 * Math.sin(radian(239.56 + 25.513099 * k))
+ 0.000023 * Math.sin(radian(331.55 + 3.592518 * k))
// added correction for quarter phases
var W
//group First, Last Quarter
var orbitalCorrectionIndex = 2-Math.abs(2-phi)
var orbitalCorrectionCoeffArray = [
// New Moon
[+0.17241, +0.00208, +0.00004, -0.40720, +0.00739, -0.00514, +0.00000, -0.00007, -0.00111, -0.00057, +0.01608, +0.01039, -0.00024, +0.00056, -0.00042, -0.00002, +0.00002, +0.00038, +0.00042, -0.00017, +0.00004, +0.00003, -0.00002, +0.00003, +0.00003, -0.00003,],
// First Quarter, Last Quarter,
[+0.17172, +0.00204, +0.00003, -0.62801, +0.00454, -0.01183, +0.00004, -0.00028* E * E, -0.00180, -0.00070, +0.00862, +0.00804, -0.00034, +0.00027, -0.00040, -0.00002, +0.00000, +0.00032, +0.00032, -0.00017, +0.00002, +0.00004, -0.00005, +0.00002, +0.00003, -0.00004,],
// New Moon
[+0.17302, +0.00209, +0.00004, -0.40614, +0.00734, -0.00515, +0.00000, -0.00007, -0.00111, -0.00057, +0.01614, +0.01043, -0.00024, +0.00056, -0.00042, -0.00002, +0.00002, +0.00038, +0.00042, -0.00017, +0.00004, +0.00003, -0.00002, +0.00003, +0.00003, -0.00003,]]
var orbitalCorrectionCoeffs =
orbitalCorrectionCoeffArray[orbitalCorrectionIndex]
var orbitalApproximation =
orbitalCorrectionCoeffs[0] * E * Math.sin(M)
+ orbitalCorrectionCoeffs[1] * E * E * Math.sin(2.0 * M)
+ orbitalCorrectionCoeffs[2] * Math.sin(3.0 * M)
+ orbitalCorrectionCoeffs[3] * Math.sin(M1)
+ orbitalCorrectionCoeffs[4] * E * Math.sin(M1 - M)
+ orbitalCorrectionCoeffs[5] * E * Math.sin(M1 + M)
+ orbitalCorrectionCoeffs[6] * Math.sin(M1 - 2.0 * M)
+ orbitalCorrectionCoeffs[7] * Math.sin(M1 + 2.0 * M)
+ orbitalCorrectionCoeffs[8] * Math.sin(M1 - 2.0 * F)
+ orbitalCorrectionCoeffs[9] * Math.sin(M1 + 2.0 * F)
+ orbitalCorrectionCoeffs[10] * Math.sin(2.0 * M1)
+ orbitalCorrectionCoeffs[11] * Math.sin(2.0 * F)
+ orbitalCorrectionCoeffs[12] * E * Math.sin(2.0 * M1 - M)
+ orbitalCorrectionCoeffs[13] * E * Math.sin(2.0 * M1 + M)
+ orbitalCorrectionCoeffs[14] * Math.sin(3.0 * M1)
+ orbitalCorrectionCoeffs[15] * Math.sin(3.0 * M1 + M)
+ orbitalCorrectionCoeffs[16] * Math.sin(4.0 * M1)
+ orbitalCorrectionCoeffs[17] * E * Math.sin(M - 2.0 * F)
+ orbitalCorrectionCoeffs[18] * E * Math.sin(M + 2.0 * F)
+ orbitalCorrectionCoeffs[19] * Math.sin(O)
+ orbitalCorrectionCoeffs[20] * Math.sin(2.0 * M1 - 2.0 * F)
+ orbitalCorrectionCoeffs[21] * Math.sin(2.0 * M1 + 2.0 * F)
+ orbitalCorrectionCoeffs[22] * Math.sin(M1 - M - 2.0 * F)
+ orbitalCorrectionCoeffs[23] * Math.sin(M1 - M + 2.0 * F)
+ orbitalCorrectionCoeffs[24] * Math.sin(M1 + M - 2.0 * F)
+ orbitalCorrectionCoeffs[25] * Math.sin(M1 + M + 2.0 * F)
// this is the first approximation. all else is for style points! (47.1)
JDE = LUNATION_EPOCH_JD + (LUNATION_CYCLE_DURATION * k)
+ T * T * (0.0001337 + T * (-0.000000150 + 0.00000000073 * T))
JDE = JDE
+ orbitalApproximation
+ planetaryArguments
// extra correction for quarter phases
if (orbitalCorrectionIndex == 1) {
W = 0.00306
- 0.00038 * E * Math.cos(M)
+ 0.00026 * Math.cos(M1)
- 0.00002 * Math.cos(M1 - M)
+ 0.00002 * Math.cos(M1 + M)
+ 0.00002 * Math.cos(2.0 * F)
if (phi == 3){
W = -W
}
JDE += W
}
return JDE
}