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")