Astronomy: Lunar Eclipse V0.2 PLATFORMS Numworks BBC microbit-serial ok
"""ALE Astronomy: Lunar Eclipse """ __version__="0.2" crid="Astronomy:Lunar Eclipse V"+__version__+""" © 2023 SteveG1CMZ""" from math import * #customise showme=True #select which edition of jmaa to use edition1=0 #jmaa edition 1 (includes 1-errata) edition2=1 #jmaa edition 2 (includes 2-errata) #end custom rads=pi/180 degs=180/pi # print("Astronomical Algorithms: Meeus") def date(JD): #07. JM: Valid for positive JD JD5=JD+0.5 Z=int(JD5); F=JD5-Z if Z<2299161: A=Z if Z>=2299161: alpha=int( (Z-1867216.25)/36524.25 ) A=Z+1+alpha-int(alpha/4) B=A+1524 C=int( (B-122.1)/365.25 ) D=int( 365.25*C ) E=int( (B-D)/30.6001 ) DD=B-D-int(30.6001*E) + F if E<14: MM=E-1 if E>=14: MM=E-13 #14 or 15 if MM>2: YYYY=C-4716 if MM<=2: YYYY=C-4715 #1 or 2 return [YYYY, MM, DD] #DD includes decimals def DTm(JD): YMD=date(JD) DD=YMD[2] DDi=int(DD) HH=(DD-DDi)*24 HHi=int(HH) MM=(HH-HHi)*60 MM=round(MM) #MMi=int(MM) #SS=(MM-MMi)*60 return [YMD[0],YMD[1],int(YMD[2]),HHi,MM] #45=47 position of the moon def position_E06(T): #45=47.6 Enhancement for varying eccentricity E45=1 -0.002516*T -0.0000074*T**2 return E45 #47=49 Mean Phases of the Moon def phases_k02(year): #47=49.2 Approximate k from real year k=(year-2000)*12.3685 return k def phases_T03(k): #47=49.3 Approximate T (centuries since 2000.0) T=k/1236.85 return T def phases_meanm01(k,T): #47=49.01 mean phases of Moon JDE=2451550.09765 +29.530588853*k +(edition1*0.0001337 +edition2*0.00015437)*T**2 -0.000000150*T**3 +0.00000000073*T**4 return JDE def phases_sma04(k,T): #47=49.4 Suns Mean Anomaly if edition1: M = 2.5534 +29.10535669*k -0.0000218*T**2 -0.00000011*T**3 if edition2: M =2.5534 + 29.10535670*k -0.0000014*T**2 -0.00000011*T**3 return reduce(M,360) def phases_mma05(k,T): #47=49.5 Moons mean anomaly MP=201.5643 +385.81693528*k +(edition1*0.0107438+edition2*0.107582)*T**2 +(edition1*0.00001239+edition2*0.00001238)*T**3 -0.0000000058*T**4 return reduce(MP,360) def phases_mal06(k,T): #47=49.6 Moons argument of latitude if edition1: F= 160.7108 +390.67050274*k -0.0016341*T**2 -0.00000227*T**3 +0.000000011*T**4 if edition2: F= 160.7108 + 390.67050284*k -0.0016118*T**2 -0.00000227*T**3 +0.000000011*T**4 return reduce(F,360) def phases_long07(k,T): #47=49.7 Long of the ascending node of lunar orbit if edition1: OMEGA=124.7746 -1.56375580*k +0.0020691*T**2 +0.00000215*T**3 if edition2: OMEGA=124.7746 -1.56375588*k +0.0020672*T**2 +0.00000215*T**3 return reduce(OMEGA,360) #52=54 eclipses def eclipses_adjdays01( lunar,M,MP,F1,A1,OMEGA,E45,indegrees): """#52=54.1 adjustment in days (Adjustments are typically hours). In 1950-2050 JM reports around 1 minute accuracy. #matches jm examples 52.c 52.d """ if indegrees: return eclipses_adjdays01(lunar, M*rads,MP*rads,F1*rads,A1*rads,OMEGA*rads,E45, False) if lunar: days= -0.4065*sin(MP) + 0.1727*E45*sin(M) else: #solar days= -0.4075*sin(MP) + 0.172*E45*sin(M) days=days + 0.0161*sin(2*MP) -0.0097*sin(2*F1) +0.0073*E45*sin(MP-M) -0.0050*E45*sin(MP+M) -0.0023*sin(MP-2*F1) +0.0021*E45*sin(2*M) +0.0012*sin(MP+2*F1) +0.0006*E45*sin(2*MP+M) -0.0004*sin(3*MP) -0.0003*E45*sin(M+2*F1) +0.0003*sin(A1) -0.0002*E45*sin(M-2*F1) -0.0002*E45*sin(2*MP-M) -0.0002*sin(OMEGA) #print("daysadj:",days) return days # def eclipses_P(M,MP,F1,E45,indegrees): #52=54.1+ if indegrees: return eclipses_P(M*rads,MP*rads,F1*rads,E45,False) P=0.2070*E45*sin(M) +0.0024*E45*sin(2*M) -0.0392*sin(MP) +0.0116*sin(2*MP) -0.0073*E45*sin(MP+M) +0.0067*E45*sin(MP-M) +0.0118*sin(2*F1) return P def eclipses_Q(M,MP,E45,indegrees): #52=54.1+ if indegrees: return eclipses_Q(M*rads,MP*rads,E45,False) Q=+5.2207 -0.0048*E45*cos(M) +0.0020*E45*cos(2*M) -0.3299*cos(MP) -0.0060*E45*cos(MP+M) +0.0041*E45*cos(MP-M) return Q def eclipses_u(M,MP,E45,indegrees): #52=54.1+ if indegrees: return eclipses_u(M*rads,MP*rads,E45,False) u=0.0059 +0.0046*E45*cos(M) -0.0182*cos(MP) +0.0004*cos(2*MP) -0.0005*cos(M+MP) return u def lunareclipse(year=None,koff=-0.5): """Look for a lunar eclipse near year. year: real (eg 1965.3 for March) koff: 0 solar 0.5 or -0.5 lunar FP(k): 0 New Moon 0.25 1stQ 0.5 Full 0.75 3rdQ k=0=>2000 Jan6 New Moon Note: JM add 0.5 but example subtracts 0.5 Accuracy: Peak Time JDE within a minute, UTC 2mins """ if year==None: year=float(input( "Lunar eclipse near YYYY.?")) if False: #show input print("Lunar eclipse near",year) #----- #47=49.02 estimate k from real year k=phases_k02(year) k=floor(k)+koff #k=-328.5 # test value #47=49.03 estimate T centuries from k T=phases_T03(k) #47.01=49.01 estimate JDE JDE=phases_meanm01(k,T) #47=49.04 Suns mean anomaly M= phases_sma04(k,T) #47=49.05 Moons mean anomaly MP=phases_mma05(k,T) #47=49.06 Moons argument of latitude F=phases_mal06(k,T) #if abs(sin(F*rads))>0.36: # print("no eclipse") #else eclipse possible... #47.07 Long of the ascending node of lunar orbit OMEGA=phases_long07(k,T) #45=47.6 enhancement E45=position_E06(T) #52=54 eclipses (note only A1 not 47=49 A1..A14) F1=F-0.02665*sin(OMEGA*rads) #degrees A1=299.77 +0.107408*k -0.009173*T**2 #degrees #52.1 days adjustment JDE0=JDE adjustdays=eclipses_adjdays01( True,M,MP,F1,A1,OMEGA,E45,True) JDE=JDE+adjustdays #JDE is now the peak time of a potential eclipse P=eclipses_P(M,MP,F1,E45,True) Q=eclipses_Q(M,MP,E45,True) F1rads=F1*rads W=abs(cos(F1rads)) GAMMA=(P*cos(F1rads)+Q*sin(F1rads)) * (1 -0.0048*W) #units:eqradius u=eclipses_u(M,MP,E45,True) #magnitude of eclipse (<0: no eclipse) #note: fraction illuminated not brightness #52=54.3 magecl_penumbral=(1.5573+u-abs(GAMMA))/0.5450 #52=54.4 magecl_umbral=(1.0128-u-abs(GAMMA))/0.5450 n=0.5458 +0.0400*cos(MP*rads) if magecl_umbral<0: #no umbral eclipse partialumbralsdur=0 #semi: minutes totalumbralsdur=0 #semi: minutes else: #reuse P,T PART=1.0128-u TOTAL=0.4678-u #sqrt needs guards when not total num=max(PART**2-GAMMA**2,0) partialumbralsdur=(60/n)* sqrt(num) num=max(TOTAL**2-GAMMA**2,0) totalumbralsdur=(60/n)* sqrt(num) if magecl_penumbral<0: #no penumbral eclipse penumbralsdur=0 #semi: minutes else: h=1.5573 +u #edition1:H edition2:h num=max(h**2-GAMMA**2,0) #ditto penumbralsdur=(60/n)*sqrt(num) #print("umbral:partsemi: ",partialumbralsdur) #print("umbral:totalsemi:",totalumbralsdur) #print("penumbral:semi: ",penumbralsdur) if magecl_umbral>=0: if magecl_umbral>=1: itsa="total umbral" else: itsa="partial umbral" else: if magecl_penumbral>=0: itsa="penumbral" else: itsa="None" JDE=None if itsa=="None": return None if True: #print print(itsa+" lunar eclipse:") print("magu:",round(magecl_umbral,3)) #magu if magecl_umbral>0 or magecl_penumbral>0: print("k:",k) print("JDE :",JDE) print(DTm(JDE)) #semidurations print("mins:", round(totalumbralsdur),round(partialumbralsdur),round(penumbralsdur),"(semi)" ) #GAMMA measures how central the eclipse is print("Gamma",round(GAMMA,4)) print("M MP F OMEGA degs:") print( round(M,3),round(MP,3),round(F,3),round(OMEGA,3)) #return mins tbd return [itsa,JDE,magecl_umbral,GAMMA] #supporting math def reduce(degrees,upto): """Return degrees in 0..upto. """ a=reduces(degrees,upto) if a<0: a+=upto return a #0..upto def reduces(degrees,upto): """Return signed reduced degrees (or radians). """ a=(abs(degrees) % upto)*copysign(1,degrees) return a #-upto..upto #examples def examples(): #print("a few examples...") lunareclipse(1973.5,-0.5) #52=54.c lunareclipse(1997+9/12,-0.5) #52=54.d lunareclipse(2023+9/12,0.5) #None def calendar(yyyy,k): #year calendar (real) print("Calendar:",yyyy) for ii in range (12): i=float(ii) lunareclipse(yyyy+(i/12),k) if showme or __name__=="__main__": print(crid) #examples() calendar(2023.75,0.5) #print("OK") #Is source complete? NO! #we are at the Numworks limit #ok