### Eclipse visualizations using Python

I have been fine tuning the scripts for calculating the eclipse fraction, so I can calculate the extreme ultraviolet radiation variation during the solar eclipse. The script can also be used for nice visualizations, like this:
 This is the fraction of sunlight during this year's highly publicized US eclipse, which is due on August 21st.

 The Svalbard eclipse of 2015.

Here's the code, for those interested. I use pyephem to calculate the location of the Moon and the Sun.
``````
#!/usr/bin/env python
#
# Calculate the totality of the eclipse.
#
import ephem
import numpy as n
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
``````
``````# numerically approximate eclipse fraction
def intersection(r0,r1,d,n_s=100):
A1=n.zeros([n_s,n_s])
A2=n.zeros([n_s,n_s])
I=n.zeros([n_s,n_s])
x=n.linspace(-2.0*r0,2.0*r0,num=n_s)
y=n.linspace(-2.0*r0,2.0*r0,num=n_s)
xx,yy=n.meshgrid(x,y)
A1[n.sqrt((xx+d)**2.0+yy**2.0) < r0]=1.0
n_sun=n.sum(A1)
A2[n.sqrt(xx**2.0+yy**2.0) < r1]=1.0
S=A1+A2
I[S>1]=1.0
eclipse=n.sum(I)/n_sun
return(eclipse)
```
```
Calculate the eclipsed totality for a grid of points:
``````
# calculate the fraction that sun is eclipsed at given altitudes, latitude, and longitude
#
# returns eclipse fraction (0..1) and time (seconds after t0)
def get_eclipse(t0,n_t=60*10,dt=60.0,alts=n.linspace(0,600e3,num=600),lats=[69.0],lons=[16.02]):
# Location
obs = ephem.Observer()
n_alts=len(alts)
n_lats=len(lats)
n_lons=len(lons)

p=n.zeros([n_t,n_alts,n_lats,n_lons])
times=n.arange(n_t)*dt
dts=[]
for ti,t in enumerate(times):
print("Time %1.2f (s)"%(t))
for ai,alt in enumerate(alts):
for lai,lat in enumerate(lats):
for loi,lon in enumerate(lons):
#obs.lon, obs.lat = '-1.268304', '51.753101'#'16.02', '78.15' # ESR
obs.lon, obs.lat = '%1.2f'%(lon), '%1.2f'%(lat) # ESR
obs.elevation=alt
obs.date= t0#(ephem.date(ephem.date(t0)+t*ephem.second))
sun, moon = ephem.Sun(), ephem.Moon()

# Output list
results=[]
seps=[]
sun.compute(obs)
moon.compute(obs)
r_sun=(sun.size/2.0)/3600.0
r_moon=(moon.size/2.0)/3600.0
s=n.degrees(ephem.separation((sun.az, sun.alt), (moon.az, moon.alt)))
percent_eclipse=0.0

if s < (r_moon+r_sun):
#                        print("eclipsed")
if s < 1e-3:
percent_eclipse=1.0
else:
percent_eclipse=intersection(r_moon,r_sun,s,n_s=100)

if n.degrees(sun.alt) <= r_sun:
if n.degrees(sun.alt) <= -r_sun:
percent_eclipse=1.0
else:
percent_eclipse=1.0-((n.degrees(sun.alt)+r_sun)/(2.0*r_sun))*(1.0-percent_eclipse)

p[ti,ai,lai,loi]=percent_eclipse
dts.append(obs.date)
return(p,times,dts)
```
```
Then a small function to plot eclipse fraction as a colormesh using basemap and matplotlib:
``````
def plot_map(p,lons,lats,t0,alt=0,lat_0=69,lon_0=16):
fig = plt.figure(figsize=(8,8))
plt.style.use('dark_background')
# create polar stereographic Basemap instance.
m = Basemap(projection='gnom',lon_0=lon_0,lat_0=lat_0, resolution='l',width=15.e6,height=15.e6)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
m.drawcountries()

lon,lat=n.meshgrid(lons,lats)
x, y = m(lon, lat) # compute map proj coordinates.
# draw filled contours.
cs = m.pcolormesh(x,y,1.0-p[0,0,:,:],vmin=0,vmax=1.0,cmap="inferno")
plt.title("Fraction of sunlight at %1.2f km\n%s (UTC)"%(alt/1e3,t0))
fname="eclipse-%f.png"%(float(t0))
plt.text(1e5,1e5,"University of Tromso",color="black")
print("saving %s"%(fname))
plt.savefig(fname)
```
```
Create images with 900 second time step:
``````
for i in range(20):
create_map(t0=ephem.date((2017,8,21,16,15,0))+ephem.second*900*i,lat_0=38.496077,lon_0=360-90.152751,nlats=300,nlons=600)

```
```