mport numpy as np
from scipy.special import ellipeinc
import matplotlib.pyplot as plt

f = 1/298.257223563
esq = f*(2-f)
e = np.sqrt(esq)

def gd(x):
    return 2*np.arctan(np.tanh(x/2))

def gdinv(x):
    return 2*np.arctanh(np.tan(x/2))

npts = 181
geodetic = np.linspace(0,np.pi/2,npts)
geocentric = np.arctan((1-f)**2 * np.tan(geodetic))
parametric = np.arctan((1-f) * np.tan(geodetic))
m = (ellipeinc(geodetic, esq) -
     esq * np.sin(geodetic) * np.cos(geodetic) / 
     np.sqrt(1 - esq * np.sin(geodetic)**2))
rectifying = np.pi/2*m/m.max()
qp = 1 + (1-esq)/e*np.arctanh(e)
q = ((1 - esq)*np.sin(geodetic)/(1 - esq*np.sin(geodetic)**2) + 
     (1 - esq)/e*np.arctanh(e*np.sin(geodetic)))
authalic = np.arcsin(q/qp)
conformal = gd(gdinv(geodetic) - e * np.arctanh(e*np.sin(geodetic)))

gddegrees = np.linspace(0,90,npts)
gdmindiff = np.zeros(gddegrees.shape)
parmindiff = 10800/np.pi*(parametric - geodetic)
autmindiff = 10800/np.pi*(authalic - geodetic)
recmindiff = 10800/np.pi*(rectifying - geodetic)
conmindiff = 10800/np.pi*(conformal - geodetic)
gcmindiff = 10800/np.pi*(geocentric - geodetic)

fig = plt.figure()
ax = plt.axes()
#in this order so conformal and geocentric get contrasting colors
ax.plot(gddegrees, conmindiff, label='conformal', dashes=[4, 4])
ax.plot(gddegrees, gcmindiff, label='geocentric', dashes=[0, 4, 4, 0])
ax.plot(gddegrees, gdmindiff, label='geodetic')
ax.plot(gddegrees, parmindiff, label='parametric')
ax.plot(gddegrees, autmindiff, label='authalic')
ax.plot(gddegrees, recmindiff, label='rectifying')
bbox = dict(color='white', alpha=0)
ax.text(45, -0.1, 'Geodetic', bbox=bbox,
        horizontalalignment='center', verticalalignment='top')
ax.text(45, min(parmindiff) + 0.2, 'Parametric', bbox=bbox, 
        horizontalalignment='center', verticalalignment='bottom')
ax.text(45, min(autmindiff) + 0.15, 'Authalic', bbox=bbox, 
        horizontalalignment='center', verticalalignment='bottom')
ax.text(45, min(recmindiff) + 0.25, 'Rectifying', bbox=bbox, 
        horizontalalignment='center', verticalalignment='bottom')
ax.text(45, min(gcmindiff) + 0.35, 'Conformal/\nGeocentric', bbox=bbox,
        horizontalalignment='center', verticalalignment='bottom')
ax.set(xlim=(0, 90),
       xlabel='Geodetic latitude (degrees)', 
       ylabel='Difference from geodetic latitude (minutes)',
       title='Auxiliary latitudes');
ax.xaxis.set_major_locator(plt.MaxNLocator(7))
plt.savefig("Auxiliary Latitudes Difference.svg")