Astronomy: Eclipse V0.2
Algorithms based on Meeus. See aeclipse.py for general description.
V0.2 fixes a bug in my implementation of solar eclipse time, and returns results. Note: It reports a magnitude of 2.0 for the solar eclipse of October 2023, whereas online sources suggest 0.9…Initial investigation suggests that the calculation used should only apply to a partial eclipse, and a different calculation should be used for annular eclipses [or NaN until new calculation included]. I am surprised that an annular eclipse is apparently not considered a type of partial eclipse. Workaround: Until fixed, ignore magnitude when the solar eclipse is not partial. Magnitudes of partial solar eclipses are ok, but not annular, hybrid or total. (The magnitude of lunar eclipses is OK). Valid solar magnitudes are accurate to 2 sig figs.
Note: To find the 2008 Feb 7 annular solar eclipse , an unusually precise date is required: Eg 2008+1.5/12 (1/12 or 2/12 will miss it). [So a smaller calendar stepsize is needed].
"""Astronomy: Eclipse """ __version__="0.2" crid="Astronomy:Eclipse V"+__version__+""" © 2023 SteveG1CMZ""" from math import * #customise showe=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. """ 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 """ #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( koff,M,MP,F1,A1,OMEGA,E45,True) JDE=JDE+adjustdays #JDE is now the peak time of a potential eclipse angled=[round(M,3),round(MP,3),round(F,3),round(OMEGA,3), round(F1,3)] #degrees 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 magecl_umbral>0 or magecl_penumbral>0: #semidurations semis = [round(totalumbralsdur),round(partialumbralsdur),round(penumbralsdur),"(semi)"] else: semis=[] #when?can this happen? return [itsa,DTm(JDE),round(magecl_umbral,3), round(GAMMA,4),angled,JDE,semis] #lunar return [itsa,DTm(JDE),round(magecl_solar,3), round(GAMMA,4),angled,JDE] #solar #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 if showe or __name__=="__main__": print(crid) #print("See aeclipseuse for usage") if True: #examples print(eclipse(2023+10/12,0)) #solar print(eclipse(2023+10/12,0.5)) #lunar #print("OK") #Is source complete? NO! #we are at the Numworks limit #ok #eclipse results: None, or #itsa: type of eclipse or None #JDE: peak time #magecl measures fraction eclipsed #GAMMA measures how central the eclipse is #angled is angles in degrees #(M,MP,F,OMEGA,F1) #lunar only: mins are semidurations