Moon Phases

Just starting out? Need help? Post your questions and find answers here.
vwidmer
Enthusiast
Enthusiast
Posts: 282
Joined: Mon Jan 20, 2014 6:32 pm

Moon Phases

Post by vwidmer »

I am trying to convert this old basic code:

Code: Select all

Declare SUB FLMOON (N!, NPH!, JD&, FRAC!)
Declare SUB CALDAT (JULIAN&, MM!, ID!, IYYY!)
Declare FUNCTION JULDAY& (IM!, ID!, IY!)

'PROGRAM D1R1
'Driver for routine FLMOON
CLS
TZONE = 5!
Dim PHASE$(4), TIMSTR$(2)
For i = 1 To 4
  Read PHASE$(i)
Next i
Data "new moon", "first quarter", "full moon", "last quarter"
For i = 1 To 2
  Read TIMSTR$(i)
Next i
Data " AM", " PM"
Print "Date of the next few phases of the moon"
Print "Enter today's date (e.g. 1,31,1982)"
INPUT IM, ID, IY
Print
TIMZON = -TZONE / 24!
'Approximate number of full moons since January 1900
N = Int(12.37 * (IY - 1900 + (IM - 0.5) / 12!))
NPH = 2
J1& = JULDAY&(IM, ID, IY)
Call FLMOON(N, NPH, J2&, FRAC)
N = Int(N + (J1& - J2&) / 28!)
Print "   Date", "  Time(EST)", " Phase"
For i = 1 To 20
  Call FLMOON(N, NPH, J2&, FRAC)
  IFRAC = CInt(24! * (FRAC + TIMZON))
  If IFRAC < 0 Then
    J2& = J2& - 1
    IFRAC = IFRAC + 24
  End If
  If IFRAC >= 12 Then
    J2& = J2& + 1
    IFRAC = IFRAC - 12
  Else
    IFRAC = IFRAC + 12
  End If
  If IFRAC > 12 Then
    IFRAC = IFRAC - 12
    ISTR = 2
  Else
    ISTR = 1
  End If
  Call CALDAT(J2&, IM, ID, IY)
  Print USING; "##"; IM;
  Print USING; "###"; ID;
  Print USING; "#####"; IY;
  Print USING; "#########"; IFRAC;
  Print TIMSTR$(ISTR); "     "; PHASE$(NPH + 1)
  If NPH = 3 Then
    NPH = 0
    N = N + 1
  Else
    NPH = NPH + 1
  End If
Next i
End

Sub FLMOON(N, NPH, JD&, FRAC)
RAD = 0.017453293
C = N + NPH / 4!
T = C / 1236.85
T2 = T ^ 2
AQ = 359.2242 + 29.105356 * C
am = 306.0253 + 385.816918 * C + 0.01073 * T2
JD& = 2415020 + 28 * N + 7 * NPH
XTRA = 0.75933 + 1.53058868 * C + (0.0001178 - 0.000000155 * T) * T2
If NPH = 0 Or NPH = 2 Then
  XTRA = XTRA + (0.1734 - 0.000393 * T) * Sin(RAD * AQ) - 0.4068 * Sin(RAD
* am)
ElseIf NPH = 1 Or NPH = 3 Then
  XTRA = XTRA + (0.1721 - 0.0004 * T) * Sin(RAD * AQ) - 0.628 * Sin(RAD *
am)
Else
  Print "NPH is unknown."
  Exit Sub
End If
If XTRA >= 0! Then
  i = Int(XTRA)
Else
  i = Int(XTRA - 1!)
End If
JD& = JD& + i
FRAC = XTRA - i
End Sub

Function JULDAY&(MM, ID, IYYY)
IGREG& = 588829
If IYYY = 0 Then Print "There is no Year Zero.": Exit Function
If IYYY < 0 Then IYYY = IYYY + 1
If MM > 2 Then
  JY = IYYY
  JM = MM + 1
Else
  JY = IYYY - 1
  JM = MM + 13
End If
JD& = Int(365.25 * JY) + Int(30.6001 * JM) + ID + 1720995
If ID + 31 * (MM + 12 * IYYY) >= IGREG& Then
  JA = Int(0.01 * JY)
  JD& = JD& + 2 - JA + Int(0.25 * JA)
End If
JULDAY& = JD&
End Function

Sub CALDAT(JULIAN&, MM, ID, IYYY)
IGREG& = 2299161
If JULIAN& >= IGREG& Then
  JALPHA& = Int(((JULIAN& - 1867216) - 0.25) / 36524.25)
  JA& = JULIAN& + 1 + JALPHA& - Int(0.25 * JALPHA&)
Else
  JA& = JULIAN&
End If
JB& = JA& + 1524
JC& = Int(6680! + ((JB& - 2439870) - 122.1) / 365.25)
JD& = 365 * JC& + Int(0.25 * JC&)
JE& = Int((JB& - JD&) / 30.6001)
ID = JB& - JD& - Int(30.6001 * JE&)
MM = JE& - 1
If MM > 12 Then MM = MM - 12
IYYY = JC& - 4715
If MM > 2 Then IYYY = IYYY - 1
If IYYY <= 0 Then IYYY = IYYY - 1
End Sub 
to RB:

Code: Select all

ImportC ""
  floor.d   (value.d)
  ceil.d    (value.d)
EndImport

Declare.d julianday(day,month,year,hour,minute,second)
Procedure.d julianday(day,month,year,hour,minute,second)
  
  Protected a.d, b.d, taf.d, mo.d, jar.d,time.d, jul_tag.d
  
  taf = day
  mo  = month
  jar = year
  
  If mo < 3
    mo  = mo + 12
    jar = jar - 1
  EndIf
  a = 0
  b = 0
  
  If jar < 0 ;year B.C. If  historical, in astronomical counting jar+1 schould be skipt
    jar + 1
  EndIf
  
  If jar + mo / 100 + taf / 10000 >= 1582.1015
    a = IntQ(jar / 100)
    b = 2 - a + IntQ(a / 4)
  EndIf
  
  time= hour / 24 + minute / 1440 + second /86400 ;convert time 0-24 in range  0.00-1.00 
  
  jul_tag = IntQ(365.25 * (jar+ 4716 )) + IntQ(30.6001 * (mo + 1 )) + b -1524.5 + taf + time
  ProcedureReturn jul_tag
  
EndProcedure

Procedure flmoon (n, nph, Array r(1))
RAD.f = #PI/180.0
c.f = n + nph / 4
t.f = c / 1236.85
t2.f = Pow(t,2)
aq.f = 359.2242 + 29.105356 * c
am.f = 306.0253 + 385.816918 * c + 0.010730 * t2
xtra = 0.75933 + 1.53058868 * c + ((0.0001178) - (1.55e-7) * t) * t2
If nph = 0 Or nph = 2
  xtra = xtra + (0.1734 - 0.000393 * t) * Sin(RAD * aq) - 0.4068 * Sin(RAD * am)
ElseIf nph = 1 Or nph = 3
  xtra = xtra + (0.1721 - 0.0004 * t) * Sin(RAD * aq) - 0.628 * Sin(RAD * am)
Else
  Debug "nph is unknown"
  End
EndIf

If xtra>=0
  xtrai = floor(xtra)
Else
  xtrai = ceil(xtra - 1)
EndIf

mt = 2415020 + 28 * n + 7 * nph

mt = jd + xtrai
frac = xtra - xtrai

jd = julianday(1,1,2015,0,0,0)

r(0)=mt
r(1)=frac

Debug Mod((jd-mt +30),30)
Debug frac

EndProcedure

Dim a(1)

year=2015
month=01
nph.i=2     ;0=new moon, 1=first quarter, 2=full, 3=last quarter

n.i = floor(12.37 * (year -1900 + ((month - 0.5)/12.0)))

For nph = 0 To 3
  Debug "NPH: " + Str(nph)
  flmoon(n, nph, a())
Next 
it seems it will only calculate the full moon as shown above for the month of JAN 2015
I dont know what some of the things are in the old basic code like the ! etc..

Any one know what I am doing wrong why I can get the other moon phases?

The dates for JAN 2015 should be:
Full Moon: 5th
New Moon: 20th
First Quarter: 27th
Last Quarter: 13th

Thanks
WARNING: I dont know what I am doing! I just put stuff here and there and sometimes like magic it works. So please improve on my code and post your changes so I can learn more. TIA
User avatar
Demivec
Addict
Addict
Posts: 4091
Joined: Mon Jul 25, 2005 3:51 pm
Location: Utah, USA

Re: Moon Phases

Post by Demivec »

Here's my initial run at translating the original code:

Code: Select all

Procedure FLMOON(N, NPH, *JD.INTEGER, *FRAC.FLOAT)
  Protected RAD.f, C, T.f, T2.f, AQ.f, AM.f, XTRA.f
  
  RAD = 0.017453293
  C = N + NPH / 4
  T = C / 1236.85
  T2 = Pow(T, 2)
  AQ = 359.2242 + 29.105356 * C
  AM = 306.0253 + 385.816918 * C + 0.01073 * T2
  *JD\i = 2415020 + 28 * N + 7 * NPH
  XTRA = 0.75933 + 1.53058868 * C + (0.0001178 - 0.000000155 * T) * T2
  
  If NPH = 0 Or NPH = 2
    XTRA = XTRA + (0.1734 - 0.000393 * T) * Sin(RAD * AQ) - 0.4068 * Sin(RAD * AM)
  ElseIf NPH = 1 Or NPH = 3
    XTRA = XTRA + (0.1721 - 0.0004 * T) * Sin(RAD * AQ) - 0.628 * Sin(RAD * AM)
  Else
    PrintN("NPH is unknown.")
    ProcedureReturn
  EndIf
  
  If XTRA >= 0
    i = Int(XTRA)
  Else
    i = Int(XTRA - 1)
  EndIf
  
  *JD\i + i
  *FRAC\f = XTRA - i
EndProcedure

Procedure JULDAY(MM, ID, IYYY)
  Protected IGREG = 588829, JD, JM, JY
  
  If IYYY = 0
    PrintN("There is no Year Zero.")
    ProcedureReturn
  EndIf
  
  If IYYY < 0
    IYYY + 1
  EndIf
  
  If MM > 2
    JY = IYYY
    JM = MM + 1
  Else
    JY = IYYY - 1
    JM = MM + 13
  EndIf
  
  JD = Int(365.25 * JY) + Int(30.6001 * JM) + ID + 1720995
  If ID + 31 * (MM + 12 * IYYY) >= IGREG
    JA = Int(0.01 * JY)
    JD = JD + 2 - JA + Int(0.25 * JA)
  EndIf
  ProcedureReturn JD
EndProcedure

Procedure CALDAT(JULIAN, *MM.INTEGER, *ID.INTEGER, *IYYY.INTEGER)
  Protected IGREG = 2299161, JALPHA, JA, JB, JC, JD, JE
  
  If JULIAN >= IGREG
    JALPHA = Int(((JULIAN - 1867216) - 0.25) / 36524.25)
    JA = JULIAN + 1 + JALPHA - Int(0.25 * JALPHA)
  Else
    JA = JULIAN
  EndIf
      
  JB = JA + 1524
  JC = Int(6680 + ((JB - 2439870) - 122.1) / 365.25)
  JD = 365 * JC + Int(0.25 * JC)
  JE = Int((JB - JD) / 30.6001)
  *ID\i = JB - JD - Int(30.6001 * JE)
  *MM\i = JE - 1
  If *MM\i > 12
    *MM\i - 12
  EndIf
  
  *IYYY\i = JC - 4715
  If *MM\i > 2
    *IYYY\i - 1
  EndIf
  
  If *IYYY\i <= 0
    *IYYY\i - 1
  EndIf
  
EndProcedure

;PROGRAM D1R1
;Driver for routine FLMOON
If OpenConsole()
  DataSection
    Data.s "new moon", "first quarter", "full moon", "last quarter"
    Data.s " AM", " PM"
  EndDataSection
  
  Define TZONE, i
  ClearConsole()
    TZONE = 0 ;enter timezone difference from GMT in hours, value is negative if east of GMT, and positive for west of GMT
  Dim PHASE$(4)
  Dim TIMSTR$(2)
  For i = 1 To 4
    Read.s PHASE$(i)
  Next
  
  For i = 1 To 2
    Read.s TIMSTR$(i)
  Next
  
  Define InputLine$, IM, ID, IY
  Repeat  
    PrintN("Date of the next few phases of the moon")
    PrintN("Enter today's date (e.g. 1,31,1982)")
    InputLine$ = Input()
  Until CountString(InputLine$, ",") = 2
  PrintN("")
  IM = Val(Trim(StringField(InputLine$, 1, ",")))
  ID = Val(Trim(StringField(InputLine$, 2, ",")))
  IY = Val(Trim(StringField(InputLine$, 3, ",")))
  
  Define TIMZON.f, N, NPH, J1, J2, FRAC.f, IFRAC, ISTR
  
  TIMZON = -TZONE / 24
  
  ;Approximate number of full moons since January 1900
  N = Int(12.37 * (IY - 1900 + (IM - 0.5) / 12))
  NPH = 2
  J1 = JULDAY(IM, ID, IY)
  FLMOON(N, NPH, @J2, @FRAC)
  N = Int(N + (J1 - J2) / 28)
  
  PrintN(Space(3) + "Date" + Space(9) + "Time(EST)" + Space(4) + "Phase")
  For i = 1 To 20
    FLMOON(N, NPH, @J2, @FRAC)
    IFRAC = Round(24 * (FRAC + TIMZON), #PB_Round_Nearest)
    If IFRAC < 0
      J2 = J2 - 1
      IFRAC = IFRAC + 24
    EndIf
    
    If IFRAC >= 12
      J2 = J2 + 1
      IFRAC = IFRAC - 12
    Else
      IFRAC = IFRAC + 12
    EndIf
    
    If IFRAC > 12
      IFRAC = IFRAC - 12
      ISTR = 2
    Else
      ISTR = 1
    EndIf
    
    CALDAT(J2, @IM, @ID, @IY)
    Print(RSet(Str(IM), 2, " "))
    Print(RSet(Str(ID), 3, " "))
    Print(RSet(Str(IY), 5, " "))
    Print(RSet(Str(IFRAC), 9, " "))
    PrintN(TIMSTR$(ISTR) + "     " + PHASE$(NPH + 1))
    
    If NPH = 3
      NPH = 0
      N + 1
    Else
      NPH + 1
    EndIf
  Next
  
  Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()
  CloseConsole()
EndIf
With the TZONE = 0 for GMT it generates these values:

Code: Select all

Date of the next few phases of the moon
Enter today's date (e.g. 1,31,1982)
1,1,2015

   Date         Time(EST)    Phase
12  6 2014       12 AM     full moon
12 13 2014        6 PM     last quarter
12 22 2014        2 AM     new moon
12 29 2014        7 AM     first quarter
 1  5 2015        2 AM     full moon
 1 12 2015        7 AM     last quarter
 1 20 2015        2 PM     new moon
 1 27 2015        5 PM     first quarter
 2  3 2015        2 PM     full moon
 2 10 2015        5 PM     last quarter
 2 19 2015        0 AM     new moon
 2 26 2015        1 AM     first quarter
 3  5 2015        0 AM     full moon
 3 12 2015        1 AM     last quarter
 3 20 2015       10 AM     new moon
 3 27 2015        8 AM     first quarter
 4  3 2015       10 AM     full moon
 4 10 2015        8 AM     last quarter
 4 18 2015        7 PM     new moon
 4 25 2015        3 PM     first quarter
Looks like the formula needs to be tuned up a little bit. The hour of each quarter is not precise enough and that wobbles the data a little as a result.
User avatar
glomph
User
User
Posts: 48
Joined: Tue Apr 27, 2010 1:43 am
Location: St. Elsewhere / Germany
Contact:

Re: Moon Phases

Post by glomph »

With so little code it is impossible to get a better result.
If you want to go deeper have a look here:
http://www.naughter.com/aa.html
vwidmer
Enthusiast
Enthusiast
Posts: 282
Joined: Mon Jan 20, 2014 6:32 pm

Re: Moon Phases

Post by vwidmer »

@Demivec: Thanks I learned a few things from your code.

@glomph: I see the same formula above around done in a few different ways. I will look at the code for what you sent and see if I can make heads and tails from it.

Thanks
WARNING: I dont know what I am doing! I just put stuff here and there and sometimes like magic it works. So please improve on my code and post your changes so I can learn more. TIA
WilliamL
Addict
Addict
Posts: 1224
Joined: Mon Aug 04, 2008 10:56 pm
Location: Seattle, USA

Re: Moon Phases

Post by WilliamL »

Here's a really simple example of finding the moon phases.

https://www.subsystems.us/uploads/9/8/9 ... nphase.pdf <if you look at the article you will see where I got my code

I've just got it working and only tried it for this date. If you find a problem I will post the correction.

Code: Select all

Procedure PDFMoonCalc(mnth,dy,yr)
    Define.f a,b,c,e,f,jd,totaldays,mooncycles,newmoonday,fullmoonday,nextnewmoonday
    Define daysinmonth,newmoonmonth,newmoonyear,fullmoonmonth,fullmoonyear,nextnewmoonmonth,nextnewmoonyear
    
    daysinmonth=DaysInMonth(newmoonmonth,yr)
    a=yr/100 : Debug a ; 20.2
    b=a/4    : Debug b ; 5.0500
    c=2-a+b  : Debug c ; -13.150
    e=365.25*(yr+4716) : Debug e ; 2460324.0
    f=30.6001*(mnth+1) : Debug f ; 275.4009
    jd=c+dy+e+f-1524.5  : Debug jd ; 2459072.75
    TotalDays=jd-2451549.5 : Debug totaldays ; 7523.25
    mooncycles=totaldays/29.53 : Debug mooncycles ; remainder is how far past new moon ; 254.766
    
    newmoonmonth=mnth : newmoonyear=yr
    ;newmoonday=dy-(29.53*(mooncycles-Int(mooncycles))) ;: Debug newmoonday ; -11.6296
    newmoonday=dy-(29.53*(Mod(mooncycles,1))) : Debug newmoonday ; -11.6296
    If newmoonday<0
        newmoonmonth-1 : If newmoonmonth<1 : newmoonmonth=12 : newmoonyear-1 : EndIf
        Debug "Subtracted month from newmoonmonth"
        daysinmonth=DaysInMonth(newmoonmonth,newmoonyear)
        newmoonday=daysinmonth+newmoonday
    EndIf
    
    fullmoonmonth=newmoonmonth : fullmoonyear=newmoonyear
    fullmoonday=newmoonday+(29.53/2) ;: Debug fullmoonday ; =34.135 
    If fullmoonday>daysinmonth
        fullmoonmonth+1 : If fullmoonmonth>12 : fullmoonmonth=1 : fullmoonyear+1 : EndIf
        Debug "Added month to fullmoonmonth"
        fullMoonday-daysinmonth ; days in old month
        daysinmonth=DaysInMonth(fullmoonmonth,fullmoonyear) ; new month
    EndIf
    
    nextnewmoonmonth=newmoonmonth : nextnewmoonyear=newmoonyear
    nextnewmoonday=newmoonday+29.53 ;: Debug nextnewmoonday ; =49 
    If nextnewmoonday>daysinmonth
        nextnewmoonmonth+1 : If nextnewmoonmonth>12 : nextnewmoonmonth=1 : nextnewmoonyear+1 : EndIf
        Debug "Added month to nextnewmoonmonth"
        nextnewmoonday-daysinmonth
    EndIf
    
    MessageRequester("Moon phases in "+Str(yr),"New moon on "+Str(newmoonmonth)+"/"+Round(NewMoonday,#PB_Round_Nearest)+"/"+Str(newmoonyear)+Chr(13)+
            "Full moon on "+Str(FullmoonMonth)+"/"+Round(FullMoonday,#PB_Round_Nearest)+"/"+Str(fullmoonyear)+Chr(13)+
            "next new moon on "+Str(nextnewmoonMonth)+"/"+Round(nextnewMoonday,#PB_Round_Nearest)+"/"+Str(nextnewmoonyear))
EndProcedure

PDFMoonCalc(8,11,2020)
Last edited by WilliamL on Sun Aug 16, 2020 6:07 pm, edited 1 time in total.
MacBook Pro-M1 (2021), Sonoma 14.4.1, PB 6.10LTS M1
vwidmer
Enthusiast
Enthusiast
Posts: 282
Joined: Mon Jan 20, 2014 6:32 pm

Re: Moon Phases

Post by vwidmer »

Thanks thats working for me nicely.
WARNING: I dont know what I am doing! I just put stuff here and there and sometimes like magic it works. So please improve on my code and post your changes so I can learn more. TIA
dige
Addict
Addict
Posts: 1256
Joined: Wed Apr 30, 2003 8:15 am
Location: Germany
Contact:

Re: Moon Phases

Post by dige »

vwidmer wrote:Thanks thats working for me nicely.
DaysInMonth(newmoonmonth,yr) is not a function, array, list, map or macro - where can I find it?
"Daddy, I'll run faster, then it is not so far..."
vwidmer
Enthusiast
Enthusiast
Posts: 282
Joined: Mon Jan 20, 2014 6:32 pm

Re: Moon Phases

Post by vwidmer »

Code: Select all

Procedure DaysInMonth(mnth,yr)
    Protected days
    days=Val(Mid(" 312831303130313130313031",mnth*2,2))
    If (yr/4=Int(yr/4) And yr/100<>Int(yr/100)) Or (yr/400=Int(yr/400)) ; is leap year
        If mnth=2 : days=29 : EndIf
    EndIf
    ProcedureReturn days
EndProcedure
Let me know if that works for you

*I just tried this again strange it shows last full moon as 2020/11/29 but it should be 2020/11/30 I think????? So may want to check
WARNING: I dont know what I am doing! I just put stuff here and there and sometimes like magic it works. So please improve on my code and post your changes so I can learn more. TIA
User avatar
Kiffi
Addict
Addict
Posts: 1361
Joined: Tue Mar 02, 2004 1:20 pm
Location: Amphibios 9

Re: Moon Phases

Post by Kiffi »

my version:

Code: Select all

Procedure DaysInMonth(Month, Year)
  
  Protected d
  
  d = Date(Year, Month, 1, 0, 0, 0) ; CurrentMonth
  d = AddDate(d, #PB_Date_Month, 1) ; NextMonth
  d = AddDate(d, #PB_Date_Day, -1)  ; minus 1 day
  
  ProcedureReturn Day(d)
  
EndProcedure
Hygge
vwidmer
Enthusiast
Enthusiast
Posts: 282
Joined: Mon Jan 20, 2014 6:32 pm

Re: Moon Phases

Post by vwidmer »

Kiffi wrote:my version:

Code: Select all

Procedure DaysInMonth(Month, Year)
  
  Protected d
  
  d = Date(Year, Month, 1, 0, 0, 0) ; CurrentMonth
  d = AddDate(d, #PB_Date_Month, 1) ; NextMonth
  d = AddDate(d, #PB_Date_Day, -1)  ; minus 1 day
  
  ProcedureReturn Day(d)
  
EndProcedure
Do you show November 30th 2020 as full moon?
WARNING: I dont know what I am doing! I just put stuff here and there and sometimes like magic it works. So please improve on my code and post your changes so I can learn more. TIA
dige
Addict
Addict
Posts: 1256
Joined: Wed Apr 30, 2003 8:15 am
Location: Germany
Contact:

Re: Moon Phases

Post by dige »

@Kiffi: smart implementation *thumbsup*
Kiffi wrote:my version:

Code: Select all

Procedure DaysInMonth(Month, Year)
  
  Protected d
  
  d = Date(Year, Month, 1, 0, 0, 0) ; CurrentMonth
  d = AddDate(d, #PB_Date_Month, 1) ; NextMonth
  d = AddDate(d, #PB_Date_Day, -1)  ; minus 1 day
  
  ProcedureReturn Day(d)
  
EndProcedure
"Daddy, I'll run faster, then it is not so far..."
eck49
Enthusiast
Enthusiast
Posts: 153
Joined: Sat Nov 14, 2020 10:08 pm
Location: England

Re: Moon Phases

Post by eck49 »

OK for the present day, but if you want dates from a long time ago watch out for the Julian -> Gregorian calendar change.
In England that cut out 11 dates in September 1752. Elsewhere it would depend when the change was made.
Ubuntu 22.04 64-bit
Purebasic 6.00 (as of 5 Sep 2022)
(native tongue: English)
Olli
Addict
Addict
Posts: 1071
Joined: Wed May 27, 2020 12:26 pm

Re: Moon Phases

Post by Olli »

Hello, actually I can see Jupiter which is near Saturne. It is like that :

Code: Select all

.
 
  °
Saturne (.)
Jupiter (°)

At this time, orientation WSW in France. It will be for Spain, Morocco and Senegal in a few minutes !

On December, the 21st, their apparent distance will reach its minimum : maybe anyone could take a photography...

Thank you, anyway for these codes !

PS : I published a code 13 years ago. Test it.
Thanks to Dobro (Zorro)
Thanks to Huitbit (he gave a link to its own elder subject)
Thanks to Kernadec (he gave a list of the main changes (12) of the calendar in the History and in the several countries.
CalamityJames
User
User
Posts: 78
Joined: Sat Mar 13, 2010 4:50 pm

Re: Moon Phases

Post by CalamityJames »

I think anyone reading this thread would be interested in the algorithms offered by "Sky and Telescope" magazine back in the 1980s. They published Applesoft Basic programs for all sorts of thing including phases of the moon. Amazingly, the source code for these programs is still available (https://skyandtelescope.org/astronomy-r ... telescope/). Just in case no-one has an Apple ][ still operating you can test these programs here: https://www.calormen.com/jsbasic/. In the 1990s I found these programs and adapted some of them into Visual Basic 6. I have repeated the exercise here, adapting the Full and New moons program with the minimum possible changes to get them to work in Pure Basic. I have left the original lines in the remarks for each Pure Basic line. The output from this program is in Julian Days (and fractions) so I merged another "Sky and Telescope" program (something they made very easy) to convert the output to something understandable. I believe the forecast time is accurate to 2 minutes or so.

One bit of the submission was not provided by Sky and Telescope. It is easy to spot.

How many people remember when any program variable could be only 2 characters?

Code: Select all

;https://skyandtelescope.org/astronomy-resources/basic-programs-from-sky-telescope/ ; Astronomy programs in Applesoft Basic
;https://www.calormen.com/jsbasic/ ;run applesoft basic programs

Global OutPut$

Procedure.s DayFractToHm(DayFract.d)
  Protected hours.i, minutes.i, NbMinutes.i
  Protected divisor.i, remainder.i
  
  NbMinutes = Round(DayFract * 60 * 24, #PB_Round_Nearest)
  divisor = 60
  hours = NbMinutes / divisor
  remainder = NbMinutes % divisor
  divisor / 60 ; seconds in a minute
  minutes = remainder / divisor
  
  ProcedureReturn RSet(Str(hours), 2, "0") + ":" + RSet(Str(minutes), 2, "0")
EndProcedure

 ;New and full moons                        ; 10 REM    NEW And FULL MOONS
                                            ; 12 REM
                                            ; 14 REM
R1.d=3.14159265/180: U.i=0                  ; 16 R1=3.14159265/180: U=0
Y.i = Val(InputRequester("New and full moons", 
 "Enter the year:", "2020"))                ; 18 INPUT "YEAR ";Y
                                            ; 20 PRINT
K0.d=Int((Y-1900)*12.3685)                  ; 22 K0=Int((Y-1900)*12.3685)
T.d=(Y-1899.5)/100                          ; 24 T=(Y-1899.5)/100
T2.d=T*T: T3.d=T*T*T                        ; 26 T2=T*T: T3=T*T*T
J0.d=2415020+29*K0.d                        ; 28 J0=2415020+29*K0
F0.d=0.0001178*T2-0.000000155*T3            ; 30 F0=0.0001178*T2-0.000000155*T3
F0=F0+0.75933+0.53058868*K0                 ; 32 F0=F0+0.75933+0.53058868*K0
F0=F0-0.000837*T-0.000335*T2                ; 34 F0=F0-0.000837*T-0.000335*T2
J.i=J+Int(F.d): F=F-Int(F)                  ; 36 J=J+Int(F): F=F-Int(F)
M0.d=K0*0.08084821133                       ; 38 M0=K0*0.08084821133

M0=360*(M0-Int(M0))+359.2242                ; 40 M0=360*(M0-Int(M0))+359.2242
M0=M0-0.0000333*T2                          ; 42 M0=M0-0.0000333*T2
M0=M0-0.00000347*T3                         ; 44 M0=M0-0.00000347*T3
M1.d=K0*0.07171366128                       ; 46 M1=K0*0.07171366128
M1=360*(M1-Int(M1))+306.0253                ; 48 M1=360*(M1-Int(M1))+306.0253
M1=M1+0.0107306*T2                          ; 50 M1=M1+0.0107306*T2
M1=M1+0.00001236*T3                         ; 52 M1=M1+0.00001236*T3
B1.d=K0*0.08519585128                       ; 54 B1=K0*0.08519585128
B1=360*(B1-Int(B1))+21.2964                 ; 56 B1=360*(B1-Int(B1))+21.2964
B1=B1-0.0016528*T2                          ; 58 B1=B1-0.0016528*T2
B1=B1-0.00000239*T3                         ; 60 B1=B1-0.00000239*T3
For K9.i=0 To 28                            ; 62 For K9=0 To 28
  J=J0+14*K9: F=F0+0.765294*K9              ; 64 J=J0+14*K9: F=F0+0.765294*K9
  K.d=K9/2                                  ; 66 K=K9/2
  M5.d=(M0+K*29.10535608)*R1                ; 68 M5=(M0+K*29.10535608)*R1
  M6.d=(M1+K*385.81691806)*R1               ; 69 M6=(M1+K*385.81691806)*R1
  B6.d=(B1+K*390.67050646)*R1               ; 70 B6=(B1+K*390.67050646)*R1
  F=F-0.4068*Sin(M6)                        ; 71 F=F-0.4068*Sin(M6)
  F=F+(0.1734-0.000393*T)*Sin(M5)           ; 72 F=F+(0.1734-0.000393*T)*Sin(M5)
  F=F+0.0161*Sin(2*M6)                      ; 73 F=F+0.0161*Sin(2*M6)
  F=F+0.0104*Sin(2*B6)                      ; 74 F=F+0.0104*Sin(2*B6)
  F=F-0.0074*Sin(M5-M6)                     ; 75 F=F-0.0074*Sin(M5-M6)
  F=F-0.0051*Sin(M5+M6)                     ; 76 F=F-0.0051*Sin(M5+M6)
  F=F+0.0021*Sin(2*M5)                      ; 77 F=F+0.0021*Sin(2*M5)
  F=F+0.0010*Sin(2*B6-M6)                   ; 78 F=F+0.0010*Sin(2*B6-M6)
  J=J+Int(F): F=F-Int(F)                    ; 82 J=J+Int(F): F=F-Int(F)
  If U=0 :OutPut$="New Moon  : ":EndIf      ; 84 If U=0 THEN PRINT " NEW MOON ";
  If U=1 :OutPut$="Full Moon: ":EndIf       ; 86 If U=1 THEN PRINT "FULL MOON ";
  OutPut$+Chr(9)+Str(J)+Chr(9)+StrD(F, 5)   ; 88 PRINT J;F
  U=U+1: If U=2: U=0 :EndIf                 ; 90 U=U+1: If U=2 THEN U=0
                                            ; 92 Next
                                            ; 94 End
                                            ; 95 REM  ------------------------
                                            ; 96 REM  APPEARED IN ASTRONOMICAL
                                            ; 97 REM  COMPUTING, SKY & TELE-
                                            ; 98 REM  SCOPE, MARCH, 1985
                                            ; 99 REM  ------------------------
  
  
  ; Julian Day to Calendar Date             ; 900 REM   JD --> CALENDAR
                                            ; 905 REM
                                            ; 910 REM INPUT "J,F ";J,F
  G.i = 1                                   ; 915 INPUT "JC (0) OR GC (1) ";G
  F=F+0.5                                   ; 920 F=F+0.5
  If F<1:Goto L935:EndIf                    ; 925 If F<1 THEN 935
  F=F-1: J=J+1                              ; 930 F=F-1: J=J+1
  L935:
  If G=1: Goto  L945:EndIf                  ; 935 If G=1 THEN 945
  A.i=J: Goto L955                          ; 940 A=J: Goto 955
  L945:
  A1.i=Int((J/36524.25)-51.12264)             ; 945 A1=Int((J/36524.25)-51.12264)
  A=J+1+A1-Int(A1/4)                        ; 950 A=J+1+A1-Int(A1/4)
  L955:
  B.i=A+1524                                ; 955 B=A+1524
  C.i=Int((B/365.25)-0.3343)                ; 960 C=Int((B/365.25)-0.3343)
  D.i=Int(365.25*C)                         ; 965 D=Int(365.25*C)
  E.d=Int((B-D)/30.61)                      ; 970 E=Int((B-D)/30.61)
  D=Int(B-D-Int(30.61*E)+F)                 ; 975 D=B-D-Int(30.61*E)+F
  M.d=E-1: Y=C-4716                         ; 980 M=E-1: Y=C-4716
  If E>13.5: M=M-12:EndIf                   ; 985 If E>13.5 THEN M=M-12
  If M<2.5: Y=Y+1:EndIf                     ; 990 If M<2.5 THEN Y=Y+1
  OutPut$ + Chr(9) + Chr(9) + RSet(Str(D), 2,"0") +"/" + RSet(Str(M), 2,"0") + "/" + RSet(Str(Y), 4,"0") + Chr(9) + "  " + DayFractToHm(F)
  Debug OutPut$                             ; 995 PRINT "DATE: ";Y;M;D
Next
End                                         ; 997 End
                                            ; 1000 REM  ------------------------
                                            ; 1001 REM  APPEARED IN ASTRONOMICAL
                                            ; 1002 REM  COMPUTING, SKY & TELE-
                                            ; 1003 REM  SCOPE, MAY, 1984
                                            ; 1004 REM  ------------------------


  
Post Reply