Celestial Programming: Moon Phases

Celestial Programming : Celestial Programming: Moon Phases

This is an implementation of the algorithm from Meeus' Astronomical Algorithms for computing the dates of the phases of the Moon. View the source code of this page for the code, the code is public domain. The main algorithm is based on "lunation cycles", and since no one knows what those are, a helper function getCycleEstimate(year,month) is available. The value returned from that truly is an estimate, so the user should compute phases for cycles both before and after the cycle returned. In one of the examples in the book, Meeus uses the middle of the month as an estimate, and it might appear this is a sufficient method of getting all of the phases in a month, but it is not. Months are longer than lunation cycles, and if you always use the middle of the month as a starting estimate, you will eventually skip a lunation cycle. This page serves as an example of how to compute phases for an entire year, but the main algorithm is the getPhaseDate() function.

Year

New   2024 Dec 01 06:22:42
First 2024 Dec 08 15:27:59
Full  2024 Dec 15 09:02:53
Last  2024 Dec 22 22:19:27
New   2024 Dec 30 22:28:00
First 2025 Jan 06 23:57:33
Full  2025 Jan 13 22:28:05
Last  2025 Jan 21 20:31:55
New   2025 Jan 29 12:37:08
First 2025 Feb 05 08:03:19
Full  2025 Feb 12 13:54:34
Last  2025 Feb 20 17:33:48
New   2025 Feb 28 00:45:54
First 2025 Mar 06 16:32:47
Full  2025 Mar 14 06:55:46
Last  2025 Mar 22 11:30:45
New   2025 Mar 29 10:58:53
First 2025 Apr 05 02:15:48
Full  2025 Apr 13 00:23:26
Last  2025 Apr 21 01:36:46
New   2025 Apr 27 19:32:19
First 2025 May 04 13:52:55
Full  2025 May 12 16:57:11
Last  2025 May 20 11:59:54
New   2025 May 27 03:03:30
First 2025 Jun 03 03:42:11
Full  2025 Jun 11 07:45:00
Last  2025 Jun 18 19:20:13
New   2025 Jun 25 10:32:41
First 2025 Jul 02 19:31:22
Full  2025 Jul 10 20:37:55
Last  2025 Jul 18 00:38:48
New   2025 Jul 24 19:12:17
First 2025 Aug 01 12:42:26
Full  2025 Aug 09 07:56:14
Last  2025 Aug 16 05:13:26
New   2025 Aug 23 06:07:37
First 2025 Aug 31 06:26:14
Full  2025 Sep 07 18:10:00
Last  2025 Sep 14 10:34:10
New   2025 Sep 21 19:55:11
First 2025 Sep 29 23:54:51
Full  2025 Oct 07 03:48:42
Last  2025 Oct 13 18:13:52
New   2025 Oct 21 12:26:18
First 2025 Oct 29 16:21:59
Full  2025 Nov 05 13:20:28
Last  2025 Nov 12 05:29:18
New   2025 Nov 20 06:48:29
First 2025 Nov 28 07:00:02
Full  2025 Dec 04 23:15:15
Last  2025 Dec 11 20:52:53
New   2025 Dec 20 01:44:36
First 2025 Dec 27 19:11:06
Full  2026 Jan 03 10:04:05
Last  2026 Jan 10 15:49:38
New   2026 Jan 18 19:53:16
First 2026 Jan 26 04:48:42
Full  2026 Feb 01 22:10:28
Last  2026 Feb 09 12:44:24



    /*
    Greg Miller gmiller@gregmiller.net
    Algorithm from Meeus Astronomical Algorithms for computing dates of moon phases
    Released as public domain
    www.celestialprogramming.com
    */
    
    function computePhases(year){
        let k=getCycleEstimate(year,0);
        let text="";
    
        for(let i=-1;i<14;i++){
            let n=julainDateToGregorian(getPhaseDate(k+i,0));
            let f=julainDateToGregorian(getPhaseDate(k+i,.25));
            let h=julainDateToGregorian(getPhaseDate(k+i,.5));
            let l=julainDateToGregorian(getPhaseDate(k+i,.75));
    
            text+=format("New  ",n);
            text+=format("First",f);
            text+=format("Full ",h);
            text+=format("Last ",l);
        }
    
        document.getElementById("output").innerText=text;
    }
    
    function mod360(f){
        let t=f%360;
        if(t<0)t+=360;
        return t;
    }
    
    /*
    Year - integer year
    Month - integer month, January = 0
    
    Gets an estimate of the cycle number to be passed into getPhaseDate().  Since this
    is just an estimate, the user should also compute phase dates for cycles before and after
    the cycle number returned from this function
    */
    
    function getCycleEstimate(year,month){
        const yearfrac=(month*30+15)/365;  //Estimate fraction of year
        let k=12.3685*((year + yearfrac) - 2000);  //49.2
        k=Math.floor(k);
    
        return k;
    }
    
    /*
    Gets the given phase within a given cycle number based on the year 2000.  Compute cycle
    estimate using getCycleEstimate() above.
    
    cycle - integer cycle number from getCycleEstimate()
    phase - 0 = new, .25 = first quarter, .5 = full, .75 = last quarter, all other values are invalid
    
    returns JD of specified phase in TDB time scale.
    */
    
    //From Meeus ch49
    function getPhaseDate(cycle,phase){
        const k=cycle+phase;
    
        const toRad=Math.PI/180;
    
        const T = k/1236.85; //49.3
        let JDE = 2451550.09766 + 29.530588861*k + 0.00015437*T*T - 0.000000150*T*T*T + 0.00000000073*T*T*T*T; //49.1
    
        const E = 1 - 0.002516*T - 0.0000074*T*T; //47.6
    
        const M = mod360(2.5534 + 29.10535670*k - 0.0000014*T*T - 0.00000011*T*T*T)*toRad; //49.4
        const Mp = mod360(201.5643 + 385.81693528*k + 0.0107582*T*T + 0.00001238*T*T*T - 0.000000058*T*T*T*T)*toRad; //49.5
        const F = mod360(160.7108 + 390.67050284*k - 0.0016118*T*T - 0.00000227*T*T*T + 0.000000011*T*T*T*T)*toRad; //49.6
        const Om = mod360(124.7746 - 1.56375588*k + 0.0020672*T*T + 0.00000215*T*T*T)*toRad; //49.7
    
        //P351-352
        const A1 = mod360(299.77 + 0.107408*k - 0.009173*T*T)*toRad;
        const A2 = mod360(251.88 + 0.016321*k)*toRad;
        const A3 = mod360(251.83 + 26.651886*k)*toRad;
        const A4 = mod360(349.42 + 36.412478*k)*toRad;
        const A5 = mod360(84.66 + 18.206239*k)*toRad;
        const A6 = mod360(141.74 + 53.303771*k)*toRad;
        const A7 = mod360(207.14 + 2.453732*k)*toRad;
        const A8 = mod360(154.84 + 7.306860*k)*toRad;
        const A9 = mod360(34.52 + 27.261239*k)*toRad;
        const A10 = mod360(207.19 + 0.121824*k)*toRad;
        const A11 = mod360(291.34 + 1.844379*k)*toRad;
        const A12 = mod360(161.72 + 24.198154*k)*toRad;
        const A13 = mod360(239.56 + 25.513099*k)*toRad;
        const A14 = mod360(331.55 + 3.592518*k)*toRad;
    
        if (phase == 0) {
            correction = 0.00002*Math.sin(4*Mp) + -0.00002*Math.sin(3*Mp + M) + -0.00002*Math.sin(Mp - M - 2*F) + 0.00003*Math.sin(Mp - M + 2*F) + -0.00003*Math.sin(Mp + M + 2*F) +
                0.00003*Math.sin(2*Mp + 2*F) + 0.00003*Math.sin(Mp + M - 2*F) + 0.00004*Math.sin(3*M) + 0.00004*Math.sin(2*Mp - 2*F) + -0.00007*Math.sin(Mp + 2*M) + -0.00017*Math.sin(Om) +
                -0.00024*E*Math.sin(2*Mp - M) + 0.00038*E*Math.sin(M - 2*F) + 0.00042*E*Math.sin(M + 2*F) + -0.00042*Math.sin(3*Mp) + 0.00056*E*Math.sin(2*Mp + M) + -0.00057*Math.sin(Mp + 2*F) +
                -0.00111*Math.sin(Mp - 2*F) + 0.00208*E*E*Math.sin(2*M) + -0.00514*E*Math.sin(Mp + M) + 0.00739*E*Math.sin(Mp - M) + 0.01039*Math.sin(2*F) + 0.01608*Math.sin(2*Mp) +
                0.17241*E*Math.sin(M) + -0.40720*Math.sin(Mp);
        } else if ((phase == 0.25) || (phase == 0.75)) {
            correction = -0.00002*Math.sin(3*Mp + M) + 0.00002*Math.sin(Mp - M + 2*F) + 0.00002*Math.sin(2*Mp - 2*F) + 0.00003*Math.sin(3*M) + 0.00003*Math.sin(Mp + M - 2*F) + 0.00004*Math.sin(Mp - 2*M) +
                -0.00004*Math.sin(Mp + M + 2*F) + 0.00004*Math.sin(2*Mp + 2*F) + -0.00005*Math.sin(Mp - M - 2*F) + -0.00017*Math.sin(Om) + 0.00027*E*Math.sin(2*Mp + M) + -0.00028*E*E*Math.sin(Mp + 2*M) +
                0.00032*E*Math.sin(M - 2*F) + 0.00032*E*Math.sin(M + 2*F) + -0.00034*E*Math.sin(2*Mp - M) + -0.00040*Math.sin(3*Mp) + -0.00070*Math.sin(Mp + 2*F) + -0.00180*Math.sin(Mp - 2*F) +
                0.00204*E*E*Math.sin(2*M) + 0.00454*E*Math.sin(Mp - M) + 0.00804*Math.sin(2*F) + 0.00862*Math.sin(2*Mp) + -0.01183*E*Math.sin(Mp + M) + 0.17172*E*Math.sin(M) + -0.62801*Math.sin(Mp);    
    
            const W = 0.00306 - 0.00038*E*Math.cos(M) + 0.00026*Math.cos(Mp) - 0.00002*Math.cos(Mp - M) + 0.00002*Math.cos(Mp + M) + 0.00002*Math.cos(2*F);
            if (phase == 0.25){
                correction += W;
            } else {
                correction -= W;
            }
    
        } else if (phase == 0.5) {
            correction = 0.00002*Math.sin(4*Mp) + -0.00002*Math.sin(3*Mp + M) + -0.00002*Math.sin(Mp - M - 2*F) + 0.00003*Math.sin(Mp - M + 2*F) + -0.00003*Math.sin(Mp + M + 2*F) + 0.00003*Math.sin(2*Mp + 2*F) +
                0.00003*Math.sin(Mp + M - 2*F) + 0.00004*Math.sin(3*M) + 0.00004*Math.sin(2*Mp - 2*F) + -0.00007*Math.sin(Mp + 2*M) + -0.00017*Math.sin(Om) + -0.00024*E*Math.sin(2*Mp - M) +
                0.00038*E*Math.sin(M - 2*F) + 0.00042*E*Math.sin(M + 2*F) + -0.00042*Math.sin(3*Mp) + 0.00056*E*Math.sin(2*Mp + M) + -0.00057*Math.sin(Mp + 2*F) + -0.00111*Math.sin(Mp - 2*F) +
                0.00209*E*E*Math.sin(2*M) + -0.00514*E*Math.sin(Mp + M) + 0.00734*E*Math.sin(Mp - M) + 0.01043*Math.sin(2*F) + 0.01614*Math.sin(2*Mp) + 0.17302*E*Math.sin(M) + -0.40614*Math.sin(Mp);
        }
    
        JDE+=correction;
    
        //Additional corrections P 252
        correction = 0.000325*Math.sin(A1) + 0.000165*Math.sin(A2) + 0.000164*Math.sin(A3) + 0.000126*Math.sin(A4) + 0.000110*Math.sin(A5) + 0.000062*Math.sin(A6) + 0.000060*Math.sin(A7) +
            0.000056*Math.sin(A8) + 0.000047*Math.sin(A9) + 0.000042*Math.sin(A10) + 0.000040*Math.sin(A11) + 0.000037*Math.sin(A12) + 0.000035*Math.sin(A13) + 0.000023*Math.sin(A14);
        
        JDE += correction;
    
        return JDE;
    }