In [ ]:

    

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


    

In [ ]:

    

HSC_telescope = batoid.Optic.fromYaml("HSC.yaml")
LSST_telescope = batoid.Optic.fromYaml("LSST_r.yaml")
DECam_telescope = batoid.Optic.fromYaml("DECam.yaml")


    

In [ ]:

    

# HSC
max_th = 0.74
pixSize = 15e-6
telescope = HSC_telescope

dthx_drx = []
dthy_dry = []
thetas = np.linspace(0., max_th, 20)
for th in thetas:
    dth_dr = batoid.psf.dthdr(telescope, np.deg2rad(th), 0, 620e-9)
    dthx_drx.append(abs(dth_dr[0,0]))
    dthy_dry.append(abs(dth_dr[1,1]))

dthx_drx = np.array(dthx_drx) * 206265 * pixSize # arcsec per pixel
dthy_dry = np.array(dthy_dry) * 206265 * pixSize # arcsec per pixel

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 4))
ax1.plot(thetas, dthx_drx, c='b', label='radial')
ax1.plot(thetas, dthy_dry, c='r', label='tangential')
ax1.set_xlabel('radius (deg)')
ax1.set_ylabel('pixel size (arcsec)')
ax1.legend(loc='lower left')
ax1.set_title("HSC")

ax2.plot(thetas, np.array(dthx_drx)/np.array(dthy_dry))
ax2.set_xlabel('radius (deg)')
ax2.set_ylabel('b/a')
ax2.set_title("HSC")

fig.tight_layout()
fig.show()


    

In [ ]:

    

# DECam
max_th = 1.1
pixSize = 15e-6
telescope = DECam_telescope

dthx_drx = []
dthy_dry = []
thetas = np.linspace(0., max_th, 20)
for th in thetas:
    dth_dr = batoid.psf.dthdr(telescope, np.deg2rad(th), 0, 620e-9)
    dthx_drx.append(abs(dth_dr[0,0]))
    dthy_dry.append(abs(dth_dr[1,1]))

dthx_drx = np.array(dthx_drx) * 206265 * pixSize # arcsec per pixel
dthy_dry = np.array(dthy_dry) * 206265 * pixSize # arcsec per pixel

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 4))
ax1.plot(thetas, dthx_drx, c='b', label='radial')
ax1.plot(thetas, dthy_dry, c='r', label='tangential')
ax1.set_xlabel('radius (deg)')
ax1.set_ylabel('pixel size (arcsec)')
ax1.legend(loc='lower left')
ax1.set_title("DECam")

ax2.plot(thetas, np.array(dthx_drx)/np.array(dthy_dry))
ax2.set_xlabel('radius (deg)')
ax2.set_ylabel('b/a')
ax2.set_title("DECam")

fig.tight_layout()
fig.show()


    

In [ ]:

    

# LSST
max_th = 1.75
pixSize = 10e-6
telescope = LSST_telescope

dthx_drx = []
dthy_dry = []
thetas = np.linspace(0., max_th, 20)
for th in thetas:
    dth_dr = batoid.psf.dthdr(telescope, np.deg2rad(th), 0, 620e-9)
    dthx_drx.append(abs(dth_dr[0,0]))
    dthy_dry.append(abs(dth_dr[1,1]))

dthx_drx = np.array(dthx_drx) * 206265 * pixSize # arcsec per pixel
dthy_dry = np.array(dthy_dry) * 206265 * pixSize # arcsec per pixel

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 4))
ax1.plot(thetas, dthx_drx, c='b', label='radial')
ax1.plot(thetas, dthy_dry, c='r', label='tangential')
ax1.set_xlabel('radius (deg)')
ax1.set_ylabel('pixel size (arcsec)')
ax1.legend(loc='lower left')
ax1.set_title("LSST")

ax2.plot(thetas, np.array(dthx_drx)/np.array(dthy_dry))
ax2.set_xlabel('radius (deg)')
ax2.set_ylabel('b/a')
ax2.set_title("LSST")

fig.tight_layout()
fig.show()


    

In [ ]:

    

# All on the same plot!
# And show residuals to third order fits
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6))

ax1.set_xlabel('radius (deg)', fontsize=18)
ax1.set_ylabel('relative pixel size on sky', fontsize=18)

ax2.set_xlabel('radius (deg)', fontsize=18)
ax2.set_ylabel('radial/tangential', fontsize=18)

fig2, (ax3, ax4) = plt.subplots(ncols=2, figsize=(12, 6))

for max_th, pixSize, telescope, color, name in zip(
    [0.75, 1.1, 1.75], 
    [15e-6, 15e-6, 10e-6], 
    [HSC_telescope, DECam_telescope, LSST_telescope],
    ['red', 'blue', 'green'],
    ['HSC', 'DECam', 'LSST'],
):

    dthx_drx = []
    dthy_dry = []
    thetas = np.linspace(0., max_th, 20)
    r = []
    for th in thetas:
        dth_dr = batoid.psf.dthdr(telescope, np.deg2rad(th), 0, 620e-9)
        dthx_drx.append(abs(dth_dr[0,0]))
        dthy_dry.append(abs(dth_dr[1,1]))
        chiefRay = batoid.Ray.fromStop(0.0, 0.0, optic=telescope, wavelength=620e-9, theta_x=np.deg2rad(th), theta_y=0)
        telescope.trace(chiefRay)
        r.append(chiefRay.x)
    r = np.array(r)
    dthx_drx = np.array(dthx_drx) * 206265 * pixSize # arcsec per pixel
    dthy_dry = np.array(dthy_dry) * 206265 * pixSize # arcsec per pixel

    ax1.plot(thetas, dthx_drx/dthx_drx[0], c=color, ls=':', label="{} radial".format(name))
    ax1.plot(thetas, dthy_dry/dthy_dry[0], c=color, label="{} tangential".format(name))
    ax2.plot(thetas, np.array(dthx_drx)/np.array(dthy_dry), c=color, label=name)

    print(telescope.name)
    print("th_r: m->rad  ", np.polyfit(r, np.deg2rad(thetas), 3))
    print("r_th: rad->m  ", np.polyfit(np.deg2rad(thetas), r, 3))
    print("dthrdr_r: m->arcsec/pix ", np.polyfit(r, dthx_drx, 3))
    print("dthtdt_r: m->arcsec/pix", np.polyfit(r, dthy_dry, 3))

    ax3.plot(
        thetas[1:], 
        r[1:]/np.poly1d(np.polyfit(np.deg2rad(thetas), r, 3))(np.deg2rad(thetas[1:])), 
        c=color,
        label="{} r_th".format(name)
    )
    ax4.plot(
        r[1:],
        np.deg2rad(thetas[1:])/np.poly1d(np.polyfit(r, np.deg2rad(thetas), 3))(r[1:]),
        ls='--',
        c=color,
        label="{} th_r".format(name)
    )
    ax4.plot(
        r[1:],
        dthx_drx[1:]/np.poly1d(np.polyfit(r, dthx_drx, 3))(r[1:]), 
        c=color,
        ls=':', label="{} dthrdr_r".format(name)
    )
    ax4.plot(
        r[1:],
        dthy_dry[1:]/np.poly1d(np.polyfit(r, dthy_dry, 3))(r[1:]), 
        c=color,
        ls='-', label="{} dthtdt_r".format(name)
    )

    
for ax in (ax1, ax2):
    ax.tick_params(axis='both', which='major', labelsize=18)
    ax.tick_params(axis='both', which='minor', labelsize=18)

ax1.legend(fontsize=14)
ax2.legend(fontsize=14)

fig.tight_layout()
fig.show()

ax3.legend(fontsize=14)
ax4.legend(fontsize=14)

fig2.tight_layout()
fig2.show()


    

In [ ]: