REM Global variables declaration
OPTION EXPLICIT
REM tz is the time zone (e.g. 2.0 for UTC+2)
REM moonphase is the percentage of the moon disk that is lit
REM moonage is the number of days since last new moon
DIM FLOAT tz,moonphase,moonage

REM Constants
CONST TO_DEG=57.2957795131
CONST TO_RAD=1.74532925199e-2

tz=2.0
ComputeMoonPhase()
PRINT "Moon: "+STR$(moonphase*100.0,2,0)+"% "+STR$(moonage,2,1)+"d"

REM Computes the number of Julian days elapsed since 2000-01-01 12:00 UTC for
REM date d$ and time t$ (leave the latter empty to get the elapsed time in full
REM days). Note: t$ is the local time (the time zone is subtracted from it).
FUNCTION JD2000(d$,t$)
LOCAL FLOAT jd
LOCAL INTEGER day,month,year,hour,min,sec
day=VAL(LEFT$(d$,2)):month=VAL(MID$(d$,4,2)):year=VAL(RIGHT$(d$,4))-2000
jd=365*year+(year+3)\4+275*month\9-((month+9)\12)*(1+((year MOD 4)+2)\3)+day-31
IF LEN(t$)=8 THEN
hour=VAL(LEFT$(t$,2))-tz:min=VAL(MID$(t$,4,2)):sec=VAL(RIGHT$(t$,4))
jd=jd+hour/24.0+min/1440.0+sec/3600.0
ENDIF
JD2000=jd
END FUNCTION

SUB ComputeMoonPhase()
LOCAL FLOAT jd,N,M,Ec,Ls,ml,MM,MN,Ev,Ae,A3,MmP,mEc,A4,lP,V,lPP,NP,y,x
LOCAL FLOAT Lm,BM,delta
REM Constants for the Sun and the Moon at julian date 1980-01-00
CONST MAfactor=0.9856473321
CONST Elonge=278.83354
CONST Elongp=282.596403
CONST Eccent=0.016718
CONST Eccent2=1.016860112
CONST Mlfactor=13.1763966
CONST Mmlong=64.975464
CONST MMfactor=0.1114041
CONST Mmlongp=349.383063
CONST Mlnode=151.950429
CONST MNfactor=0.0529539
CONST Evfactor=1.2739
CONST Aefactor=0.1858
CONST A3factor=0.37
CONST mEcfactor=6.2886
CONST A4factor=0.214
CONST Vfactor=0.6583
CONST NPfactor=0.16
CONST CosMinc=0.995970321:REM COS(TO_RAD*5.145396)
CONST SinMinc=8.9683442e-2:REM SIN(TO_RAD*5.145396)
CONST Synmonth=29.53058868
REM Was 1e-6 in initial algorithm, but we only have simple precision floats
CONST EPSILON=1e-5
REM 7306=days elapsed between 2000-01-01 and 1980-01-00=1979-12-31
jd=JD2000(DATE$,TIME$)+7306
N=MAfactor*jd:N=N-N\360*360.0
M=N+Elonge-Elongp:M=M-M\360*360.0:M=TO_RAD*M
Ec=M
DO
delta=Ec-Eccent*SIN(Ec)-M
Ec=Ec-delta/(1-Eccent*COS(Ec))
LOOP UNTIL ABS(delta)<=EPSILON
Ec=Eccent2*TAN(Ec/2.0):Ec=2*TO_DEG*ATN(Ec)
Ls=Ec+Elongp:Ls=Ls-Ls\360*360.0
ml=Mlfactor*jd+Mmlong:ml=ml-ml\360*360.0
MM=ml-MMfactor*jd-Mmlongp:MM=MM-MM\360*360.0
MN=Mlnode-MNfactor*jd
Ev=Evfactor*SIN(TO_RAD*(2*(ml-Ls)-MM))
Ae=Aefactor*SIN(M):A3=A3factor*SIN(M)
MmP=MM+Ev-Ae-A3
mEc=mEcfactor*SIN(TO_RAD*MmP):A4=A4factor*SIN(TO_RAD*2*MmP)
lP=ml+Ev+mEc-Ae+A4
V=Vfactor*SIN(TO_RAD*2*(lP-Ls)):lPP=lP+V
NP=MN-NPfactor*SIN(M)
y=SIN(TO_RAD*(lPP-NP))*CosMinc:x=COS(TO_RAD*(lPP-NP))
Lm=TO_DEG*ATAN2(y,x)+NP
BM=TO_DEG*ASIN(SIN(TO_RAD*(lPP-NP))*SinMinc)
moonage=lPP-Ls
moonage=moonage-moonage\360*360.0:IF moonage<0 THEN moonage=moonage+360.0
moonphase=0.5-COS(TO_RAD*moonage)/2.0
moonage=Synmonth*moonage/360.0
END SUB

Function MoonPhase(y As Integer,m As Integer,d As Integer) As Integer
'calculates the moon phase (0-7), accurate to 1 segment.
'0 => new moon, 4 => full moon, 7= last part of third quarter ( before new)

Local Float jd

If m<3 Then y=y-1:m=m+12

m=m+1

jd=(((365.25*y)+(30.6*m)+d)-694039.09)/29.53 'total days divided by the moon cycle
MoonPhase=Fix((jd-Int(jd))*8) Mod 7
End Function

> a=0.5
> ? a
0.5
>
> a=694039.09
> ? a
694039

option explicit
option default none
dim integer k, a%
dim phase$(7)=("New", "Waxing Crescent","First Quarter","Waxing Gibbous", "Full", "Waning Gibbous", "Third Quarter", "Waning Crescent")
for k=1 to 31:? "2018 May ",k,str$(moon_phase(2018,5,k,0,a%),3,2)+"% ",phase$(a%):next k
end
function Julian(year as integer,month as integer,day as float) as float
' Returns the number of julian days for the specified day.
local integer a,b=0,c,e
if month < 3) then
year=year-1
month =month+12
endif
if (year > 1582 OR (year = 1582 and month>10) or (year = 1582 and month=10 and day > 15)) then
a=year/100
b=2-a+a/4
endif
c = 365.25*year
e = 30.6001*(month+1)
Julian= b+c+e+day+1720992.5
end function

function sun_position(j as float) as float
local float n,x,e,l,dl,v
local float SMALL_FLOAT
local integer i
SMALL_FLOAT=2.718281828^(-12)
n=360/365.2422*j
i=n\360
n=n-i*360.0
x=n-3.762863
if x<0 then x = x+ 360
x=rad(x)
e = x
do
dl = e - 0.016718 * sin(e) - x
e = e - dl/ (1 - 0.016718 * cos(e))
loop while abs(dl)>=SMALL_FLOAT
v=360/PI*atn(1.01686011182*tan(e/2))
l=v+282.596403
i=l\360
l=l-i*360.0
sun_position= l
end function

function moon_position(j as float, ls as float) as float

local float ms,l,mm,n,ev,sms,ae,ec
local integer i

ms = 0.985647332099 * j - 3.762863
if ms < 0 then ms =ms + 360.0
l = 13.176396 * j + 64.975464
i = l \ 360
l = l - i * 360.0
if l < 0 then l = l + 360.0
mm = l - 0.1114041 * j - 349.383063
i = mm \ 360
mm = mm - i * 360.0
n = 151.950429 - 0.0529539 * j
i = n \ 360
n = n- i * 360.0
ev = 1.2739 * sin(RAD((2*(l-ls)-mm)))
sms = sin(RAD(ms))
ae = 0.1858*sms
mm =mm+ ev-ae- 0.37*sms
ec = 6.2886*sin(RAD(mm))
l =l + ev + ec -ae + 0.214 * sin(RAD(2*mm))
l= 0.6583*sin(RAD(2*(l-ls)))+l
moon_position= l
end function

function moon_phase(year as integer,month as integer,day as integer, hour as float, ip as integer) as float


' Calculates more accurately than Moon_phase , the phase of the moon at
' the given epoch.
' returns the moon phase as a real number (0-1)
'

local float j= Julian(year,month,day+hour/24.0)-2444238.5
local float ls = sun_position(j)
local float lm = moon_position(j, ls)

local float t = lm - ls
if t < 0 then t = t + 360
ip = INT((t + 22.5)/45) AND 7
moon_phase = ((1.0 - cos(RAD(lm - ls)))/2)*100
end function

if moonage > 2.75 < 4 then

if moonage > 2.75 and moonage < 4 then

Select Case moonage
Case 2.75 To 4
...

Select Case moonage
Case 2.75 To 4
x$="Waxing Crescent"
Case ...
Case ...
etc...
End Select

Text 280,310,x$, LB,4, 1, RGB(white)

if moonage > 24 < 26 then
Text 280,310, "Waning Crescent" , LB,4, 1, RGB(white)
endif

if moonage > 26 < 28. then
Text 280,310, "Waning Crescent" , LB,4, 1, RGB(white)
endif
[/code]

what if moonage is exactly 26?

Edited by CaptainBoing 2018-05-06

' return days since 01-01-0001 without Julian calendar correction
FUNCTION Days(d$)
LOCAL d,m,y,a
d = VAL(MID$(d$,1,2))
m = VAL(MID$(d$,4,2))
y = VAL(MID$(d$,7,4))
a = INT((14-m)/12)
m = m+12*a
y = y-a
Days = 365*y+INT(y/4)-INT(y/100)+INT(y/400)+INT((153*m-457)/5)+d-306
END FUNCTION

' returns the moon phase for phase 0<=p<=4
FUNCTION MoonPhase$(p)
IF p < 0.1 THEN
MoonPhase$ = "New"
ELSEIF p < 0.9 THEN
MoonPhase$ = "Waxing Crescent"
ELSEIF p < 1.1 THEN
MoonPhase$ = "First Quarter"
ELSEIF p < 1.9 THEN
MoonPhase$ = "Waxing Gibbous"
ELSEIF p < 2.1 THEN
MoonPhase$ = "Full"
ELSEIF p < 2.9 THEN
MoonPhase$ = "Waning Gibbous"
ELSEIF p < 3.1 THEN
MoonPhase$ = "Third Quarter"
ELSEIF p < 3.9 THEN
MoonPhase$ = "Waning Crescent"
ELSE
MoonPhase$ = "New"
ENDIF
END FUNCTION

d = Days(d$)+19
p = d/29.530588
p = p-INT(p)
phase = 4*p
PRINT MoonPhase$(phase)