Saturday, November 18, 2017

#1 2013-10-11 04:56:19 am

McUsrII
Member
Registered: 2012-11-21
Posts: 3046
Website

Accurate Table of Moonphases

Hello.

The script below produces accurate tables of Moon Phases, besides tell you the current phase of the Moon.

Enjoy! smile

Edit

You need Satimage.osax in order to run this!

Applescript:

-- © 2013 McUsr, sustainable share of code kindly provided by Nigel Garvey.
-- It is in public domain Licenced under GLPGL 3.0, which means
-- that if you use this then you must provide the source code,
-- and licence your product too under GLPGL 3.0.

script MoonPhases
   property scriptTitle : "Moon Phases"
   
   property JD0 : («data isot303030312D31312D3234543132» as date) ¬
       - ((«data isot343731342D31312D3234543132» as date) - ¬
       ((«data isot303030312D31312D3234543132» as date) ¬
           - 366 * days))
   -- Nigel Garvey [url]http://macscripter.net/viewtopic.php?pid=167391#p167391[/url] post # 6
   
   property Debug : true
   property degrad : pi / 180
   
   property movingStates : {¬
       "Waning Crescent", ¬
       "Waxing Crescent", ¬
       "Waxing Gibbous", ¬
       "Waning Gibbous"}
   
   property exactStates : {¬
       "New Moon", ¬
       "First Quarter", ¬
       "Full Moon", ¬
       "Third Quarter"}
   
   property aNewMoon : 1
   property tFirstQuarter : 2
   property aFullMoon : 3
   property tThirdQuarter : 4
   
   property monthNames : {¬
       "January", ¬
       "February", ¬
       "Mars", ¬
       "April", ¬
       "May", ¬
       "June", ¬
       "July", ¬
       "August", ¬
       "September", ¬
       "October", ¬
       "November", ¬
       "December"}
   
   
   on createMoonPhasesTable()
       local dyr, startY, endY, MoonPhases
       
       set {startY, endY} to getRangeOfYears((current date)'s year, "Decade", scriptTitle)
       if endY - startY > 0 then
           _MAKEDOC("Lunar phases from " & startY & " to " & endY & linefeed & linefeed & linefeed, 500, 800)
       else
           _MAKEDOC("Lunar phases for " & startY & linefeed & linefeed & linefeed, 500, 800)
       end if
       
       repeat with dyr from startY to endY
           
           repeat with mnth from 1 to 12
               _PRINT(item mnth of monthNames & " " & dyr & linefeed)
               set MoonPhases to moonphasesForMonth(dyr, mnth)
               repeat with aPhase in MoonPhases
                   local aJde, yr_, mnt_, dy_, hr_, min_, sec_
                   
                   set {yr_, mnt_, dy_} to gregorianCalendarTriplet from item 1 of aPhase
                   set dateStr to formatDate(yr_, mnt_, dy_)
                   set {hr_, min_, sec_} to timeTriplet from item 1 of aPhase
                   set timeStr to formatClockTime(hr_, min_, sec_)
                   _PRINT(dateStr & " : " & timeStr & tab & item (item 2 of aPhase) of exactStates & return)
               end repeat
               _PRINT(linefeed)
           end repeat
           _PRINT(linefeed)
       end repeat
   end createMoonPhasesTable
   
   
   on run
       global am_pm_clock, us_citizen -- "System variables" for configuring the locale. (Date/clock).
       -- resets whatever values, so that any script copying will work with the receivers locale.
       set am_pm_clock to undef()
       set us_citizen to undef()
       
       
       moonPhaseOfToday()
       createMoonPhasesTable()
   end run
   
   on undef()
       return
   end undef
   
   
   on moonPhaseOfToday()
       script o
           property L : {}
       end script
       local YR, mnth, dy, curJDe
       tell (current date)
           set curJDe to (JulianDayNumber of me for it) + 0.5 div 1
           set {YR, mnth, dy} to {its year, its month as integer, its day}
       end tell
       
       set o's L to moonphasesForMonth(YR, mnth)
       local phase
       set phase to ""
       repeat with i from 1 to (count o's L)
           if curJDe < (item 1 of item i of o's L) div 1 then
               set phase to item (item 2 of item i of o's L) of movingStates
               exit repeat
           else if curJDe = (item 1 of item i of o's L) div 1 then
               set phase to item (item 2 of item i of o's L) of exactStates
               exit repeat
           end if
       end repeat
       if phase = "" then
           set phase to item (((item 2 of item (count o's L) of o's L) + 1) mod 4) of movingStates
       end if
       local msg
       #        set msg to formatDate(YR, mnth, dy) & return & return & "The current Phase is " & phase & "."
       set msg to formatDate(YR, mnth, dy) & return & return & phase & "."
       
       tell application (path to frontmost application as text)
           display dialog msg with title "Todays Moon Phase:" buttons {"Quit", "Make Table"} cancel button 1 default button 2 with icon note
       end tell
   end moonPhaseOfToday
   
   
   to sort_items from alist
       -- Kai
       script o
           property L : alist
       end script
       repeat with i from 2 to count o's L
           set V to o's L's item i
           repeat with i from (i - 1) to 1 by -1
               tell o's L's item i to if item 1 of V < item 1 of it then
                   set o's L's item (i + 1) to it
               else
                   set o's L's item (i + 1) to V
                   exit repeat
               end if
           end repeat
           if i is 1 and item 1 of V < item 1 of o's L's beginning then set o's L's item 1 to V
       end repeat
   end sort_items
   
   
   on moonphasesForMonth(YR, mnth)
       script o
           property n : {}
       end script
       
       local tyr, tmnth
       # checks for previous last quarter        
       if mnth = 1 then
           set {tyr, tmnth} to {YR - 1, 11.5}
       else
           set {tyr, tmnth} to {YR, mnth - 1.5}
       end if
       
       local test, JD_adjusted, dyr, mt, dy, L
       
       set {test, L} to {newMoon(tyr, tmnth), {}}
       set JD_adjusted to adjustJDNforLocalTZ(test)
       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
       if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, aNewMoon}, JD_adjusted}
       
       set test to firstQuarter(tyr, tmnth)
       set JD_adjusted to adjustJDNforLocalTZ(test)
       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
       if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, tFirstQuarter}, JD_adjusted}
       
       
       set test to fullMoon(tyr, tmnth)
       set JD_adjusted to adjustJDNforLocalTZ(test)
       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
       if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, aFullMoon}, JD_adjusted}
       
       set test to lastQuarter(tyr, tmnth)
       set JD_adjusted to adjustJDNforLocalTZ(test)
       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
       if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, tThirdQuarter}, JD_adjusted}
       
       # and if we have an extra new moon!
       set test to newMoon(YR, mnth)
       set JD_adjusted to adjustJDNforLocalTZ(test)
       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
       if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aNewMoon}, JD_adjusted}
       
       if length of L < 5 then
           # There can't be more than five moonphases of a month, so the very moment we have five
           # - We bail out!
           set test to firstQuarter(YR, mnth)
           set JD_adjusted to adjustJDNforLocalTZ(test)
           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
           if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tFirstQuarter}, JD_adjusted}
           
           if length of L < 5 then
               set test to fullMoon(YR, mnth)
               set JD_adjusted to adjustJDNforLocalTZ(test)
               set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
               if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aFullMoon}, JD_adjusted}
               
               if length of L < 5 then
                   set test to lastQuarter(YR, mnth)
                   set JD_adjusted to adjustJDNforLocalTZ(test)
                   set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                   if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tThirdQuarter}, JD_adjusted}
                   
                   if length of L < 5 then
                       
                       # Before someone starts to think that I should have coalesced
                       # the four nearly identical handlers and iterated over the handler here
                       # instead of having it spread out boringly one by oneinto one, I'd like
                       # to retort that making code shorter, by adding if tests, and loops only
                       # makes it easier to maintain, not faster. It will also add to overall
                       # complexity here.
                       
                       # the phases that should be there, but not necessarily are!
                       set JD_adjusted to findNewMoonForMonth(YR, mnth)
                       if JD_adjusted ≠ 0 then
                           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                           if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aNewMoon}, JD_adjusted}
                       end if
                       if length of L < 5 then
                           set JD_adjusted to findFirstQuarterForMonth(YR, mnth)
                           if JD_adjusted ≠ 0 then
                               set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                               if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tFirstQuarter}, JD_adjusted}
                           end if
                           if length of L < 5 then
                               set JD_adjusted to findFullMoonForMonth(YR, mnth)
                               if JD_adjusted ≠ 0 then
                                   set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                                   if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aFullMoon}, JD_adjusted}
                               end if
                               if length of L < 5 then
                                   set JD_adjusted to findlastQuarterForMonth(YR, mnth)
                                   if JD_adjusted ≠ 0 then
                                       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                                       if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tThirdQuarter}, JD_adjusted}
                                   end if
                                   
                                   # additional phases that eventually may be there
                                   if length of L < 5 then
                                       if mnth = 12 then
                                           set {tyr, tmnth} to {YR + 1, 11}
                                       else
                                           set {tyr, tmnth} to {YR, mnth + 1}
                                       end if
                                       
                                       set test to newMoon(tyr, tmnth)
                                       set JD_adjusted to adjustJDNforLocalTZ(test)
                                       set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                                       if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, aNewMoon}
                                       
                                       if length of L < 5 then
                                           set test to firstQuarter(tyr, tmnth)
                                           set JD_adjusted to adjustJDNforLocalTZ(test)
                                           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                                           if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, tFirstQuarter}
                                           
                                           if length of L < 5 then
                                               set test to fullMoon(tyr, tmnth)
                                               set JD_adjusted to adjustJDNforLocalTZ(test)
                                               set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                                               if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, aFullMoon}
                                               
                                               if length of L < 5 then
                                                   set test to lastQuarter(tyr, tmnth)
                                                   set JD_adjusted to adjustJDNforLocalTZ(test)
                                                   set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
                                                   if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, tThirdQuarter}
                                               end if
                                           end if
                                       end if
                                   end if
                               end if
                           end if
                       end if
                   end if
               end if
           end if
       end if
       
       sort_items from L
       -- Kai Edwards!
       return L
   end moonphasesForMonth
   
   
   
   on findNewMoonForMonth(YR, mnth)
       local JD, origMonth, dyr, mt, dy, forwardDir
       set {origMonth, forwardDir} to {mnth, false}
       set mnth to mnth - 0.75 # probability!!
       repeat while true
           set JD to newMoon(YR, mnth)
           set JD_adjusted to adjustJDNforLocalTZ(JD)
           
           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
           if mt = origMonth then
               exit repeat
           else if (mnth div 1 + 1) = 1 and mt = 12 then
               # the month we got returned was lesser than the
               # the one we looked for.
               set mnth to 1.25
           else if (mnth div 1 + 1) = 12 and mt = 1 then
               # the month we got returned was greater then the
               # one we looked for.
               set mnth to 10.75
           else if mt = 12 and origMonth = 1 then
               # the month we got returned was greater than the
               # one we looked for.
               set forwardDir to true
               set mnth to mnth + 0.125
           else if mt = 1 and origMonth = 12 then
               # the month we got returned was greater
               set mnth to mnth - 0.25
           else if origMonth < 12 and mt < origMonth then
               # the month we got returned was lesser than
               # we looked for
               set forwardDir to true
               set mnth to mnth + 0.125
           else
               # if we have been forward, and then thrown backwards
               # then the cycle simply doesn't exist for this month!
               if forwardDir then return 0
               # the month was greater than the one we looked for.
               set mnth to mnth - 0.25
           end if
       end repeat
       return JD_adjusted
   end findNewMoonForMonth
   
   
   on newMoon(YR, mnt)
       local k, t, JD, M_, MPrime, F_, Omega, E_
       set k to approxK(YR + (mnt / 12)) div 1 # new moon
       set t to centurieJuliae(k)
       set JD to meanPhaseOfMoon(t, k)
       set JD to computeAdditionalCorrectionsAllPhases(JD, t, k)
       set M_ to sunsMeanAnamoly(t, k)
       set MPrime to moonsMeanAnamoly(t, k)
       set F_ to moonsArgOfLatitude(t, k)
       set Omega to longOfAscendingNodeOfLunarOrbit(t, k)
       set E_ to eccentricityOfEarthsOrbit(t)
       
       return JD + (-0.4072 * (sin MPrime) + ¬
           0.17241 * E_ * (sin M_) + ¬
           0.01608 * (sin (2 * MPrime)) + ¬
           0.01039 * (sin (2 * F_)) + ¬
           0.00739 * E_ * (sin (MPrime - M_)) + ¬
           -0.00514 * E_ * (sin (MPrime + M_)) + ¬
           0.00208 * E_ * E_ * (sin (2 * M_)) + ¬
           -0.00111 * (sin (MPrime - 2 * F_)) + ¬
           -5.7E-4 * (sin (MPrime + 2 * F_)) + ¬
           5.6E-4 * E_ * (sin (2 * MPrime + M_)) + ¬
           -4.2E-4 * (sin (3 * MPrime)) + ¬
           4.2E-4 * E_ * (sin (M_ + 2 * F_)) + ¬
           3.8E-4 * E_ * (sin (M_ - 2 * F_)) + ¬
           -2.4E-4 * E_ * (sin (2 * MPrime - M_)) + ¬
           -1.7E-4 * (sin Omega) + ¬
           -7.0E-5 * (sin (MPrime + 2 * M_)) + ¬
           4.0E-5 * (sin (2 * MPrime - 2 * F_)) + ¬
           4.0E-5 * (sin (3 * M_)) + ¬
           3.0E-5 * (sin (MPrime + M_ - 2 * F_)) + ¬
           3.0E-5 * (sin (2 * MPrime + 2 * F_)) + ¬
           -3.0E-5 * (sin (MPrime + M_ + 2 * F_)) + ¬
           3.0E-5 * (sin (MPrime - M_ + 2 * F_)) + ¬
           -2.0E-5 * (sin (MPrime - M_ - 2 * F_)) + ¬
           -2.0E-5 * (sin (3 * MPrime + M_)) + ¬
           2.0E-5 * (sin (4 * MPrime)))
       
   end newMoon
   
   # Before someone starts to think that I should have coalesced
   # the four nearly identical handlers into one I'd like to retort
   # That making code shorter, by adding if tests, only makes it
   # easier to maintain, not faster.
   
   on findFirstQuarterForMonth(YR, mnth)
       local JD, origMonth, dyr, mt, dy, forwardDir
       set {origMonth, forwardDir} to {mnth, false}
       set mnth to mnth - 0.75
       repeat while true
           set JD to firstQuarter(YR, mnth)
           set JD_adjusted to adjustJDNforLocalTZ(JD)
           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
           if mt = origMonth then
               exit repeat
           else if (mnth div 1 + 1) = 1 and mt = 12 then
               # the month we got returned was lesser than the
               # the one we looked for.
               set mnth to 1.25
           else if (mnth div 1 + 1) = 12 and mt = 1 then
               # the month we got returned was greater then the
               # one we looked for.
               set mnth to 10.75
           else if mt = 12 and origMonth = 1 then
               # the month we got returned was greater than the
               # one we looked for.
               set forwardDir to true
               set mnth to mnth + 0.125
           else if mt = 1 and origMonth = 12 then
               # the month we got returned was greater
               set mnth to mnth - 0.25
           else if origMonth < 12 and mt < origMonth then
               # the month we got returned was lesser than
               # we looked for
               set forwardDir to true
               set mnth to mnth + 0.125
           else
               # if we have been forward, and then thrown backwards
               # then the cycle simply doesn't exist for this month!
               if forwardDir then return 0
               # the month was greater than the one we looked for.
               set mnth to mnth - 0.25
           end if
       end repeat
       return JD_adjusted
   end findFirstQuarterForMonth
   
   on firstQuarter(YR, mnt)
       local JD, W_
       set {JD, W_} to commonQuarters(YR, mnt, 0.25)
       return JD + W_
   end firstQuarter
   
   on commonQuarters(YR, mnt, factor)
       local k, t, JD, M_, MPrime, F_, Omega, E_, W_
       set k to approxK(YR + (mnt / 12)) div 1 + factor # Last Quater
       set t to centurieJuliae(k)
       set JD to meanPhaseOfMoon(t, k)
       set JD to computeAdditionalCorrectionsAllPhases(JD, t, k)
       set M_ to sunsMeanAnamoly(t, k)
       set MPrime to moonsMeanAnamoly(t, k)
       set F_ to moonsArgOfLatitude(t, k)
       set Omega to longOfAscendingNodeOfLunarOrbit(t, k)
       set E_ to eccentricityOfEarthsOrbit(t)
       set W_ to 0.00306 - 3.8E-4 * E_ * (cos M_) + 2.6E-4 * (cos MPrime) - 2.0E-5 * (cos (MPrime - M_)) + 2.0E-5 * (cos (MPrime + M_)) + 2.0E-5 * (cos (2 * F_))
       
       return {JD - 0.62801 * (sin MPrime) + ¬
           0.17172 * E_ * (sin M_) + ¬
           -0.01183 * E_ * (sin (MPrime + M_)) + ¬
           0.00862 * (sin (2 * MPrime)) + ¬
           0.00804 * (sin (2 * F_)) + ¬
           0.00454 * E_ * (sin (MPrime - M_)) + ¬
           0.00204 * E_ * E_ * (sin (2 * M_)) + ¬
           -0.0018 * (sin (MPrime - (2 * F_))) + ¬
           -7.0E-4 * (sin (MPrime + (2 * F_))) + ¬
           -4.0E-4 * (sin (3 * MPrime)) + ¬
           -3.4E-4 * E_ * (sin (2 * MPrime - M_)) + ¬
           3.2E-4 * E_ * (sin (M_ + (2 * F_))) + ¬
           3.2E-4 * E_ * (sin (M_ - (2 * F_))) + ¬
           -2.8E-4 * E_ * E_ * (sin (MPrime + (2 * M_))) + ¬
           2.7E-4 * E_ * (sin ((2 * MPrime) + M_)) + ¬
           -1.7E-4 * (sin Omega) + ¬
           -5.0E-5 * (sin (MPrime - M_ - (2 * F_))) + ¬
           4.0E-5 * (sin ((2 * MPrime) + (2 * F_))) + ¬
           -4.0E-5 * (sin (MPrime + M_ + (2 * F_))) + ¬
           4.0E-5 * (sin (MPrime - (2 * M_))) + ¬
           3.0E-5 * (sin (MPrime + M_ - (2 * F_))) + ¬
           3.0E-5 * (sin (3 * M_)) + ¬
           2.0E-5 * (sin ((2 * MPrime) - (2 * F_))) + ¬
           2.0E-5 * (sin (MPrime - M_ + (2 * F_))) + ¬
           -2.0E-5 * (sin ((3 * MPrime) + M_)), W_}
   end commonQuarters
   
   on findFullMoonForMonth(YR, mnth)
       local JD, origMonth, dyr, mt, dy, forwardDir
       set {origMonth, forwardDir} to {mnth, false}
       set mnth to mnth - 0.75
       repeat while true
           set JD to fullMoon(YR, mnth)
           set JD_adjusted to adjustJDNforLocalTZ(JD)
           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
           if mt = origMonth then
               exit repeat
           else if (mnth div 1 + 1) = 1 and mt = 12 then
               # the month we got returned was lesser than the
               # the one we looked for.
               set mnth to 1.25
           else if (mnth div 1 + 1) = 12 and mt = 1 then
               # the month we got returned was greater then the
               # one we looked for.
               set mnth to 10.75
           else if mt = 12 and origMonth = 1 then
               # the month we got returned was greater than the
               # one we looked for.
               set forwardDir to true
               set mnth to mnth + 0.125
           else if mt = 1 and origMonth = 12 then
               # the month we got returned was greater
               set mnth to mnth - 0.25
           else if origMonth < 12 and mt < origMonth then
               # the month we got returned was lesser than
               # we looked for
               set forwardDir to true
               set mnth to mnth + 0.125
           else
               # if we have been forward, and then thrown backwards
               # then the cycle simply doesn't exist for this month!
               if forwardDir then return 0
               # the month was greater than the one we looked for.
               set mnth to mnth - 0.25
           end if
       end repeat
       return JD_adjusted
   end findFullMoonForMonth
   
   on fullMoon(YR, mnt)
       local k, t, JD, M_, MPrime, F_, Omega, E_
       set k to approxK(YR + (mnt / 12)) div 1 + 0.5 # full Moon
       set t to centurieJuliae(k)
       set JD to meanPhaseOfMoon(t, k)
       set JD to computeAdditionalCorrectionsAllPhases(JD, t, k)
       set M_ to sunsMeanAnamoly(t, k)
       set MPrime to moonsMeanAnamoly(t, k)
       set F_ to moonsArgOfLatitude(t, k)
       set Omega to longOfAscendingNodeOfLunarOrbit(t, k)
       set E_ to eccentricityOfEarthsOrbit(t)
       
       return JD - 0.40614 * (sin MPrime) + ¬
           0.17302 * E_ * (sin M_) + ¬
           0.01614 * (sin (2 * MPrime)) + ¬
           0.01043 * (sin (2 * F_)) + ¬
           0.00734 * E_ * (sin (MPrime - M_)) + ¬
           -0.00514 * E_ * (sin (MPrime + M_)) + ¬
           0.00209 * E_ * E_ * (sin (2 * M_)) + ¬
           -0.00111 * (sin (MPrime - 2 * F_)) + ¬
           -5.7E-4 * (sin (MPrime + 2 * F_)) + ¬
           5.6E-4 * E_ * (sin (2 * MPrime + M_)) + ¬
           -4.2E-4 * (sin (3 * MPrime)) + ¬
           4.2E-4 * E_ * (sin (M_ + 2 * F_)) + ¬
           3.8E-4 * E_ * (sin (M_ - 2 * F_)) + ¬
           -2.4E-4 * E_ * (sin (2 * MPrime - M_)) + ¬
           -1.7E-4 * (sin Omega) + ¬
           -7.0E-5 * (sin (MPrime + 2 * M_)) + ¬
           4.0E-5 * (sin (2 * MPrime - 2 * F_)) + ¬
           4.0E-5 * (sin (3 * M_)) + ¬
           3.0E-5 * (sin (MPrime + M_ - 2 * F_)) + ¬
           3.0E-5 * (sin (2 * MPrime + 2 * F_)) + ¬
           -3.0E-5 * (sin (MPrime + M_ + 2 * F_)) + ¬
           3.0E-5 * (sin (MPrime - M_ + 2 * F_)) + ¬
           -2.0E-5 * (sin (MPrime - M_ - 2 * F_)) + ¬
           -2.0E-5 * (sin (3 * MPrime + M_)) + ¬
           2.0E-5 * (sin (4 * MPrime))
   end fullMoon
   
   on findlastQuarterForMonth(YR, mnth)
       local JD, origMonth, dyr, mt, dy, forwardDir
       set {origMonth, forwardDir} to {mnth, false}
       set mnth to mnth - 0.75
       repeat while true
           set JD to lastQuarter(YR, mnth)
           set JD_adjusted to adjustJDNforLocalTZ(JD)
           set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
           if mt = origMonth then
               exit repeat
           else if (mnth div 1 + 1) = 1 and mt = 12 then
               # the month we got returned was lesser than the
               # the one we looked for.
               set mnth to 1.25
           else if (mnth div 1 + 1) = 12 and mt = 1 then
               # the month we got returned was greater then the
               # one we looked for.
               set mnth to 10.75
           else if mt = 12 and origMonth = 1 then
               # the month we got returned was greater than the
               # one we looked for.
               set forwardDir to true
               set mnth to mnth + 0.125
           else if mt = 1 and origMonth = 12 then
               # the month we got returned was greater
               set mnth to mnth - 0.25
           else if origMonth < 12 and mt < origMonth then
               # the month we got returned was lesser than
               # we looked for
               set forwardDir to true
               set mnth to mnth + 0.125
           else
               # if we have been forward, and then thrown backwards
               if forwardDir then return 0
               set mnth to mnth - 0.25
           end if
           
       end repeat
       return JD_adjusted
   end findlastQuarterForMonth
   
   on lastQuarter(YR, mnt)
       local JD, W_
       set {JD, W_} to commonQuarters(YR, mnt, 0.75)
       return JD - W_
   end lastQuarter
   
   on computeAdditionalCorrectionsAllPhases(JDE, t, k)
       return (JDE + 3.25E-4 * (sin (map360(299.77 + 0.107408 * k - 0.009173 * t * t) * degrad)) + ¬
           1.65E-4 * (sin (map360(251.88 + 0.016321 * k) * degrad)) + ¬
           1.64E-4 * (sin (map360(251.83 + 26.651886 * k) * degrad)) + ¬
           1.26E-4 * (sin (map360(349.42 + 36.412478 * k) * degrad)) + ¬
           1.1E-4 * (sin (map360(84.66 + 18.206239 * k) * degrad)) + ¬
           6.2E-5 * (sin (map360(141.74 + 53.303771 * k) * degrad)) + ¬
           6.0E-5 * (sin (map360(207.14 + 2.453732 * k) * degrad)) + ¬
           5.6E-5 * (sin (map360(154.84 + 7.30686 * k) * degrad)) + ¬
           4.7E-5 * (sin (map360(34.52 + 27.261239 * k) * degrad)) + ¬
           4.2E-5 * (sin (map360(207.19 + 0.121824 * k) * degrad)) + ¬
           4.0E-5 * (sin (map360(291.34 + 1.844379 * k) * degrad)) + ¬
           3.7E-5 * (sin (map360(161.72 + 24.198154 * k) * degrad)) + ¬
           3.5E-5 * (sin (map360(239.56 + 25.513099 * k) * degrad)) + ¬
           2.3E-5 * (sin (map360(331.55 + 3.592518 * k) * degrad)))
   end computeAdditionalCorrectionsAllPhases
   
   
   on eccentricityOfEarthsOrbit(t)
       # (Meeus 45.6)
       return (1 + (t * (-0.002516 + (t * -7.4E-6))))
   end eccentricityOfEarthsOrbit
   
   on longOfAscendingNodeOfLunarOrbit(k, t)
       # Meeus ( 47.7)
       return (map360(124.7746 - (k * 1.5637558) + (t * t * (0.0020691 + (t * 2.15E-6)))) * degrad)
   end longOfAscendingNodeOfLunarOrbit
   
   on moonsArgOfLatitude(t, k)
       # meeus ( 47.6 )
       return (map360((160.7108 + (k * 390.67050274)) + (t * t * (-0.0016341 + (t * (-2.27E-6 + (t * 1.1E-8)))))) * degrad)
   end moonsArgOfLatitude
   
   on moonsMeanAnamoly(t, k)
       # Meeus (47.5)
       return (map360(201.5643 + (k * 385.81693528) + (t * t * (0.0107438 + (t * 1.239E-5 - (t * 5.8E-8))))) * degrad)
   end moonsMeanAnamoly
   
   on sunsMeanAnamoly(t, k)
       # Meeus (47.4)
       return (map360(2.5534 + (k * 29.10535669) + (t * t * (-2.18E-5 + (t * -1.1E-7)))) * degrad)
   end sunsMeanAnamoly
   
   on meanPhaseOfMoon(t, k)
       # Meeus (47.1)
       return 2.45155009765E+6 + (29.530588853 * k) + (t * (t * (1.337E-4 + (t * (-1.5E-7 + (t * 7.3E-10))))))
   end meanPhaseOfMoon
   
   on centurieJuliae(k)
       # Meeus (47.3)
       return (k / 1236.85)
   end centurieJuliae
   
   on approxK(decyear)
       # Meeus (47.2)
       set a to ((decyear - 2000) * 12.3685)
       return round ((decyear - 2000) * 12.3685) rounding to nearest
   end approxK
   
   # ----
   
   on map360(aDecNumber)
       set aDecNumber to aDecNumber mod 360
       if aDecNumber < 0 then set aDecNumber to aDecNumber + 360
       return aDecNumber
   end map360
   
   # end script
   
   on GregorianDate for theJDN
       -- Added 2013.10.7
       -- Nigel Garvey [url]http://macscripter.net/viewtopic.php?pid=167391#p167391[/url] post # 6
       return JD0 + theJDN * days
   end GregorianDate
   
   on TZtoTZ(TZ1date, TZ1, TZ2)
       # Nigel Garvey
       return (do shell script ("eraTime=$(TZ=" & TZ1 & " date -jf '%FT%T' '" & (TZ1date as «class isot» as string) & "' '+%s') ; TZ=" & TZ2 & " date -r \"$eraTime\" '+%FT%T'") as «class isot») as date
   end TZtoTZ
   
   
   on JulianDayNumber for theDate
       -- Added 2013.10.7
       -- Nigel Garvey [url]http://macscripter.net/viewtopic.php?pid=167391#p167391[/url] post # 6
       return (theDate - JD0) div days
   end JulianDayNumber
   
   
   on timeTriplet from JDN
       # This returns as standard time triplet, where the first point of time of a day is 00:00
       # and the last point of a time is 23:59:59.9999999 and so on..
       local timeFraction, _hours, _minutes, _seconds
       set timeFraction to (JDN mod 1) * 24
       set _hours to (timeFraction div 1 as integer)
       set timeFraction to (timeFraction - _hours)
       set _minutes to timeFraction * 60 div 1
       set _seconds to (timeFraction - (_minutes / 60)) * 3600
       return {_hours, _minutes, _seconds}
   end timeTriplet
   
   on gregorianCalendarTriplet from JDN
       # This is an algorithm by Richards to convert a Julian Day Number, J, to a date in the Gregorian calendar
       # (proleptic, when applicable). Richards does not state which dates the algorithm is valid for.
       # it seems to work for dates after, BUT NOT INCLUDING 15th of oct 1582.
       # I have extended it to return a hour, min, sec triplet for the fractional part of the JDN.
       local y, m, n, R, p, V, s, w, b, C, F, g, h, D, mt, YR, JD
       set y to 4716
       set j to 1401
       set m to 2
       set n to 12
       set R to 4
       set p to 1461
       set V to 3
       set u to 5
       set s to 153
       set w to 2
       set b to 274277
       set C to -38
       
       set JD to JDN div 1
       set F to JD + j + (((4 * JD + b) div 146097) * 3) div 4 + C
       set e to R * F + V
       set g to (e mod p) div R
       set h to u * g + w
       set D to (h mod s) div u + 1
       set mt to (h div s + m) mod n + 1
       set YR to e div p - y + (n + m - mt) div n
       return {YR, mt, D}
   end gregorianCalendarTriplet
   
   on civilTimeToGmtFromLocalTZInSeconds(anASDate)
       # Nigel Garvey and Adam Bell made this possible!
       return (anASDate - TZtoTZ(anASDate, (do shell script "readlink /etc/localtime |sed -n 's_/[^/]*/[^/]*/[^/]*/\\([^/]*\\)/\\([^/]*\\).*_\\1/\\2_p'"), "GMT"))
   end civilTimeToGmtFromLocalTZInSeconds
   
   on JDN(_year, _month, _day, _hour, _min, _sec)
       # From Wikipedia
       # The Julian Date (JD) of any instant is the Julian day number
       # for the preceding noon plus the fraction of the day since
       # that instant. Julian Dates are expressed as a Julian day
       # number with a decimal fraction added. For example, The Julian
       # Date for 00:30:00.0 1 January 2013 is 2456293.520833.
       local a, y, m, _JDN, time_fraction
       set a to ((14 - _month) div 12)
       set y to _year + 4800 - a
       set m to _month + (12 * a) - 3
       # Gregorian Calendar starts October 15, 1582
       if _year > 1582 then
           set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
       else if _year = 1582 and _month > 10 then
           set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
       else if _year = 1582 and _month = 10 and _day ≥ 15 then
           set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
       else
           # Julian calendar before oct 15, 1582
           set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - 32083
       end if
       set _hour to (_hour - 12) / 24
       if _min ≠ 0 then
           set _min to _min / 1440 # minutes per day
       end if
       
       if _sec ≠ 0 then
           set _sec to _sec / 86400 # seconds per day
       end if
       
       set _JDN to _JDN + _hour + _min + _sec
       
       return _JDN
   end JDN
   
   on adjustJDNforLocalTZ(JDE)
       # it adds the time difference to GMT, and adds the 12 hour, to askew the day to the civil
       # time system, this may also increment the date, which is perfectly fine, since a Julian day
       # starts 12 hours before a civil day.
       return (JDE + ((civilTimeToGmtFromLocalTZInSeconds(GregorianDate for (JDE + 0.5))) / 86400) + 0.5)
   end adjustJDNforLocalTZ
   
   on _PRINT(txt)
       # Nigel Garvey made the original version.
       local a
       tell application "TextEdit"
           if not my Debug then
               activate
           else
               do shell script "open -b \"com.apple.textedit\""
           end if
           if ((count its documents) is 0) then make new document
           set a to (get path of its front document)
           try
               if a is "" then set a to missing value
           on error
               set a to missing value
           end try
           if a is not missing value then make new document
           make new paragraph at end of text of front document with data (txt & linefeed) with properties {font:"Andale Mono", size:12}
       end tell
       return true
   end _PRINT
   
   on _MAKEDOC(heading, _width, _height)
       tell application "TextEdit"
           tell (make new document at front)
               set text of it to heading
               set font of every character to "Andale Mono"
               set size of every character to 13
           end tell
           tell its front window
               #    set {item 3 of its bounds, item 4 of its bounds} to {800, 300}
               set _bounds to bounds
               set item 3 of _bounds to (item 1 of _bounds) + _width
               set item 4 of _bounds to (item 2 of _bounds) + _height
               set bounds to _bounds
               
           end tell
       end tell
   end _MAKEDOC
   
   on formatDate(YR, mnt, _day)
       global us_citizen
       local a
       try
           set a to us_citizen
       on error
           set us_citizen to getUsCitizenShip()
       end try
       
       if us_citizen then
           return (pad(mnt) & "/" & pad(_day) & "/" & YR)
       else
           return (pad(_day) & "/" & pad(mnt) & "/" & YR)
       end if
   end formatDate
   
   on setClockFormat(m)
       local clock_format
       considering case
           tell (current date)'s time string
               if it contains "AM" or it contains "PM" then
                   set clock_format to m's big_AMPM
               else if it contains "Am" or it contains "Pm" then
                   set clock_format to m's ususal_AmPm
               else if it contains "am" or it contains "pm" then
                   set clock_format to m's small_ampm
               else
                   set clock_format to m's hour24
               end if
           end tell
       end considering
       return clock_format
   end setClockFormat
   
   on formatClockTime(_hr, _min, _sec)
       script o
           property big_AMPM : 3
           property ususal_AmPm : 2
           property small_ampm : 1
           property hour24 : 0
       end script
       global am_pm_clock
       local a
       try
           set a to am_pm_clock
       on error
           set am_pm_clock to setClockFormat(o)
       end try
       
       if am_pm_clock = o's hour24 then
           return (pad((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)))
       else if am_pm_clock = o's big_AMPM then
           if (_hr div 12) > 0 then
               return (padHourWithBlank((_hr div 12)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " PM")
           else
               return (padHourWithBlank((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " AM")
           end if
       else if am_pm_clock = o's ususal_AmPm then
           if (_hr div 12) > 0 then
               return (padHourWithBlank((_hr div 12)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " Pm")
           else
               return (padHourWithBlank((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " Am")
           end if
       else # am_pm_clock = small_ampm
           if (_hr div 12) > 0 then
               return (padHourWithBlank((_hr div 12)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " pm")
           else
               return (padHourWithBlank((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " am")
           end if
       end if
   end formatClockTime
   
   on padHourWithBlank(n)
       -- "Special Handler" for padding the 12 hour clock with blanks.
       return (text -2 thru -1 of (" " & n as text))
   end padHourWithBlank
   
   to getUsCitizenShip()
       try
           set localeString to (do shell script "defaults read .GlobalPreferences AppleLocale 2>/dev/null || echo en_US |tr '
' ' '"
)
       on error e number n
           return true
       end try
       if localeString is "en_US" then return true
       return false
   end getUsCitizenShip
   
   on pad(n)
       -- 22. September 2013, improved pad, for the cases where n < 0 and n > 100.
       -- Thanks to Yvan Koenig.
       -- Fixed a bug
       if n < -100 then
           return n as text
       else if n < 0 then
           return ("-" & text -2 thru -1 of ("0" & ((n as number) * -1) as text))
           --    error "The pad handler assumes positive values" number 5000
       else if n < 100 then
           return text 2 thru -1 of ((100 + n) as text)
       else
           return n as text
       end if
   end pad
   
   -- Ask the user for the range of dates to be covered.
   -- from a start date, up to, and including enddate.
   on getRangeOfYears(initialYear, initialRange, scriptTitle)
       -- Added 2013.10.7, used in Equation Of Time, high accuracy..
       -- Nigel Garvey [url]http://macscripter.net/viewtopic.php?id=41273[/url]
       -- customized a little by McUsr (All faults are mine!),
       -- so it is able to deliver just one date denoting a range of 1 day.
       -- There are preset ranges of Month, Week, and Year, based on the initalDate.
       local anicon, theMessage, StartYear, EndYear
       set anicon to note
       set theMessage to "Enter a range of years or a single year."
       set StartYear to initialYear
       repeat
           if initialRange = "Century" then
               set EndYear to StartYear + 100
           else if initialRange = "SemiCentury" then
               set EndYear to StartYear + 50
           else if initialRange = "Decade" then
               set EndYear to StartYear + 10
           else if initialRange = "Year" then
               set EndYear to StartYear + 5
               # We'll add one day for a leap year
           end if
           local d1, d2, yearRange
           set d1 to StartYear as text
           set d2 to EndYear as text
           
           
           tell application (path to frontmost application as text)
               set yearRange to text returned of (display dialog theMessage default answer d1 & " - " & d2 with title scriptTitle buttons {"Quit", "Ok"} cancel button 1 default button 2 with icon anicon)
               
           end tell
           local tids, gotOne, yearRangeStart, yearRangeEnd
           try
               set {tids, AppleScript's text item delimiters, gotOne} to {AppleScript's text item delimiters, "-", false}
               set {yearRangeStart, gotOne} to {text item 1 of yearRange as number, true}
               
               set yearRangeEnd to text item 2 of yearRange as number
               set AppleScript's text item delimiters to tids
               exit repeat
           on error
               if gotOne and length of text items of yearRange = 1 then
                   set yearRangeEnd to text item 1 of yearRange as number
                   set AppleScript's text item delimiters to tids
                   exit repeat
               else if gotOne then
                   set today to yearRangeStart
               else
                   set today to initialYear
               end if
               set AppleScript's text item delimiters to tids
               if theMessage does not contain "wasn't" then
                   set theMessage to "The input wasn't valid!" & return & return & theMessage
                   set anicon to caution
               end if
           end try
       end repeat
       return {yearRangeStart, yearRangeEnd}
   end getRangeOfYears
end script

tell MoonPhases to run

Last edited by McUsrII (2013-10-13 03:48:51 pm)


Filed under: Astronomy, Moon

Offline

 

#2 2013-10-11 07:16:36 am

McUsrII
Member
Registered: 2012-11-21
Posts: 3046
Website

Re: Accurate Table of Moonphases

Hello.

I added a property for setting a 12 hour am/pm clock, for those of you that are accustomed to that particular time format. I currently don't know how to easily deduct that info, so you'll have to set the property to true in order to turn the feature on!

I also neglected to say in the original post that the times of the lunar phases are returned according to your Mac's date and time settings, in the preferences.

Offline

 

#3 2013-10-11 09:58:59 am

McUsrII
Member
Registered: 2012-11-21
Posts: 3046
Website

Re: Accurate Table of Moonphases

Hello.

I added functionality to automatically detect clock format: 24 hour, AM/PM, Am/Pm, and am/pm (for prosperity).

Offline

 

#4 2013-10-12 06:03:08 am

McUsrII
Member
Registered: 2012-11-21
Posts: 3046
Website

Re: Accurate Table of Moonphases

Hello.

I changed the detection of clock format that is am/pm or a 24 hour clock to be a private subsystem, not relying on global properties, as I deemed that the soundest way to implement it. Not that it is much better if the user just shares a copy of his script, as the previous clock format will still be in effect until the script is edited, but the use of this  design principle makes it easier to implement the functionality (less editing), and steal. smile

The script is still in post # 1, and you still need Satimage.osax.

Last edited by McUsrII (2013-10-12 06:05:45 am)

Offline

 

#5 2013-10-13 03:52:24 pm

McUsrII
Member
Registered: 2012-11-21
Posts: 3046
Website

Re: Accurate Table of Moonphases

Hello.

I just added forced initialization of the variables that controls the format of clock, and date according to locale.

I do that with an undef() handler that just has a return statement in it, effectively undefining the variable that is assigned the result.

It is a funny quirk here: It works if you do it in a single statement, but at least Script Debugger barks, when I try to do it simultaneously to several variables at once with a  list statement. smile

Offline

 

#6 2015-05-20 03:02:51 pm

McUsrII
Member
Registered: 2012-11-21
Posts: 3046
Website

Re: Accurate Table of Moonphases

Hello.

Here is another script, that shows moonphases for a month at a time. Which should compute a little bit faster than the script above, you really need to have Satimage.osax installed, or provide your handlers for trigs and so on.

Applescript:

# © 2013 McUsr. You may not post this anywhere else. [url]http://macscripter.net/viewtopic.php?pid=180859#p180859[/url]
# Revised after advice by Nigel Garvey,
# New Moon January the 6 year 2000
# [url]http://kalender-365.de/manekalender.php?yy=2000[/url]
# [url]http://www.calendar-365.com/moon/moon-phases.html[/url]
# [url]http://en.wikipedia.org/wiki/Lunar_phase[/url]
# ----
# [url]http://en.wikipedia.org/wiki/Full_moon_cycle[/url]
# [url]http://en.wikipedia.org/wiki/New_moon[/url]
# [url]http://en.wikipedia.org/wiki/Full_moon[/url]

property fullMoonsPerFullmoonCycle : 14
property daysPerFullmoonCycle : 411.78443
# what kind of moon cycle to follow, a synodic month contains the
# phases cf. link.
property synodicMonth : 29.530588853
property EpochJD2000 : 2.451545E+6

global w, exc
set TGMT to ((time to GMT) / hours) / 24
# log "TGMT : " & TGMT

# ENTER YEAR AND MONTH FOR THE MOON HERE
set sYear to 2015
set sMonth to 5
set m to fullMoonNumberForCalendricMonth(sYear, sMonth, 1)

set JD to JDN(sYear, sMonth, 1, 0, 0, 0)
# log "JD = " & JD

set JDE to calculateJDE(JD, sYear, sMonth)
# log "JDE = " & JDE

set t to (T_inJulianCenturies for JDE)
# log "T = " & T

set MAS to (MeanAnomalyForSunInRad for t)
# log "MAS = " & MAS
set MAM to (MeanAnomalyForMoonInRad for t)
# log "MAM = " & MAM
set exc to (MoonsAmplitude for t)
# log "E = " & exc
# Correction of quarters
set w to 0.0036 - (3.8E-4 * (MoonsAmplitude for t) * (cos MAS)) + (2.6E-4 * (cos MAM)) - (2.0E-5 * (cos (MAM - MAS))) + (2.0E-5 * (cos (MAS + MAM))) + (2.0E-5 * (cos (2 * (MoonsArgOfLatitude for t))))
# log " W = " & w
# NEW MOON
#log "moon number = " & m
repeat while true
   set dnr to meanSyzygyForMoonNr(m, 5.597661) # New
   #    log "mean syzygy = " & dnr
   # Calculations in order to find the true syzygy
   
   set tdnr to (trueSyzygyForNewMoonNr(dnr, MAM, MAS) + (EpochJD2000) - 1)
   set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) + 0.5 - TGMT
   #    log "true syzygy = " & tdnr
   #    log "tdnr - Jd = " & (tdnr - JD)
   if tdnr - JD > 30 then
       set m to m - 1
   else if JD - tdnr > 30 then
       set m to m + 1
   else
       exit repeat
   end if
end repeat
set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to "Calendar date for New Moon = Yr " & yr & " mt = " & MT & " d = " & d & return

# FIRST QUARTER

set dnr to meanSyzygyForMoonNr(m, (12.980308 + (0.5 / exc))) # FirstQuarter
# log "mean syzygy = " & dnr
# Calculations in order to find the true syzygy

set tdnr to (trueSyzygyForFirstQuarterMoonNr(dnr, MAM, MAS) + EpochJD2000 - 1)
set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) # + 0.5 - TGMT
# log "true syzygy = " & tdnr
set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to MTEXT & "Calendar date for First Quarter = Yr " & yr & " mt = " & MT & " d = " & d & return

# FULL MOON

set dnr to meanSyzygyForMoonNr(m, 20.362955) # Full
# log "mean syzygy = " & dnr
# Calculations in order to find the true syzygy

set tdnr to (trueSyzygyForFullMoonNr(dnr, MAM, MAS) + EpochJD2000 - 1)
set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) # + 0.5 - TGMT
# log "true syzygy = " & tdnr


set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to MTEXT & "Calendar date for Full Moon = Yr " & yr & " mt = " & MT & " d = " & d & return


# LAST QUARTER

set dnr to meanSyzygyForMoonNr(m, (27.745602 + (1 / exc))) # LastQuarter
# log "mean syzygy = " & dnr
# Calculations in order to find the true syzygy

set tdnr to (trueSyzygyForLastQuarterMoonNr(dnr, MAM, MAS) + EpochJD2000 - 1)
set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) # + 0.5 - TGMT
# log "true syzygy = " & tdnr

set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to MTEXT & "Calendar date for Last Quarter = Yr " & yr & " mt = " & MT & " d = " & d & return

tell application "TextEdit"
   make new document
   set text of its document 1 to MTEXT
end tell


to MeanAnomalyForSunInRad for t
   # Meeus (45.3)
   return (((357.52772 + (3.599905034E+4 * t) + (1.603E-4 * (t ^ 2)) + ((t ^ 3) / 300000)) mod 360) * pi / 180)
end MeanAnomalyForSunInRad

to MeanAnomalyForMoonInRad for t
   # Meeus ( 45.4)
   return (((134.96298 + (4.77198867398E+5 * t) + (0.0086972 * (t ^ 2)) + ((t ^ 3) / 56250)) mod 360) * pi / 180)
end MeanAnomalyForMoonInRad

to MoonsArgOfLatitude for t
   # Meeus (45.5)
   return (((93.2720993 + (4.832020175273E+5 * t) - (0.0034029 * (t ^ 2)) - ((t ^ 3) / 3526000) + ((t ^ 4) / 8.6331E+8)) mod 360) * pi / 180)
end MoonsArgOfLatitude

to MoonsAmplitude for t
   # Meeus 45.6
   return (1 - (0.002516 * t) - (7.4E-7 * (t ^ 2)))
end MoonsAmplitude
to calculateJDE(JD, yr, mnt)
   # Meus p. 131, but the deltaT is from Nasa
   
   return (JD + ((deltaT for (decimalYear(yr, mnt))) / 86400))
end calculateJDE

# returns the time in Julian Centuries
# JDE is a JD extrepolated to Dynamical time.

to T_inJulianCenturies for aJDE
   return ((aJDE - 2451545) / 36525)
end T_inJulianCenturies

to fullMoonNumberForCalendricMonth(y, m, d)
   # The Calendric Month is a Gregorian Month.
   
   # get todays julian day number smart to start with the first of the
   # month, maybe to add the coeffecient of Terrestial Time too.
   set curJD to JDN(y, m, d, 0, 0, 0)
   set JDSinceEpoch to curJD - EpochJD2000
   #    log " Julian days since epoch: " & JDSinceEpoch
   set fullMoonCycles to JDSinceEpoch div daysPerFullmoonCycle
   #    log "fullMoonCycles since EpochJD2000: " & fullMoonCycles
   #        set daysLeftOver to JDSinceEpoch mod fullMoonCycles
   set daysLeftOver to JDSinceEpoch mod daysPerFullmoonCycle
   #    log "daysLeftOver: " & daysLeftOver
   set extraMoonPhases to daysLeftOver div synodicMonth
   #    log "extraMoonPhases: " & extraMoonPhases
   # calulates the moonNumber we are after:
   set moonNumber to (fullMoonCycles * fullMoonsPerFullmoonCycle) + extraMoonPhases
   # So, now we know the moon number we are looking for the full-moon for.
   # I kind of know the current moons "age" now, if it should matter later
   return moonNumber
   
end fullMoonNumberForCalendricMonth

# Simplified: thanks to Nigel Garvey.
to meanSyzygyForMoonNr(moonNumber, primeFactor)
   local d, corrD
   # [url]http://en.wikipedia.org/wiki/Full_moon[/url]
   # Returns the mean syzygy, which must be corrected.
   set d to (primeFactor + (29.530588861 * moonNumber) + ((102.026 * 1.0E-11) * (moonNumber ^ 2)))
   #    log "d: " & d
   set corrD to d + (-7.39E-4 - (235 * 1.0E-11) * (moonNumber ^ 2))
   #    log "corrD: " & corrD
   return corrD
end meanSyzygyForMoonNr


to trueSyzygyForNewMoonNr(meansyzygy, MAM, MAS)
   # instant dynamical time, (21.1) Meeus
   # (45.6)
   return (meansyzygy + (-0.4072 * (sin MAM)) + (0.01608 * (sin (2 * MAM))) + (0.17241 * (sin MAS)))
end trueSyzygyForNewMoonNr

to trueSyzygyForFirstQuarterMoonNr(meansyzygy, MAM, MAS)
   # instant dynamical time, (21.1) Meeus
   # (45.6)
   return (meansyzygy + (-0.62801 * (sin MAM)) + (0.00862 * (sin (2 * MAM))) + (0.17172 * (sin MAS)) + w)
end trueSyzygyForFirstQuarterMoonNr

to trueSyzygyForFullMoonNr(meansyzygy, MAM, MAS)
   # instant dynamical time, (21.1) Meeus
   # (45.6)
   return (meansyzygy + (-0.40614 * (sin MAM)) + (0.01614 * (sin (2 * MAM))) + (0.17302 * (sin MAS)))
end trueSyzygyForFullMoonNr

to trueSyzygyForLastQuarterMoonNr(meansyzygy, MAM, MAS)
   # instant dynamical time, (21.1) Meeus
   # (45.6)
   return (meansyzygy + (-0.62801 * (sin MAM)) + (0.00862 * (sin (2 * MAM))) + (0.17172 * (sin MAS)) - w)
end trueSyzygyForLastQuarterMoonNr

on timeInJulianCenturies(JDE)
   # returns the time in Julian Centuries
   # JDE is a JD extrepolated to Dynamical time.
end timeInJulianCenturies

# aDecimalYear = year + (month - 0.5) / 12
to decimalYear(yr, mnth)
   return (yr + (mnth - 0.5) / 12)
end decimalYear

to deltaT for aDecimalYear
   # [url]http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html[/url]
   # simplifed, halved the contents of the if tests,
   # and speeded up the multiplication of powers thanks to Nigel Garvey
   local u
   if aDecimalYear < -500 then
       set u to (aDecimalYear - 1820) / 100
       return (-20 + 32 * (u ^ 2))
   else if aDecimalYear < 500 then
       set u to aDecimalYear / 100
       return (1.05836E+4 - (1014.41 * u) + (33.78311 * (u ^ 2)) - (5.952053 * (u ^ 3)) - (0.1798452 * (u ^ 4)) + (0.022174192 * (u ^ 5)) + (0.0090316521 * (u ^ 6)))
   else if aDecimalYear < 1600 then
       set u to (aDecimalYear - 1000) / 100
       return (1574.2 - (556.01 * u) + (71.23472 * (u ^ 2)) + (0.319781 * (u ^ 3)) - (0.8503463 * (u ^ 4)) - (0.005050998 * (u ^ 5)) + (0.0083572073 * (u ^ 6)))
   else if aDecimalYear < 1700 then
       set u to aDecimalYear - 1600
       return (120 - (0.9808 * u) - (0.01532 * (u ^ 2)) + ((u ^ 3) / 7129))
   else if aDecimalYear < 1800 then
       set u to aDecimalYear - 1700
       return (8.83 + (0.1603 * u) - (0.0059285 * (u ^ 2)) + (1.3336E-4 * (u ^ 3)) - ((u ^ 4) / 1174000))
   else if aDecimalYear < 1860 then
       set u to aDecimalYear - 1800
       return (13.72 - (0.332447 * u) + (0.0068612 * (u ^ 2)) + (0.0041116 * (u ^ 3)) - (3.7436E-4 * (u ^ 4)) + (1.21272E-5 * (u ^ 5)) - (1.69E-7 * (u ^ 6)) + (8.75E-10 * (u ^ 7)))
   else if aDecimalYear < 1900 then
       set u to aDecimalYear - 1860
       return (7.62 + (0.5737 * u) - (0.251574 * (u ^ 2)) + (0.01680668 * (u ^ 3)) - (4.473624E-4 * (u ^ 4)) + ((u ^ 5) / 233174))
   else if aDecimalYear < 1920 then
       set u to aDecimalYear - 1900
       return (-2.79 + (1.494119 * u) - (0.0598939 * (u ^ 2)) + (0.0061966 * (u ^ 3)) - (1.97E-4 * (u ^ 4)))
   else if aDecimalYear < 1941 then
       set u to aDecimalYear - 1920
       return (21.2 + (0.84493 * u) - (0.0761 * (u ^ 2)) + (0.0020936 * (u ^ 3)))
   else if aDecimalYear < 1961 then
       set u to aDecimalYear - 1950
       return (29.07 + (0.407 * u) - ((u ^ 2) / 233) + ((u ^ 3) / 2547))
   else if aDecimalYear < 1986 then
       set u to aDecimalYear - 1975
       return (45.45 + (1.067 * u) - ((u ^ 2) / 260) - ((u ^ 3) / 718))
   else if aDecimalYear < 2005 then
       set u to aDecimalYear - 2000
       return (63.86 + (0.3345 * u) - (0.060374 * (u ^ 2)) + (0.0017275 * (u ^ 3)) + (6.51814E-4 * (u ^ 4)) + (2.373599E-5 * (u ^ 5)))
   else if aDecimalYear < 2050 then
       set u to aDecimalYear - 2000 # (correct!)
       return (62.92 + (0.32217 * u) + (0.005589 * (u ^ 2)))
   else if aDecimalYear < 2150 then
       return (-20 + (32 * (((aDecimalYear - 1820) / 100) ^ 2 - (0.5628 * (2150 - aDecimalYear)))))
   else
       set u to (aDecimalYear - 1820) / 100
       return (-20 + (32 * (u ^ 2)))
   end if
end deltaT

to JDN(_year, _month, _day, _hour, _min, _sec)
   # The Julian Date (JD) of any instant is the Julian day number
   # for the preceding noon plus the fraction of the day since
   # that instant. Julian Dates are expressed as a Julian day
   # number with a decimal fraction added. For example, The Julian
   # Date for 00:30:00.0 1 January 2013 is 2456293.520833.
   local a, y, m, _JDN, time_fraction
   set a to ((14 - _month) div 12)
   set y to _year + 4800 - a
   set m to _month + (12 * a) - 3
   # Gregorian Calendar starts October 15, 1582
   if _year > 1582 then
       set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
   else if _year = 1582 and _month > 10 then
       set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
   else if _year = 1582 and _month = 10 and _day ≥ 15 then
       set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
   else
       # Julian calendar before oct 15, 1582
       set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - 32083
   end if
   set _hour to (_hour - 12) / 24
   if _min ≠ 0 then
       set _min to _min / 1440 # minutes per day
   end if
   
   if _sec ≠ 0 then
       set _sec to _sec / 86400 # seconds per day
   end if
   
   set _JDN to _JDN + _hour + _min + _sec
   
   return _JDN
end JDN

to gregorianCalendarTriplet from JDN
   # This is an algorithm by Richards to convert a Julian Day Number, J, to a date in the Gregorian calendar
   # (proleptic, when applicable). Richards does not state which dates the algorithm is valid for.
   # it seems to work for dates after, BUT NOT INCLUDING 15th of oct 1582.
   # I have extended it to return a hour, min, sec triplet for the fractional part of the JDN.
   local y, m, n, r, p, v, s, w, b, c, f, g, h, d, MT, yr, JD
   set y to 4716
   set j to 1401
   set m to 2
   set n to 12
   set r to 4
   set p to 1461
   set v to 3
   set u to 5
   set s to 153
   set w to 2
   set b to 274277
   set c to -38
   
   set JD to JDN as integer
   set f to JD + j + (((4 * JD + b) div 146097) * 3) div 4 + c
   set e to r * f + v
   set g to (e mod p) div r
   set h to u * g + w
   set d to (h mod s) div u + 1
   set MT to (h div s + m) mod n + 1
   set yr to e div p - y + (n + m - MT) div n
   return {yr, MT, d}
end gregorianCalendarTriplet

Last edited by McUsrII (2015-05-20 04:17:43 pm)


Filed under: MoonPhases

Offline

 

Board footer

Powered by FluxBB

RSS (new topics) RSS (active topics)