Skip to content

Instantly share code, notes, and snippets.

@erikson1970
Created September 18, 2017 20:03
Show Gist options
  • Save erikson1970/915c20abbc1b1b282a35bf7d92a49b12 to your computer and use it in GitHub Desktop.
Save erikson1970/915c20abbc1b1b282a35bf7d92a49b12 to your computer and use it in GitHub Desktop.
Find solar eclipses in the future by looking for matching the Declination and Right-Ascension of the sun and moon
import ephem
import math
def FindEclipes():
DATE='2017/8/21 '
pi=3.14159276
ME = ephem.Moon('2017/8/21 16:48:32')
SE = ephem.Sun('2017/8/21')
InAPass=False
oo=ephem.Observer()
oo.lat='42.293010'
oo.lon='-71.326204'
DDD=[0,31,28,31,30,31,30,31,31,30,31,30,31]
MIN_DIST=1.
MIN_DATE=-1
MIN_TIME=-1
AllMEs=[]
AllSEs=[]
AllDISTs=[]
BigIncrement=31
LittleIncrement=11
#LAST_DIST
for YY in range(2017,2022,1):
for MM in range(1,13):
DM=1+DDD[MM]
if YY % 4 == 0:
DM +=1
for DD in range(1,DM):
DATE='%d/%d/%d'%(YY,MM,DD)
HM=0
Increment=BigIncrement
while HM<1440:
ME.compute('%s %02d:%02d:00' %(DATE,HM / 60,HM % 60))
SE.compute('%s %02d:%02d:00' %(DATE,HM / 60,HM % 60))
DIST=math.sqrt((ME.ra.real-SE.ra.real)**2.+(ME.dec.real-SE.dec.real)**2.)*180./pi
if DIST<1.:
if InAPass:
MEs.append((ME.ra.real,ME.dec.real))
SEs.append((SE.ra.real,SE.dec.real))
DISTs.append((YY*10**6+MM*10**4+DD*100 + HM/100.,DIST))
else:
#Pass Starting
InAPass=True
Increment=LittleIncrement
MEs=[(ME.ra.real,ME.dec.real)]
SEs=[(SE.ra.real,SE.dec.real)]
DISTs=[(YY*10**6+MM*10**4+DD*100 + HM/100.,DIST)]
if HM%5<1:
print "%s %02d:%02d:00 %03.3f | MOON: RA %03.3f / el %03.3f"%(DATE,HM / 60,HM % 60,DIST,ME.ra.real,ME.dec.real)
else:
if InAPass:
#Pass Complete
InAPass=False
Increment=BigIncrement
AllMEs.append(MEs)
AllSEs.append(SEs)
AllDISTs.append(DISTs)
HM+=Increment
return AllMEs,AllSEs,AllDISTs
if __name__ == "__main__":
%matplotlib inline
import matplotlib.pyplot as plt
AllMEs,AllSEs,AllDISTs=FindEclipes()
fig, ax = plt.subplots(len(AllMEs), 2, figsize=(20,5*len(AllMEs)))
for i in range(len(AllMEs)):
xx,yy=zip(*AllMEs[i])
ax[i,0].plot(xx,yy,'rx',label='Moon')
xx,yy=zip(*AllSEs[i])
ax[i,0].plot(xx,yy,'bo',label='Sun')
xx,yy=zip(*AllDISTs[i])
ax[i,1].plot(xx,yy,'x')
ax[i,0].set_title("RA vs DEC")
ax[i,1].set_title("Distance (deg) vs Time")
#print AllMEs
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment