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")
    # add colorbar.
    cbar = m.colorbar(cs,location='bottom',pad="5%")
    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)