Astronomy: Eclipse V0.1
Algorithms based on Jean Meeus
Calculate eclipses (solar eclipse and lunar eclipse).
Call eclipse, or DTm for JDE to date.
See aeclipseuse for example calls.
Known bugs V0.1: “adjustdays(True “ always calculates adjustment for lunar eclipse. Workaround: “adjustdays(koff” will calculate the correct adjustment. (For the Oct 2023 solar eclipse, the mistake is less than 2 minutes)
This source is at the limit for Numworks source size.
Known bugs v0.1 The adjust JDE time calculation always performs the lunar eclipse adjustment, so my JDE time for solar eclipses will be incorrect. Workaround: change adjustdays(True to adjustdays(False when a solar eclipse is to be calculated V0.2 has been released which fixes that.
"""Astronomy: Eclipse """ __version__="0.1" crid="Astronomy:Eclipse V"+__version__+""" © 2023 SteveG1CMZ""" from math import * #customise #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. """ 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) 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 eclipses_solar_type(GAMMA,u): if abs(GAMMA)>1.5433+u: return None itsa="?";#print(GAMMA,u) if -0.9972<GAMMA<0.9972: itsa="central" if u<0: itsa="total central" if u>0.0047: itsa="annular central" if 0<u<0.0047: w=0.00462*sqrt(1-GAMMA**2) if u<w: itsa="annular-total (hybrid) central" else: itsa="annular central" if 0.9972<abs(GAMMA)<1.5433+u: itsa="non-central partial" #usually partial if 0.9972<abs(GAMMA)<0.9972+abs(u): #non-central total or annular itsa="non-central" #indistinguishable return itsa+" solar eclipse" def eclipse(year=None,koff=-0.5): """eclipse near year. year: real (eg 1965.25 for March) koff: 0 solar 0.5 or -0.5 lunar """ if year==None: year=float(input("Eclipse near YYYY.?")) #47=49.02 estimate k from real year k=phases_k02(year) k=floor(k)+koff #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 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) if koff==0: #solar itsa=eclipses_solar_type(GAMMA,u) #52=54.2 [maybe limit to <=1?] magecl_solar=(1.5433+u-abs(GAMMA))/ (0.5461+2*u) #only when partial #print(itsa,magecl_solar,GAMMA) if not(koff==0): #lunar -0.5 or 0.5 #magnitude of eclipse (<0: no eclipse) #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) 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 itsa=itsa+" lunar eclipse" if True: #print print(itsa) 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 [itsa,DTm(JDE),magecl_umbral,GAMMA] return [itsa,DTm(JDE),magecl_solar,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 print(crid) #print("See aeclipseuse for usage") #print("OK") #Is source complete? NO! #we are at the Numworks limit #ok