11
Vote

How to calculate a real moonphase

A small note about this script:
This is an original script from Dirk Malorny written in javascript !
I've translated it into LUA for all who wants to play with real moonphases (newmoon, halfmoon, fullmoon, age of moon ....)

So, have fun with this script

epoch     = 2444238.5

 

elonge    = 278.833540      -- Ecliptic longitude of the Sun at epoch 1980.0
elongp    = 282.596403      -- Ecliptic longitude of the Sun at perigee
eccent    = 0.016718        -- Eccentricity of Earth's orbit
sunsmax   = 1.49598022e8    -- Semi-major axis of Earth's orbit in km
sunangsiz = 0.533128        -- Sun's angular size, degrees, at semi-major axis distance

mmlong    = 64.975464       -- Moon's mean longitude at the epoch
mmlongp   = 349.383063      -- Mean longitude of the perigee at the epoch
mlnode    = 151.950429      -- Mean longitude of the node at the epoch
minc      = 5.145396        -- Inclination of the Moon's orbit
mecc      = 0.054900        -- Eccentricity of the Moon's orbit
mangsiz   = 0.5181          -- Moon's angular size at distance a from Earth
msmax     = 384401.0        -- Semi-major axis of Moon's orbit in km
mparallax = 0.9507          -- Parallax at distance a from Earth
synmonth  = 29.53058868     -- Synodic month (new Moon to new Moon)

earthrad  = 6378.16         -- Radius of Earth in kilometres

-- ***** math functions ***** --
function abs(x)
    if x < 0 then
        return x * -1
    else
        return x
    end
end

function torad(x)
    return x * (math.pi / 180)
end

function todeg(x)
    return x * (180 / math.pi)
end

function fixangle(x)
    return x - (360 * (math.floor(x / 360)))
end

function dsin(x)
    return math.sin(torad(x))
end

function dcos(x)
    return math.cos(torad(x))
end
function round(num, numDecimalPlaces)
    mult = 10 ^ (numDecimalPlaces or 0)

    return math.floor(num * mult + 0.5) / mult
end
-- ***** math functions END ***** --

function ucttoj(year, mon, mday, hour, min, sec)
    y = year
    m = mon
    if (m > 2) then
        m = m - 3
    else
        m = m + 9
        y = y - 1
    end

    c = y / 100
    y = y - (100 * c)

    jd = math.floor((mday + (c * 146097) / 4 + (y * 1461) / 4 + (m * 153 + 2) / 5 + 1721119)) - 0.5
    jt = (sec + 60 * (min + 60 * hour)) / 86400

    return jd + jt
end

function show(year, mon, mday, hour, min, sec)
    s = ''
    if (mday < 10) then s = s .. '0' end
    s = s .. mday .. '.'
    if (mon < 10) then s = s .. '0' end
    s = s .. mon .. '.' .. year .. ' - '
    if (hour < 10) then s = s .. '0' end
    s = s .. hour .. ':'
    if (min < 10) then s = s .. '0' end
    s = s .. min .. ':'
    if (sec < 10) then s = s .. '0' end
    s = s .. sec

 
   return s
end

function showjulian(td)

    td = td + 0.5

    j  = math.floor(td)

    j  = j - 1721119.0

    y  = math.floor(((4 * j) - 1) / 146097.0)

    j  = (j * 4.0) - (1.0 + (146097.0 * y))

    d  = math.floor(j / 4.0)

    j  = math.floor(((4.0 * d) + 3.0) / 1461.0)

    d  = ((4.0 * d) + 3.0) - (1461.0 * j)

    d  = math.floor((d + 4.0) / 4.0)

    m  = math.floor(((5.0 * d) - 3) / 153.0)

    d  = (5.0 * d) - (3.0 + (153.0 * m))

    d  = math.floor((d + 5.0) / 5.0)

    y  = (100.0 * y) + j

    

    if m < 10.0 then

       m = m + 3

    else

        m = m - 9

        y = y + 1

    end

    

    ij = ((td - math.floor(td)) * 86400.0)

    h  = math.floor(ij / 3600)

    mi = math.floor((ij / 60) % 60)

    s  = math.floor(ij % 60)

    

    return show(y, m, d, h, mi, s)

end

 

function meanphase(sdate, k)

    t   = (sdate - 2415020.0) / 36525

    t2  = t * t

    t3  = t2 * t

    nt3 = 2415020.75933 + synmonth * k

            + 0.0001178 * t2

            - 0.000000155 * t3

            + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2)

    

    return nt3

end

 

function kepler(m, ecc)

    m = torad(m)

    e = m

    epsilon = 1e-6

    

    repeat

        delta = e - ecc * math.sin(e) - m

        e = e - delta / (1 - ecc * math.cos(e))

    until abs(delta) > epsilon

    

    return e

end

 

function phase(pdate)

    day        = pdate - epoch

    n          = fixangle((360 / 365.2422) * day)

    m          = fixangle(n + elonge - elongp)

    ec         = kepler(m, eccent)

    ec         = math.sqrt((1 + eccent) / (1 - eccent)) * math.tan(ec / 2)

    ec         = 2 * todeg(math.atan(ec))

    lambdasun  = fixangle(ec + elongp)

    f          = ((1 + eccent * dcos(ec)) / (1 - eccent * eccent))

    sundist    = sunsmax / f

    sunang     = f * sunangsiz

    

    m1         = fixangle(13.1763966 * day * mmlong)

    mm         = fixangle(m1 - 0.1114041 * day * mmlong)

    mn         = fixangle(mlnode - 0.0529539 * day)

    ev         = 1.2739 * dsin(2 * (m1 - lambdasun) - mn)

    ae         = 0.1858 * dsin(m)

    a3         = 0.37 * dsin(m)

    mmp        = mm + ev - ae - a3

    mec        = 6.2886 * dsin(mmp)

    a4         = 0.214 * dsin(2 * mmp)

    lp         = m1 + ev + mec - ae +a4

    v          = 0.6583 * dsin(2 * (lp - lambdasun))

    lpp        = lp + v

    np         = mn - 0.16 * dsin(m)

    y          = dsin(lpp - np) * dcos(minc)

    x          = dcos(lpp - np)

    lambdamoon = todeg(math.atan2(y, x)) + np

    betam      = todeg(math.asin(dsin(lpp - np) * dsin(minc)))

    

    mage       = lpp - lambdasun

    moonphase  = (1 - dcos(mage)) / 360

    

    moondist  = (msmax * (1 - mecc * mecc)) / (1 + mecc * dcos(mmp + mec))

    moondfrac = moondist / msmax

    moonang   = mangsiz / moondfrac

    moonpar   = mparallax / moondfrac

    moonage   = synmonth * fixangle(mage) / 360

    

    return fixangle(mage) / 360

end

 

function truephase(k, phase)

    apcor = 0

    

    k      = k + phase

    t      = k / 1236.85

    t2     = t * t

    t3     = t2 * t

    pt     = 2415020.75933

            + synmonth * k

            + 0.0001178 * t2

            - 0.000000155 * t3

            + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2)

    m      = 359.2242

            + 29.10535608 * k

            - 0.0000333 * t2

            - 0.00000347 * t3

    mprime = 306.0253

            + 385.81691806 * k

            + 0.0107306 * t2

            + 0.00001236 * t3

    f = 21.2964

            + 390.67050646 * k

            - 0.0016528 * t2

            - 0.00000239 * t3;

    if phase < 0.01 or abs(phase - 0.5) < 0.01 then

        pt = pt + (0.1734 - 0.000393 * t) * dsin(m)

                + 0.0021 * dsin(2 * m)

                - 0.4068 * dsin(mprime)

                + 0.0161 * dsin(2 * mprime)

                - 0.0004 * dsin(3 * mprime)

                + 0.0104 * dsin(2 * f)

                - 0.0051 * dsin(m + mprime)

                - 0.0074 * dsin(m - mprime)

                + 0.0004 * dsin(2 * f + m)

                - 0.0004 * dsin(2 * f - m)

                - 0.0006 * dsin(2 * f + mprime)

                + 0.0010 * dsin(2 * f - mprime)

                + 0.0005 * dsin(m + 2 * mprime)

        apcor = 1

    elseif abs(phase - 0.25) < 0.01 or abs(phase - 0.75) < 0.01 then

         pt = pt + (0.1721 - 0.0004 * t) * dsin(m)

                + 0.0021 * dsin(2 * m)

                - 0.6280 * dsin(mprime)

                + 0.0089 * dsin(2 * mprime)

                - 0.0004 * dsin(3 * mprime)

                + 0.0079 * dsin(2 * f)

                - 0.0119 * dsin(m + mprime)

                - 0.0047 * dsin(m - mprime)

                + 0.0003 * dsin(2 * f + m)

                - 0.0004 * dsin(2 * f - m)

                - 0.0006 * dsin(2 * f + mprime)

                + 0.0021 * dsin(2 * f - mprime)

                + 0.0003 * dsin(m + 2 * mprime)

                + 0.0004 * dsin(m - 2 * mprime)

                - 0.0003 * dsin(2 * m + mprime);

        

        if phase < 0.5 then

            pt = pt + 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime)

        else

            pt = pt + -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime)

            apcor = 1

        end

    end

    

    return pt

end

 

function phasehunt(sdate)

    adate = sdate - 45

    td    = adate + 0.5

    j     = math.floor(td)

    j     = j - 1721119.0

    yy    = math.floor(((4 * j) - 1) / 146097.0)

    j     = (j * 4.0) - (1.0 + (146097.0 * yy))

    dd    = math.floor(j / 4.0)

    j     = math.floor(((4.0 * dd) + 3.0) / 1461.0)

    dd    = ((4.0 * dd) + 3.0) - (1461.0 * j)

    dd    = math.floor((dd + 4.0) / 4.0)

    mm    = math.floor(((5.0 * dd) - 3) / 153.0)

    dd    = (5.0 * dd) - (3.0 + (153.0 * mm))

    dd    = math.floor((dd + 5.0) / 5.0)

    yy    = (100.0 * yy) + j

 

    if (mm < 10.0) then

        mm = mm + 3

    else

        mm = mm - 9

        yy = yy + 1

    end

    

    k1 = math.floor((yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685)

    nt1 = meanphase(adate, k1)

    adate = nt1

    

    found = 0

    while found == 0 do

        adate = adate + synmonth

        k2    = k1 + 1

        nt2   = meanphase(adate, k2)

        if nt1 <= sdate and nt2 > sdate then

            found = 1

        else

            nt1 = nt2

            k1  = k2

       end

    end

    

    phases0 = truephase(k1, 0.0)

    phases1 = truephase(k1, 0.25)

    phases2 = truephase(k1, 0.5)

    phases3 = truephase(k1, 0.75)

    phases4 = truephase(k1 + 1, 0.0)

end

 

function moon_phase(y,m,d,h,mi,s)

    year    = y

    month   = m

    day     = d

    hours   = h

    minutes = mi

    seconds = s

    if year < 1000 then year = year + 1904 end

    

    if hours > 23 then

        hours = hours - 24

    end

    if hours < 0 then

        hours = hours + 24

        day = day - 1

    end

    

    julian = ucttoj(year, month, day, hours, minutes, seconds)

    img = round(99 * phase(julian) + 1)

    

    mad = math.floor(moonage)

    mah = math.floor(24 * (moonage - mad))

    mam = math.floor(60 * (24 * (moonage - mad) - mah))

    mas = math.floor(60 * (60 * (24 * (moonage - mad) - mah) - mam))

    

    phasehunt(julian + 0.5)

    

    -- give out some text (use this vars in Watchmaker)

    var_ms_t1 = showjulian(phases0)     -- last newmoon

    var_ms_t2 = showjulian(phases1)     -- first quater

    var_ms_t3 = showjulian(phases2)     -- fullmoon

    var_ms_t4 = showjulian(phases3)     -- last quater

    var_ms_t5 = showjulian(phases4)     -- next newmoon

end

 

-- start script here and change year, month, day, hour, minut, second
moon_phase(2018, 2, 1, 0, 0, 0)

Published by Sabine Mayerhagen on 22 April 2017WatchMaker Tips & Tricks Posted

This Website and It's Resources are Powered by:
  
WatchMaker Beautiful Watches for Android Wear

No smartwatch required!! Best clock wallpaper for Android!
  
WatchMaker Live WallPaper

Introducing the Perfect Companion App for WatchMaker...
  
GearsPro (OBD 2 & Car)