In [ ]:
import batoid
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
import ipywidgets as widgets
import galsim
%matplotlib inline

In [ ]:
def makeTelescope(L, F, dz):
    """
    Parameters
    ----------
    L : float
        focal length in meters
    F : float
        F-number
    dz : float
        Focal plane offset        
    """
    R = 2*L # radius of curvature    
    telescope = batoid.CompoundOptic(
        items = [
            batoid.Mirror(
                batoid.Paraboloid(R),
                name="Mirror"
            ),
            batoid.Detector(
                batoid.Plane(),
                name="Detector",
                coordSys=batoid.CoordSys().shiftGlobal([0,0,L+dz])
            )
        ]
    )
    telescope.backDist = 1.1*L
    telescope.pupilSize = L/F
    telescope.sphereRadius = L
    telescope.stopSurface = batoid.Interface(batoid.Plane())
    return telescope

In [ ]:
def a4coef(L, F, dz, wavelength):
    R = L
    alpha = dz/R
    term1 = alpha*R
    term1 /= 4*np.sqrt(3)*(1-alpha)
    term1 /= (2*F)**2
    term2 = alpha*(1+alpha+alpha**2)*R
    term2 /= 16*np.sqrt(3)*(1-alpha)**3
    term2 /= (2*F)**4
    return (term1+term2)/wavelength

def a11coef(L, F, dz, wavelength):
    R = L
    alpha = dz/R
    result = alpha*(1+alpha*alpha**2)*R
    result /= 48*np.sqrt(5)*(1-alpha)**3
    result /= (2*F)**4
    return -result/wavelength

In [ ]:
@interact(
    L = widgets.FloatSlider(min=1.0, max=20.0, step=0.1, value=10.0,
                            description="L (m)"),
    F = widgets.FloatSlider(min=1.0, max=10.0, step=0.05, value=1.25,
                            description="F/#"),
    dz = widgets.FloatSlider(min=-1000, max=1000, step=10, value=0.0,
                             description="dz ($\\mu m$)"),
    theta_x = widgets.FloatSlider(min=-1.75, max=1.75, step=0.05, value=0.0,
                                  description="$\\theta_x (deg)$"),
    theta_y = widgets.FloatSlider(min=-1.75, max=1.75, step=0.05, value=0.0,
                                  description="$\\theta_y (deg)$")
)
def zernike(L, F, dz, theta_x, theta_y):
    telescope = makeTelescope(L, F, dz*1e-6)
    wavelength = 750e-9
    z = batoid.analysis.zernike(
        telescope, np.deg2rad(theta_x), np.deg2rad(theta_y), 
        wavelength, jmax=22, nx=128
    )
    for i in range(1, len(z)//2+1):
        print("{:6d}   {:6.3f}      {:6d}  {:6.3f}".format(i, z[i], i+11, z[i+11]))    
        
    print("a4 analytic: {:6.3f}".format(a4coef(L, F, dz*1e-6, wavelength)))
    print("a11 analytic: {:6.3f}".format(a11coef(L, F, dz*1e-6, wavelength)))

In [ ]: