Header

/////////////////////////////////////////////////////////////////////////
///
/// moon.h
///
///
/// Created March 31, 2004 by Don Dugger
///
///
///
/////////////////////////////////////////////////////////////////////////
///
/// <PRE>
/// THE SOFTWARE IS PROVIDED ~AS IS~, WITHOUT WARRANTY OF ANY KIND,
/// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
/// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
/// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
/// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
/// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
/// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
/// DEALINGS IN THE SOFTWARE.
/// </PRE>
///
/////////////////////////////////////////////////////////////////////////

#ifndef _moon_h_
#define _moon_h_

/////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <math.h>

/////////////////////////////////////////////////////////////////////////
/// <STRONG>
/// Classes
/// </STRONG>
/// <p>
///
///   This class is an all static.
/// </p><p>
///   This group of methods deals with the phases of the moon, it also computes<BR>
///   Julian dates, due to the fact it was needed for the moon calculations. To make it<BR>
///   simpler to work with the Julian calculation there left in rather than creating a<BR>
///   separate class.
/// </p><p>
///   Some functions are taken from "Numerical Recipes in C - The Art of Scientific Computing"
///   Second Edition, translated with changes for C++ and help clarify the code.
/// <BR>
///   ISBN 0-521-43108-5
/// <BR>
///   (Copyright &copy; Cambridge University Press 1988, 1992)
/// </p><p>
///   \section lunarphase  Lunar Phases
/// </p><p>
///    The program defines the phase of the moon as ia floating point number between 0 and 4.
///    The following values apply:
///    - 0.0 = New Moon
///    - 1.0 = 1st Quarter
///    - 2.0 = Full Moon
///    - 3.0 = Last Quarter
/// </p>
///@brief Calculate the phase of the moon.
///
/////////////////////////////////////////////////////////////////////////
class Moon {
    /// Radians to degree conversion
    static const double RAD;
    /// We will call this the Gregorian Interval.
    /// The Gregorian Calendar was adopted in October 15, 1582.
    static const long IGREG;
    // The mean lunar cycle
    static const double MEAN_LUNAR_CYCLE;

public:

    Moon() {};
    ~Moon() {};

    /////////////////////////////////////////////////////////////////////////
    /// phase() calculates the phase of the moon at
    /// noon of the date given.
    ///<p> <i>See above</i> @ref lunarphase
    ///@param mon The month.
    ///@param day The day of the month.
    ///@param year The year.
    ///@return  The Lunar Phase <i>see above</i> \ref lunarphase
    static double phase(int mon,int day,int year);

    /////////////////////////////////////////////////////////////////////////
    /// nextphase() calculates the date & time of next phase of the moon,
    /// that is the next hole number, given the start date.
    ///@param mon   The month.
    ///@param day   The day of the month.
    ///@param year  The year.
    ///@param hr    The hour.
    ///@param min   The minute.
    ///@param nph   The Lunar Phase <i>See above</i> \ref lunarphase
    ///@return  Note that all the return values are references.
    static void nextphase(int& mon,int& day,int& year,int& hr,int& min,int nph);

    /////////////////////////////////////////////////////////////////////////
    /// flmoon() calculates the Julian date of nth lunar cycle and nph phase
    /// of the moon since Jan 1900.
    /// n is the number of lunar cycle since Jan 1900.
    /// This is the original function from the book.
    /// It's been modified to work in C++.
    /// I also rewrote it to be little clearer. And added
    /// comments where I understood what was going on and 
    /// increased the accuracy of a few constants.
    ///@param n     The number of lunar cycles.
    ///@param nph   The lunar phase.
    ///@param jd    The Julian date number.
    ///@param frac  The fractional part of the Julian date number.
    ///@return      The full Julian date of the nth lunar cycle plus nph phase
    ///             Note that the integer and fractional parts are
    ///             returned by references.
    static double flmoon(int n,int nph,long& jd,double& frac);

    /////////////////////////////////////////////////////////////////////////
    /// Julian2Calendar() computes the Julian date
    ///@param jd    The Julian date number.
    ///@param mon   The month.
    ///@param day   The day of the month.
    ///@param year  The year.
    ///@return  Note that all the return values are references.
    static void Julian2Calendar(long jd,int& mon,int& day,int& year);

    /////////////////////////////////////////////////////////////////////////
    /// Calendar2Julian() computes the date given a Julian date
    ///@param mon   The month.
    ///@param day   The day of the month.
    ///@param year  The year.
    static long Calendar2Julian(int mon,int day,int year);

    /////////////////////////////////////////////////////////////////////////
    /// isDST() returns true if the given date is day light savings time
    ///     There are two versions, one which takes the Julian date and
    ///     one that takes the month, day and year.
    ///@param   jd  The Julian date
    static bool isDST(long jd);
    /////////////////////////////////////////////////////////////////////////
    ///@param mon   The month.
    ///@param day   The day of the month.
    ///@param year  The year.
    static bool isDST(int mon,int day,int year);
    
};


/////////////////////////////////////////////////////////////////////////

#endif /// _moon_h_

                

Implementation

///////////////////////////////////////////////////////////////
///
///  moon.cpp
/// 
/// 
///  Created March 31, 2004 by Don Dugger
/// 
///
///////////////////////////////////////////////////////////////
///
/// THE SOFTWARE IS PROVIDED ~AS IS~, WITHOUT WARRANTY OF ANY KIND,
/// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
/// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
/// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
/// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
/// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
/// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
/// DEALINGS IN THE SOFTWARE.
///
///////////////////////////////////////////////////////////////

#include "moon.h"

///////////////////////////////////////////////////////////////
const double Moon::RAD = (3.1415926535897/180.0);
const long   Moon::IGREG = (15+31L*(10+12L*1582));
const double Moon::MEAN_LUNAR_CYCLE = 29.53058868;
///////////////////////////////////////////////////////////////

double Moon::phase(int mon,int day,int year)
{
    // the new moon Julian date
    long    nm_jd = 0;
    double  nm_frac = 0.0;
    // the 1st quarter Julian date
    long    fq_jd = 0;
    double  fq_frac = 0.0;
    // the full moon Julian date
    long    fm_jd = 0;
    double  fm_frac = 0.0;
    // the last quarter Julian date
    long    lq_jd = 0;
    double  lq_frac = 0.0;
    // The Julian date entered
    long    jd = 0;
    double  frac = 0.0;
    // The lunar phase
    double  ph = 0.0;
    // The portion of a lunar phase which elapse per day
    double phase_per_day = 1/(MEAN_LUNAR_CYCLE/4);
    // 12.37 is the number of lunar cycles per year
    // n then is the number of lunar cycles since 1900 and the year month entered
    int n = (int)(12.37*((year-1900)+((mon-0.5)/12.0)));

    // Get the Julian Date and make sure date entered is valid
    if(!(jd = Calendar2Julian(mon,day,year))) return(-1.0);
    // Now get the first date of the new moon of year mouth entered
    flmoon(n,0,nm_jd,nm_frac);
    // If it after the entered date backup one lunar cycle
    if(jd < nm_jd) flmoon(--n,0,nm_jd,nm_frac);
    // Now find the 1st quarter
    flmoon(n,1,fq_jd,fq_frac);
    // If it's after or on the day entered 
    if(fq_jd >= jd) {
        // Now calculate the phase which we now know is between the
        // new moon and first quarter 
        if(abs((int)jd - nm_jd) < abs((int)jd - fq_jd)) {
            // We are closer to the new moon so no phase offset needed 
            ph = (double)(((double)(jd) - ((double)nm_jd+nm_frac)) * phase_per_day);
            if(ph < 0) ph = 4.0 + ph;
            return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4
        } else {
            // we are closer to the 1st quarter so add it in
            ph = 1.0+(double)(((double)(jd) - ((double)fq_jd+fq_frac)) * phase_per_day);
            return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4
        }
    }
    // Now find the full moon quarter
    flmoon(n,2,fm_jd,fm_frac);
    // If it's after or on the day entered 
    if(fm_jd >= jd) {
        // Now calculate the phase which we now know is between the
        // first quarter and the full moon 
        if(abs((int)jd - fq_jd) < abs((int)jd - fm_jd)) {
            // we are closer to the 1st quarter so add it in
            ph = 1.0+(double)(((double)(jd) - ((double)fq_jd+fq_frac)) * phase_per_day);
            return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4
        } else {
            // we are closer to the full moon so add it in
            ph = 2.0+(double)(((double)(jd) - ((double)fm_jd+fm_frac)) * phase_per_day);
            return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4
        }
    }
    // Now find the last quarter
    flmoon(n,3,lq_jd,lq_frac);
    // If it's after or on the day entered 
    if(lq_jd >= jd) {
        // Now calculate the phase which we now know is between the
        // full moon and last quarter
        if(abs((int)jd - fm_jd) < abs((int)jd - lq_jd)) {
            // we are closer to the full moon so add it in
            ph = (2.0+(double)(((double)(jd) - ((double)fm_jd+fm_frac)) * phase_per_day));
            return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4
        } else {
            // we are closer to the last quarter so add it in
            ph = (3.0+(double)(((double)(jd) - ((double)lq_jd+lq_frac)) * phase_per_day));
            return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4
        }
    }
    // Now find the new moon quarter
    flmoon(++n,0,nm_jd,nm_frac);
    // Now calculate the phase which we now know is after the last quarter
    ph = (4.0+(double)(((double)(jd) - ((double)nm_jd+nm_frac)) * phase_per_day));
    return(ph > 4.0 ? ph - 4.0 : ph); // Keep the value between 0-4

} // End of phase()

///////////////////////////////////////////////////////////////

void Moon::nextphase(int& mon,int& day,int& year,int& hr,int& min,int nph)
{
    long jd = 0;
    double djd = 0.0;
    long ph_jd = 0;
    double frac = 0.0;
    int n = (int)(12.37*((year-1900)+((mon-0.5)/12.0)));
    double intpart = 0.0;

    jd = Calendar2Julian(mon,day,year);
    djd = flmoon(n,nph,ph_jd,frac) + 0.5;
    while(jd < ph_jd) {
        djd = flmoon(--n,nph,ph_jd,frac) + 0.5;
    }
    while(jd > ph_jd) {
        djd = flmoon(++n,nph,ph_jd,frac) + 0.5;
    }

    struct tm *tmpt;
    time_t e;

    e = time(NULL);
    tmpt = localtime(&e);
    djd += ((double)tmpt->tm_gmtoff * (1.0/86400.0));
    if(tmpt->tm_isdst) djd -= (1.0/24.0);
    ph_jd = (long)floor(djd);
    if(isDST(jd)) {
        djd += (1.0/24.0);
        ph_jd = (long)floor(djd);
    }
    frac = djd - floor(djd);
    
    Julian2Calendar(ph_jd,mon,day,year);
    frac *= 24.0;
    hr = (int)(frac >= 0.0 ? floor(frac) : ceil(frac-1.0));
    frac -= (double)hr;
    min = (int)floor(60*frac);

} // End of nextphase()

///////////////////////////////////////////////////////////////

double Moon::flmoon(int n,int nph,long& jd,double& frac)
{
    int int_part = 0;
    double lunar_cycles = n + ( nph / 4.0 ); // the total number lunar cycle
                                             // nph = 4 is one lunar cycle 
    double t = lunar_cycles / 1236.85;       // 1236.85 is the number of lunar cycle per Julian Century
    double t2 = t * t;                       // Square for frequent use

    // Sun's mean anomaly
    double sun_anomaly = 359.2242 + ( 29.10535608 * lunar_cycles );

    // Moon's mean anomaly
    double moon_anomaly = 306.0253 + ( 385.81691806 * lunar_cycles ) + ( 0.010730 * t2 );

    // Not sure of what's up here, but two of the numbers very interesting
    // Notice that 2415020.75933 is used as a reference epochs for mean time of lunar cycle which is 1900 Jan 1 at 6:13 AM
    // And 29.53058868 is the "Synodic month" or the mean time in days between new moons.
    double xtra = 0.75933 + ( 1.53058868 * lunar_cycles ) + ( ( ( 1.78e-4 ) - ( 1.55e-7 ) * t ) * t2 );

    // Ok 2415020 is Julian date 1899 at noon + .75933 is the Julian date of a new moon
    // 28 + 1.53058868 is mean time of a lunar cycle ,
    // and 7 close to 7.38264717 a lunar phase but unlike
    // the first two the fractional part dose not seem to be accounted for. 
    jd = 2415020 + ( 28L * n ) + ( 7L * nph );

    // Looks like this is all being done to adjust the variations in the mean lunar cycle of 29.53058868 
    // Which from what I understand is due the moons orbit in relationship to the earths orbit around the sun. 
    if(nph == 0 || nph == 2) // New or Full - sun earth moon inline
        xtra += ( ( 0.1734-3.93e-4 * t ) * sin(RAD*sun_anomaly) ) - ( 0.4068 * sin(RAD*moon_anomaly) );
    else // 1st quarter last quarter - moon 90 degrees left or right of the earth and suns line 
        xtra += ( ( 0.1721-4.0e-4 * t ) * sin(RAD*sun_anomaly) ) - ( 0.6280 * sin(RAD*moon_anomaly) );

    // This parts easy just put the integer and fractional parts right
    int_part = (int)( xtra >= 0.0 ? floor(xtra) : ceil(xtra-1.0) );
    jd += int_part;
    frac = xtra - int_part;
    return double(jd + frac);

} // End of flmoon()


///////////////////////////////////////////////////////////////

void Moon::Julian2Calendar(long jd,int& mon,int& day,int& year)
{
    const long GREG = 2299161L;
    long a = 0L;
    long b = 0L;
    long c = 0L;
    long d = 0L;
    long e = 0L;

    if(jd >= GREG) {
        a = (long)((((jd-1867216)-0.25)/36524.25));
        a = (long)(jd+1+(long)a-(long)(0.25*(long)a));
    } else {
        a = jd;
    }
    b = a+1524;
    c = (long)(6680+((b-2439870)-122.1)/365.25);
    d = 365*c+(long)(0.25*c);
    e = (long)((b-d)/30.6001);
    day = (int)(b-d-(long)(30.6001*e));
    mon = (int)(e-1);
    if(mon > 12) mon -= 12;
    year = (int)(c-4715);
    if(mon > 2) --year;
    if(year <= 0) --year;

} // End of Julian2Calendar()

///////////////////////////////////////////////////////////////

long Moon::Calendar2Julian(int mon,int day,int year)
{
    long jul;
    int ja,jy=year,jm;

    if(jy == 0) return(0);
    if(jy < 0) ++jy;
    if(mon > 2) jm=mon+1;
    else {
        --jy;
        jm=mon+13;
    }
    jul = (long)(floor(365.25*jy)+floor(30.6001*jm)+day+1720995);
    if(day+31L*(mon+12L*year) >= IGREG) {
        ja=(int)(0.01*jy);
        jul += 2-ja+(int)(0.25*ja);
    }
    return(jul);

} // End of Calendar2Julian()

///////////////////////////////////////////////////////////////

bool Moon::isDST(long jd)
{
    int m,d,y;

    Julian2Calendar(jd,m,d,y);
    return(isDST(m,d,y));

} // End of isDST()

///////////////////////////////////////////////////////////////

bool Moon::isDST(int mon,int day,int year)
{
    struct tm tmp;
    struct tm *tmpt;
    time_t e;

    e = time(NULL);
    tmpt = localtime(&e);
    tmp.tm_gmtoff = tmpt->tm_gmtoff;
    tmp.tm_mon = mon - 1;
    tmp.tm_mday = day;
    tmp.tm_year = year - 1900;
    tmp.tm_sec = 0;
    tmp.tm_min = 0;
    tmp.tm_hour = 3;
    tmp.tm_isdst = -1;
    e = mktime(&tmp);
    tmpt = localtime(&e);
    if(tmpt->tm_isdst) return(true);
    return(false);

} // End of isDST()

///////////////////////////////////////////////////////////////
///////////////////// END OF FILE /////////////////////////////
///////////////////////////////////////////////////////////////
                

phase() Example

Using clang++:

clang++ -g -std=c++2b -c phase.cpp -o phase.o -L./ -lmoon
                

Using g++:

g++ -g -c phase.cpp -o phase.o -L./ -lmoon
                

To execute:

phase 2020-01-12
                

///////////////////////////////////////////////////////////////////////
// 
//  phase.cpp
// 
//  Written by: Don Dugger
//        Date: 2025-02-17
//
//  Copyright (C) 2025 Acme Software Works
//
///////////////////////////////////////////////////////////////////////
///
/// <PRE>
/// THE SOFTWARE IS PROVIDED ~AS IS~, WITHOUT WARRANTY OF ANY KIND,
/// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
/// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
/// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
/// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
/// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
/// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
/// DEALINGS IN THE SOFTWARE.
/// </PRE>
///
///////////////////////////////////////////////////////////////////////
 
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string>
#include <iostream>
#include <limits.h>
#include <vector>
#include "moon.h"
 
///////////////////////////////////////////////////////////////////////
using namespace std;
///////////////////////////////////////////////////////////////////////
//
// Split a string at the delimiter
//
vector<string> split(vector<string>& list, string str, string delimiter) 
{
    size_t start = 0;
    size_t end;
    size_t size = delimiter.size();
    string segment;

    while ((end = str.find(delimiter, start)) != string::npos) { // look for the delimiter
        segment = str.substr(start, end - start); // Extract the segment
        list.push_back(segment);                  // Save the segment
        start = end + size;                       // Skip the delimiter to the nex segment
    }
    // Save the last segment
    list.push_back(str.substr(start));
    return list;

} // End of split ()


///////////////////////////////////////////////////////////////////////
//
// main
//
int main(int argc, char *argv[])
{
    // Make sure the date is on the command line.
    if ( argc != 2 ) {
        cerr << "Usage: phase yyyy-mm-dd" << endl;
        exit(1);
    }

    // A vector to hold the parts of the date
    vector<string> date_parts;
    // Split the date at the character "-"
    date_parts = split(date_parts,argv[1],"-");

    // Make sure date was close to the correct format
    if ( date_parts.size() != 3 ) {
        cerr << "Usage: phase yyyy-mm-dd" << endl;
        exit(1);
    }

    // Variables to hold the date
    int year, month, day;

    // Convert the parts to integers 
    year  = (int)stol(date_parts[0].c_str());
    month = (int)stol(date_parts[1].c_str());
    day   = (int)stol(date_parts[2].c_str());
    
    // And output the phase
    cout << Moon::phase(month,day,year) << endl;

    return 0;

} // End of main()

///////////////////////////////////////////////////////////////////////
//////////////////////// End of File //////////////////////////////////
///////////////////////////////////////////////////////////////////////

                

nextphase() Example

Using clang++:

clang++ -g -std=c++2b -c nextphase.cpp -o nextphase.o -L./ -lmoon
                

Using clang++:

g++ -g -c nextphase.cpp -o nextphase.o -L./ -lmoon
                

To execute: (the tailing 2 is for a full moon)

nextphase 2020-01-12 2 
                

///////////////////////////////////////////////////////////////////////
// 
//  phase.cpp
// 
//  Written by: Don Dugger
//        Date: 2025-02-17
//
//  Copyright (C) 2025 Acme Software Works
//
///////////////////////////////////////////////////////////////////////
///
/// <PRE>
/// THE SOFTWARE IS PROVIDED ~AS IS~, WITHOUT WARRANTY OF ANY KIND,
/// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
/// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
/// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
/// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
/// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
/// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
/// DEALINGS IN THE SOFTWARE.
/// </PRE>
///
///////////////////////////////////////////////////////////////////////
 
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string>
#include <iostream>
#include <limits.h>
#include <vector>
#include <iomanip>
#include "moon.h"
 
///////////////////////////////////////////////////////////////////////
using namespace std;
///////////////////////////////////////////////////////////////////////
//
// Split a string at the delimiter
//
vector<string> split(vector<string>& list, string str, string delimiter) 
{
    size_t start = 0;
    size_t end;
    size_t size = delimiter.size();
    string segment;

    while ((end = str.find(delimiter, start)) != string::npos) { // look for the delimiter
        segment = str.substr(start, end - start); // Extract the segment
        list.push_back(segment);                  // Save the segment
        start = end + size;                       // Skip the delimiter to the nex segment
    }
    // Save the last segment
    list.push_back(str.substr(start));
    return list;

} // End of split ()


///////////////////////////////////////////////////////////////////////
//
// main
//
int main(int argc, char *argv[])
{

    // A vector to hold the parts of the date
    vector<string> date_parts;
    // Split the date at the character "-"
    date_parts = split(date_parts,argv[1],"-");
    // The phase and convert it to an integer
    int phase = (int)strtol(argv[2],NULL,10);

    // Variables to hold the date
    int year, month, day;
    // Set the time to midnight
    int hour = 0;
    int minute = 0;

    // Convert the parts to integers 
    year  = (int)strtol(date_parts[0].c_str(),NULL,10);
    month = (int)strtol(date_parts[1].c_str(),NULL,10);
    day   = (int)strtol(date_parts[2].c_str(),NULL,10);
    
    // Get the date of the next phase
    //   Note it is return by modifying the date values
    //   it was handed.
    Moon::nextphase(month,day,year,hour,minute,phase);

    // And output the date of the phase
    cout << year << "-" << setw(2) << setfill('0') << month;
    cout << "-" << setw(2) <<  setfill('0') << day << endl;

    return 0;

} // End of main()

///////////////////////////////////////////////////////////////////////
//////////////////////// End of File //////////////////////////////////
///////////////////////////////////////////////////////////////////////

                

///////////////////////////////////////////////////////////////////////////////////////////////
//
//  moonphase.js
// 
// Written by: Don Dugger
//       Date: 2024-04-29
//
// Copyright (C) 2024 Acme Software Works, Inc.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// 
// 1. Redistributions of source code must retain the above copyright notice,
//    this list of conditions and the following disclaimer.
// 
// 2. Redistributions in binary form must reproduce the above copyright notice,
//    this list of conditions and the following disclaimer in the documentation
//    and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
// GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
// OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
///////////////////////////////////////////////////////////////////////////////////////////////

class MoonPhase {

    ///////////////////////////////////////////////////////////////////////////////////////////////
    RAD = (3.1415926535897/180.0);
    IGREG = (15+31*(10+12*1582));
    MEAN_LUNAR_CYCLE = 29.53058868;
    ///////////////////////////////////////////////////////////////////////////////////////////////

    flmoon(n,nph)
    {
        var jd = 0;
        var int_part = 0;
        var lunar_cycles = n + ( nph / 4.0 ); // the total number lunar cycle
                                              // nph = 4 is one lunar cycle 
        var t = lunar_cycles / 1236.85;       // 1236.85 is the number of lunar cycle per Julian Century
        var t2 = t * t;                       // Square for frequent use

        // Sun's mean anomaly
        var sun_anomaly = 359.2242 + ( 29.10535608 * lunar_cycles );

        // Moon's mean anomaly
        var moon_anomaly = 306.0253 + ( 385.81691806 * lunar_cycles ) + ( 0.010730 * t2 );

        // Not sure of what's up here, but two of the numbers very interesting
        // Notice that 2415020.75933 is used a referce epochs for mean time of lunar cycle which is 1900 Jan 1 at 6:13 AM
        // And 29.53058868 is the "Synodic month" or the mean time in days between new moons.
        var xtra = 0.75933 + ( 1.53058868 * lunar_cycles ) + ( ( ( 1.78e-4 ) - ( 1.55e-7 ) * t ) * t2 );

        // Ok 2415020 is Julian date 1899 at noon + .75933 is the Julian date of a new moon
        // 28 + 1.53058868 is mean time of a lunar cycle ,
        // and 7 close to 7.38264717 a lunar phase but unlike
        // the first two the tractional part dose not seem to be accounted for. 
        jd = 2415020 + ( 28 * n ) + ( 7 * nph );

        // Looks like this is all being done to adjust the variations in the maen lunar cycle of 29.53058868 
        // Which from what I understand is due the moons orbit in relationship to the earths orbit around the sun. 
        if(nph == 0 || nph == 2) // New or Full - sun earth moon inline
            xtra += ( ( 0.1734-3.93e-4 * t ) * Math.sin(this.RAD*sun_anomaly) ) - ( 0.4068 * Math.sin(this.RAD*moon_anomaly) );
        else // 1st quarder last quarter - moon 90 deegree left or right of the earth and suns line 
            xtra += ( ( 0.1721-4.0e-4 * t ) * Math.sin(this.RAD*sun_anomaly) ) - ( 0.6280 * Math.sin(this.RAD*moon_anomaly) );

        // This parts easy just put the integer and fractional parts right
        int_part = parseInt( xtra >= 0.0 ? Math.floor(xtra) : Math.ceil(xtra-1.0) );
        jd += int_part;
        var frac = xtra - int_part;
        return (jd + frac);

    } // End of flmoon


    ///////////////////////////////////////////////////////////////

    nextphase(year,mon,day,hr,min,nph)
    {
        var jd = 0;
        var djd = 0.0;
        var ph_jd = 0;
        var frac = 0.0;
        // 12.37 is the number of lunar cycles per year
        // n then is the number of lunar cycles since 1900 and the year month entered
        var n = parseInt(12.37*((year-1900)+((mon-0.5)/12.0)));
        var intpart = 0.0;

        // Get the starting Julian date 
        jd = this.Calendar2Julian(year,mon,day);
        djd = this.flmoon(n,nph) + 0.5;
        while(jd < djd) {
            djd = this.flmoon(--n,nph) + 0.5;
        }
        while(jd > djd) {
            djd = this.flmoon(++n,nph) + 0.5;
        }


        var now = new Date();
        djd -= (now.getTimezoneOffset() * (1.0/1440.0));
        if(this.isDST(now)) djd -= (1.0/24.0);
        ph_jd = parseInt(Math.floor(djd));
        if(this.isDST(this.Julian2Calendar(jd))) {
            djd += (1.0/24.0);
            ph_jd = parseInt(Math.floor(djd));
        }
        frac = djd - Math.floor(djd);
        
        var date = this.Julian2Calendar(ph_jd);
        frac *= 24.0;
        hr = parseInt(frac >= 0.0 ? Math.floor(frac) : Math.ceil(frac-1.0));
        frac -= hr;
        min = parseInt(Math.floor(60*frac));
        date.setHours(hr);
        date.setMinutes(min);
        return date;

    } // End of nextphase()

    ///////////////////////////////////////////////////////////////
    //
    isDST(date) {
        let jan = new Date(date.getFullYear(), 0, 1).getTimezoneOffset();
        let jul = new Date(date.getFullYear(), 6, 1).getTimezoneOffset();
        return Math.max(jan, jul) !== date.getTimezoneOffset();    
    }

    ///////////////////////////////////////////////////////////////

    Julian2Calendar(jd)
    {
        const GREG = 2299161;
        var a = 0;
        var b = 0;
        var c = 0;
        var d = 0;
        var e = 0;
        var year = 0;
        var month = 0;
        var day = 0;

        if(jd >= GREG) {
            a = parseInt((((jd-1867216)-0.25)/36524.25));
            a = parseInt(jd+1+parseInt(a)-parseInt(0.25*parseInt(a)));
        } else {
            a = jd;
        }
        b = a+1524;
        c = parseInt(6680+((b-2439870)-122.1)/365.25);
        d = 365*c+parseInt(0.25*c);
        e = parseInt((b-d)/30.6001);
        day = parseInt((b-d-parseInt(30.6001*e)));
        month = parseInt(e-1);
        if(month > 12) month -= 12;
        year = parseInt(c-4715);
        if(month > 2) --year;
        if(year <= 0) --year;
        return new Date(year,month-1,day);

    } // End of Julian2Calendar

    ///////////////////////////////////////////////////////////////

    Calendar2Julian(year,month,day)
    {
        var jul;
        var ja,jy=year,jm;

        if(jy == 0) return(0);
        if(jy < 0) ++jy;
        if(month > 2) jm=month+1;
        else {
            --jy;
            jm=month+13;
        }
        jul = parseInt(Math.floor(365.25*jy)+Math.floor(30.6001*jm)+day+1720995);
        if(day+31*(month+12*year) >= this.IGREG) {
            ja=parseInt(0.01*jy);
            jul += 2-ja+parseInt(0.25*ja);
        }
        return(jul);

    } // End of Calendar2Julian

} // End of class MoonPhase

///////////////////////////////////////////////////////////////
///////////////////// END OF FILE /////////////////////////////
///////////////////////////////////////////////////////////////
                

clang++ -g -std=c++2b -c moon.cpp -o moon.o
ar -r libmoon.a moon.o

g++ -g -c moon.cpp -o moon.o
ar -r libmoon.a moon.o

The Methods

static double phase(int mon,int day,int year);

static void nextphase(int& mon,int& day,int& year,int& hr,int& min,int nph);

static double flmoon(int n,int nph,long& jd,double& frac);

static void Julian2Calendar(long jd,int& mon,int& day,int& year);

static long Calendar2Julian(int mon,int day,int year);

static bool isDST(long jd);

static bool isDST(int mon,int day,int year);