In [ ]:

    

import batoid
import os
import yaml
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
import ipywidgets as widgets
%matplotlib inline


    

In [ ]:

    

D = 0.1
m0 = batoid.ConstMedium(1.000277)  # air
w, n = np.genfromtxt(os.path.join(batoid.datadir, "media", "silica_dispersion.txt")).T
w *= 1e-6
m1 = batoid.TableMedium(batoid.Table(w, n, 'linear'))
wavelength = 500e-9
n0 = m0.getN(wavelength)
n1 = m1.getN(wavelength)
R1 = 0.5  # convex entrance surface
R2 = -0.5  # convex exit surface
d = 0.01
# Lens maker equation:
f = 1./((n1-n0)*(1./R1 - 1./R2 + (n1-n0)*d/R1/R2))


    

In [ ]:

    

def makeSystem(focus):
    systemStr = """
    opticalSystem:
      type: CompoundOptic
      inMedium: 1.000277
      items:
        -
          type: RefractiveInterface
          surface:
            type: Sphere
            R: {R1}
          obscuration:
            type: ClearCircle
            radius: {Do2}
          inMedium: 1.000277
          outMedium: &silica
            type: SellmeierMedium
            B1: 0.6961663
            B2: 0.4079426
            B3: 0.8974794
            C1: 0.00467914825849
            C2: 0.013512063073959999
            C3: 97.93400253792099
          name: L1
        -
          type: RefractiveInterface
          surface:
            type: Sphere
            R: {R2}
          obscuration:
            type: ClearCircle
            radius: {Do2}
          coordSys:
            z: {d}
          inMedium: *silica
          outMedium: 1.000277
          name: L2
        -
          type: Detector
          surface:
            type: Plane
          obscuration:
            type: ClearCircle
            radius: {Do2}
          coordSys:
            z: {f}
          inMedium: 1.000277          
          name: D
    """.format(**dict(R1=R1, R2=R2, Do2=D/2, d=d, f=f+focus))
    config = yaml.safe_load(systemStr)
    system = batoid.parse.parse_optic(config['opticalSystem'])
    return system


    

In [ ]:

    

@interact(theta_x=widgets.FloatSlider(min=-5,max=5,step=1,value=1.1, continuous_update=False),
          theta_y=widgets.FloatSlider(min=-5,max=5,step=1,value=2.2, continuous_update=False),
          wavelength=widgets.FloatSlider(min=0.4,max=0.7,step=0.05,value=0.5, continuous_update=False),
          focus=widgets.FloatSlider(min=-20, max=20, step=1,value=0.0, continuous_update=False))
def spot(theta_x, theta_y, wavelength, focus):
    """Display a spot diagram for a Newtonian telescope.

    @param theta_x  Field angle in degrees
    @param theta_y  Field angle in degrees
    @param wavelength in microns
    @param focus    Defocus distance in mm
    """
    system = makeSystem(focus*1e-3)
    dirCos = batoid.utils.fieldToDirCos(np.deg2rad(theta_x), np.deg2rad(theta_y))
    # Flip the z-cosine since we constructed our telescope to point down
    # instead of up, and hence need the initial rays to point up instead
    # of down.
    dirCos = dirCos[0:2]+(-dirCos[2],)
    rays = batoid.RayVector.asPolar(
        backDist=1.0, nrad=16, naz=64, wavelength=wavelength*1e-6,
        dirCos=dirCos,
        outer=D/2*0.999, medium=batoid.ConstMedium(1.000277)
    )
    
    system.trace(rays)
    rays.trimVignetted()
    plt.scatter(1e3*(rays.x-np.mean(rays.x)), 1e3*(rays.y-np.mean(rays.y)), s=1)
    plt.xlim(-1, 1)
    plt.ylim(-1, 1)
    plt.xlabel("mm")
    plt.ylabel("mm")
    plt.title(r"$\theta_x = {:4.2f}\,,\theta_y = {:4.2f}\,,\lambda = {:4g}\,, f={:4.2f}$".format(theta_x, theta_y, wavelength, focus))
    plt.axes().set_aspect('equal')


    

In [ ]: