Skip to content

Instantly share code, notes, and snippets.

@jmeyers314

jmeyers314/Batoid ImSim integration.ipynb

Last active October 30, 2019 17:32
Show Gist options
  • Save jmeyers314/3c5e69ada9d2276fac401eefa8b42f9f to your computer and use it in GitHub Desktop.
Save jmeyers314/3c5e69ada9d2276fac401eefa8b42f9f to your computer and use it in GitHub Desktop.
Download ZIP
Batoid ImSim integration
Display the source blob
Display the rendered blob
Raw
Batoid ImSim integration.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
BatoidWCS.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-11T20:25:58.314624Z",
"start_time": "2019-10-11T20:25:53.884963Z"
}
},
"outputs": [],
"source": [
"import lsst.afw.cameraGeom as cameraGeom\n",
"from lsst.obs.lsst import LsstCamMapper\n",
"import lsst.obs.lsst.cameraTransforms as cameraTransforms\n",
"\n",
"import batoid\n",
"import galsim\n",
"\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-11T20:26:04.503467Z",
"start_time": "2019-10-11T20:25:58.316967Z"
}
},
"outputs": [],
"source": [
"camera = LsstCamMapper().camera"
]
},
{
"cell_type": "code",
"execution_count": 1028,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-11T20:26:04.520761Z",
"start_time": "2019-10-11T20:26:04.508135Z"
}
},
"outputs": [],
"source": [
"class BatoidWCS(galsim.wcs.CelestialWCS):\n",
" \"\"\"\n",
" Parameters\n",
" ----------\n",
" boresight : galsim.CelestialCoord\n",
" rotSkyPos : galsim.Angle\n",
" rotTelPos : galsim.Angle\n",
" fiducial_telescope : batoid.Optic\n",
" Telescope before application of rotator\n",
" wavelength : float\n",
" Nanometers\n",
" camera : lsst.afw.cameraGeom.camera.camera.Camera\n",
" det : lsst.afw.cameraGeom.detector.detector.Detector\n",
" \"\"\"\n",
" def __init__(self, boresight, rotSkyPos, rotTelPos, fiducial_telescope, wavelength, camera, det):\n",
" self.boresight = boresight\n",
" self.rotSkyPos = rotSkyPos\n",
" self.rotTelPos = rotTelPos\n",
" self.fiducial_telescope = fiducial_telescope\n",
" self.wavelength = wavelength\n",
" self.camera = camera\n",
" self.det = det\n",
" \n",
" # Initialize RaDec -> field\n",
" q = rotTelPos - rotSkyPos\n",
" cq, sq = np.cos(q), np.sin(q) # +q or -q here?\n",
" jac = [cq, -sq, sq, cq]\n",
" affine = galsim.AffineTransform(*jac)\n",
" self._radec2field = galsim.TanWCS(affine, boresight, units=galsim.radians)\n",
"\n",
" # Initialize field -> focal\n",
" _pz = fiducial_telescope.stopSurface.surface.sag(0.0, 0.0)\n",
" transform = batoid.CoordTransform(fiducial_telescope.stopSurface.coordSys, batoid.globalCoordSys)\n",
" _, _, self._pz = transform.applyForward(0, 0, _pz)\n",
" self._refractive_index = fiducial_telescope.inMedium.getN(wavelength*1e-9)\n",
" self._telescope = fiducial_telescope.withLocallyRotatedOptic(\"LSSTCamera\", batoid.RotZ(-self.rotTelPos))\n",
"\n",
" # Initialize focal -> pixel\n",
" lct = cameraTransforms.LsstCameraTransforms(camera)\n",
" center = 509*4 + 0.5, 2000.5\n",
" # Focal plane coords of center of CCD\n",
" _x, _y = lct.ccdPixelToFocalMm(*center, det.getName())\n",
" self._focal_x = _x*1e-3 # meters\n",
" self._focal_y = _y*1e-3 # meters\n",
" \n",
" # Also calculate field angle of CCD center\n",
" self._field_x, self._field_y = self._focalToField(self._focal_x, self._focal_y)\n",
" \n",
" # Required by galsim.BaseWCS\n",
" self._color = None\n",
" self.origin = galsim.PositionD(0, 0)\n",
"\n",
" def _fieldToFocal(self, field_x, field_y):\n",
" ray = batoid.Ray.fromStop(\n",
" 0.0, 0.0,\n",
" optic=self._telescope,\n",
" wavelength=self.wavelength*1e-9,\n",
" theta_x=field_x, theta_y=field_y, projection='gnomonic' \n",
" )\n",
" self._telescope.traceInPlace(ray)\n",
" return ray.x, ray.y\n",
" \n",
" def _fieldToFocalMany(self, field_x, field_y):\n",
" # Setup rays to trace with batoid. For each input coordinate, we'll \n",
" # trace exactly 1 ray, the chief ray. Any displacements between the \n",
" # mean optics/atmosphere and the chief ray can then be considered \n",
" # part of the PSF.\n",
" n = len(field_x)\n",
" dirCos = np.array(batoid.utils.fieldToDirCos(field_x, field_y, projection='gnomonic'))\n",
" w = np.empty(n); w.fill(self.wavelength*1e-9)\n",
" pz = np.empty(n); pz.fill(self._pz)\n",
" rays = batoid.RayVector.fromArrays(\n",
" np.zeros(n), np.zeros(n), pz, # Chief ray -> start at center of the entrance pupil.\n",
" *(dirCos/self._refractive_index), # velocity = dirCos adjusted for refractive index\n",
" np.zeros(n), # doesn't matter what we set time to, so just use 0.0\n",
" w # wavelength\n",
" )\n",
" self._telescope.traceInPlace(rays)\n",
" return rays.x, rays.y\n",
" \n",
" def _focalToField(self, x, y):\n",
" from scipy.optimize import least_squares\n",
"\n",
" def _resid(args):\n",
" fx1, fy1 = args\n",
" x1, y1 = self._fieldToFocal(fx1, fy1)\n",
" return np.array([x1-x, y1-y])\n",
"\n",
" result = least_squares(_resid, np.array([0.0, 0.0]))\n",
" return result.x[0], result.x[1]\n",
"\n",
" def _focalToFieldMany(self, x, y):\n",
" if not hasattr(self, '_xinterp'):\n",
" from scipy.interpolate import RectBivariateSpline\n",
" # radius of circle circumscribing CCD, plus some fudge\n",
" half_length = 2000*10e-6*np.sqrt(2)*1.2\n",
" N = 4 # seems to be sufficient, kinda surprising...\n",
" focal_xs = self._focal_x + half_length*np.linspace(-1,1,N)\n",
" focal_ys = self._focal_y + half_length*np.linspace(-1,1,N)\n",
" field_xs = np.empty((N, N), dtype=float)\n",
" field_ys = np.empty((N, N), dtype=float)\n",
" for i in range(N):\n",
" for j in range(N):\n",
" field_xs[i,j], field_ys[i,j] = self._focalToField(focal_xs[i], focal_ys[j])\n",
" self._xinterp = RectBivariateSpline(focal_xs, focal_ys, field_xs)\n",
" self._yinterp = RectBivariateSpline(focal_xs, focal_ys, field_ys)\n",
" return self._xinterp(x, y, grid=False), self._yinterp(x, y, grid=False)\n",
" \n",
" def _xy(self, ra, dec, color=None):\n",
" \"\"\"\n",
" Parameters\n",
" ----------\n",
" ra : array_like\n",
" Right ascension in radians\n",
" dec : array_like\n",
" Declination in radians\n",
" \n",
" Returns\n",
" -------\n",
" xpix, ypix : array_like\n",
" \"\"\"\n",
" # Abuse GalSim's TanWCS to convert ra/dec into field angles \n",
" # centered at ra/dec and aligned with alt/az.\n",
" field_x, field_y = self._radec2field.radecToxy(ra, dec, 'rad')\n",
"\n",
" # Field -> Focal via batoid.\n",
" focal_x, focal_y = self._fieldToFocalMany(field_x, field_y)\n",
" \n",
" # For the moment, focal -> pixel is just an offset and a scaling by the pixel size.\n",
" pixSize = 10e-6 # m\n",
" pixel_x = (focal_x - self._focal_x)/pixSize + (509*4 + 0.5)\n",
" pixel_y = (focal_y - self._focal_y)/pixSize + 2000.5\n",
" return pixel_x, pixel_y\n",
" \n",
" def _radec(self, x, y, color=None):\n",
" # reverse of _xy.\n",
"\n",
" # pixel -> focal\n",
" pixSize = 10e-6\n",
" focal_x = pixSize*(x-(509*4+0.5)) + self._focal_x\n",
" focal_y = pixSize*(y-2000.5) + self._focal_y\n",
"\n",
" # focal -> field using interpolating function derived from field->focal\n",
" field_x, field_y = self._focalToFieldMany(focal_x, focal_y)\n",
"\n",
" # field -> ra, dec using GalSim TanWCS\n",
" return self._radec2field.xyToradec(field_x, field_y, units='rad')"
]
},
{
"cell_type": "code",
"execution_count": 1038,
"metadata": {
"ExecuteTime": {
"end_time": "2019-10-11T20:26:04.610210Z",
"start_time": "2019-10-11T20:26:04.547881Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"R02_S21\n"
]
}
],
"source": [
"# Construct a WCS for some detector\n",
"det = camera[16]\n",
"print(det.getName())\n",
"\n",
"boresight = galsim.CelestialCoord(ra=30*galsim.degrees, dec=-10*galsim.degrees)\n",
"rotTelPos = 35.0*galsim.degrees\n",
"rotSkyPos = 10.0*galsim.degrees\n",
"\n",
"wcs = BatoidWCS(\n",
" boresight, rotSkyPos, rotTelPos,\n",
" batoid.Optic.fromYaml(\"LSST_r.yaml\"),\n",
" 620, # nm\n",
" camera, det\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 1039,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4.4306630586632416e-05 pixel\n"
]
}
],
"source": [
"# Test focal<->field roundtrip for individual points\n",
"# Mostly probing the accuracy of the scipy least_squares inversion. \n",
"errs = []\n",
"for _ in range(1000):\n",
" xy = np.random.uniform(-0.3, 0.3, size=2)\n",
" w = np.hypot(*xy)\n",
" if w > 0.3:\n",
" continue\n",
" errs.append((wcs._fieldToFocal(*wcs._focalToField(*xy)) - xy)/10e-6)\n",
"print(np.max(errs), \"pixel\")"
]
},
{
"cell_type": "code",
"execution_count": 1041,
"metadata": {},
"outputs": [],
"source": [
"# Roundtrip field -> focal -> field for many points at once.\n",
"# This uses an interpolation trained only in the region of the particular sensor, so \n",
"# don't wander too far away from that.\n",
"fieldx = wcs._field_x + np.deg2rad(2000*0.2/3600)*np.linspace(-1, 1, 800)\n",
"fieldy = wcs._field_y + np.deg2rad(2000*0.2/3600)*np.linspace(-1, 1, 800)\n",
"fieldx, fieldy = np.meshgrid(fieldx, fieldy)\n",
"fieldx = fieldx.ravel()\n",
"fieldy = fieldy.ravel()\n",
"\n",
"focalx, focaly = wcs._fieldToFocalMany(fieldx, fieldy)\n",
"fieldx1, fieldy1 = wcs._focalToFieldMany(focalx, focaly)"
]
},
{
"cell_type": "code",
"execution_count": 1042,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1.8315910536664856e-06 arcsec\n",
"-2.623459060261082e-06 arcsec\n",
"2.7497054732218163e-06 arcsec\n",
"3.9374238444648215e-06 arcsec\n"
]
}
],
"source": [
"print(3600*np.rad2deg(np.median(fieldx1-fieldx)), \"arcsec\")\n",
"print(3600*np.rad2deg(np.median(fieldy1-fieldy)), \"arcsec\")\n",
"\n",
"print(3600*np.rad2deg(np.median(np.abs(fieldx1-fieldx))), \"arcsec\")\n",
"print(3600*np.rad2deg(np.median(np.abs(fieldy1-fieldy))), \"arcsec\")"
]
},
{
"cell_type": "code",
"execution_count": 1043,
"metadata": {},
"outputs": [],
"source": [
"# Roundtrip focal -> field -> focal\n",
"focalx = wcs._focal_x + np.linspace(-2000*10e-6, 2000*10e-6, 80)\n",
"focaly = wcs._focal_y + np.linspace(-2000*10e-6, 2000*10e-6, 80)\n",
"focalx, focaly = np.meshgrid(focalx, focaly)\n",
"focalx = focalx.ravel()\n",
"focaly = focaly.ravel()\n",
"\n",
"fieldx, fieldy = wcs._focalToFieldMany(focalx, focaly)\n",
"focalx1, focaly1 = wcs._fieldToFocalMany(fieldx, fieldy)"
]
},
{
"cell_type": "code",
"execution_count": 1044,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-2.983737389106178e-08 pixel\n",
"-9.86531689672887e-06 pixel\n",
"2.320837749411608e-07 pixel\n",
"2.4014325250565346e-05 pixel\n"
]
}
],
"source": [
"print(np.median(focalx1-focalx)/10e-6, \"pixel\")\n",
"print(np.median(focaly1-focaly)/10e-6, \"pixel\")\n",
"\n",
"print(np.median(np.abs(focalx1-focalx))/10e-6, \"pixel\")\n",
"print(np.median(np.abs(focaly1-focaly))/10e-6, \"pixel\")"
]
},
{
"cell_type": "code",
"execution_count": 1046,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"R01_S00 7.2544181594384275e-06 arcsec\n",
"R01_S01 1.433698188485654e-05 arcsec\n",
"R01_S02 1.0243720362907117e-05 arcsec\n",
"R01_S10 1.1135579939677904e-05 arcsec\n",
"R01_S11 6.403266541792628e-06 arcsec\n",
"R01_S12 4.3049784043895906e-06 arcsec\n",
"R01_S20 5.1013095839763515e-06 arcsec\n",
"R01_S21 2.938842067336334e-06 arcsec\n",
"R01_S22 2.499024072583775e-06 arcsec\n",
"R02_S00 7.734938030505143e-06 arcsec\n",
"R02_S01 6.839891060860329e-06 arcsec\n",
"R02_S02 6.315694818984418e-06 arcsec\n",
"R02_S10 3.7584564577118143e-06 arcsec\n",
"R02_S11 3.874487684417425e-06 arcsec\n",
"R02_S12 4.002426728944234e-06 arcsec\n",
"R02_S20 2.515082693170484e-06 arcsec\n",
"R02_S21 2.7497054732218163e-06 arcsec\n",
"R02_S22 2.940118920891674e-06 arcsec\n",
"R03_S00 6.270739447745473e-06 arcsec\n",
"R03_S01 4.728063033819002e-06 arcsec\n",
"R03_S02 4.066500866608575e-06 arcsec\n",
"R03_S10 3.801634932564093e-06 arcsec\n",
"R03_S11 3.450567518303577e-06 arcsec\n",
"R03_S12 3.5868656098312167e-06 arcsec\n",
"R03_S20 2.6899331798830387e-06 arcsec\n",
"R03_S21 2.8450070620697365e-06 arcsec\n",
"R03_S22 3.27941683418751e-06 arcsec\n",
"R10_S00 7.98434114633824e-06 arcsec\n",
"R10_S01 1.1894259235862296e-05 arcsec\n",
"R10_S02 5.5671659450104484e-06 arcsec\n",
"R10_S10 1.2412012349821431e-05 arcsec\n",
"R10_S11 7.963518253623366e-06 arcsec\n",
"R10_S12 4.0037586856094305e-06 arcsec\n",
"R10_S20 1.347066975575221e-05 arcsec\n",
"R10_S21 6.1102042301456175e-06 arcsec\n",
"R10_S22 3.748152891793479e-06 arcsec\n",
"R11_S00 3.054649840197105e-06 arcsec\n",
"R11_S01 2.3180289629505335e-06 arcsec\n",
"R11_S02 2.0360890709030687e-06 arcsec\n",
"R11_S10 2.8091171811128166e-06 arcsec\n",
"R11_S11 2.169762953697558e-06 arcsec\n",
"R11_S12 1.7243557503247962e-06 arcsec\n",
"R11_S20 2.7216649676467637e-06 arcsec\n",
"R11_S21 2.017735262665882e-06 arcsec\n",
"R11_S22 1.4528247999479564e-06 arcsec\n",
"R12_S00 1.8810879588135593e-06 arcsec\n",
"R12_S01 2.0238459828615106e-06 arcsec\n",
"R12_S02 2.1910730079972984e-06 arcsec\n",
"R12_S10 1.4105771740586084e-06 arcsec\n",
"R12_S11 1.4753980214786736e-06 arcsec\n",
"R12_S12 1.526639798802465e-06 arcsec\n",
"R12_S20 1.0214354341290098e-06 arcsec\n",
"R12_S21 9.762393893233948e-07 arcsec\n",
"R12_S22 1.051031037128247e-06 arcsec\n",
"R13_S00 2.137109351734044e-06 arcsec\n",
"R13_S01 2.5018915811707152e-06 arcsec\n",
"R13_S02 3.101745197957306e-06 arcsec\n",
"R13_S10 1.757518139160202e-06 arcsec\n",
"R13_S11 2.2671283150251795e-06 arcsec\n",
"R13_S12 3.013809937749006e-06 arcsec\n",
"R13_S20 1.497208431695592e-06 arcsec\n",
"R13_S21 2.188673428578635e-06 arcsec\n",
"R13_S22 3.1698466483820334e-06 arcsec\n",
"R14_S00 4.2217800436607286e-06 arcsec\n",
"R14_S01 6.832533185538698e-06 arcsec\n",
"R14_S02 5.945353258263084e-06 arcsec\n",
"R14_S10 4.323315753256302e-06 arcsec\n",
"R14_S11 6.079756899450768e-06 arcsec\n",
"R14_S12 1.1456109720338699e-05 arcsec\n",
"R14_S20 4.307878831717484e-06 arcsec\n",
"R14_S21 5.943440393163785e-06 arcsec\n",
"R14_S22 1.0007420337222983e-05 arcsec\n",
"R20_S00 1.0828928986452875e-05 arcsec\n",
"R20_S01 5.506583653847402e-06 arcsec\n",
"R20_S02 3.7189901076206055e-06 arcsec\n",
"R20_S10 9.785994930888448e-06 arcsec\n",
"R20_S11 5.572951413732964e-06 arcsec\n",
"R20_S12 3.939686471186748e-06 arcsec\n",
"R20_S20 9.589423889598343e-06 arcsec\n",
"R20_S21 5.745308931990021e-06 arcsec\n",
"R20_S22 4.172441649407592e-06 arcsec\n",
"R21_S00 2.7216953817009056e-06 arcsec\n",
"R21_S01 1.9736783583675165e-06 arcsec\n",
"R21_S02 1.306465662436552e-06 arcsec\n",
"R21_S10 2.8951590769977083e-06 arcsec\n",
"R21_S11 2.1056587150649266e-06 arcsec\n",
"R21_S12 1.3914768796984484e-06 arcsec\n",
"R21_S20 3.146439992314739e-06 arcsec\n",
"R21_S21 2.324068836290633e-06 arcsec\n",
"R21_S22 1.5178732607820305e-06 arcsec\n",
"R22_S00 7.1918190586734e-07 arcsec\n",
"R22_S01 4.883078592584228e-07 arcsec\n",
"R22_S02 7.426067209140189e-07 arcsec\n",
"R22_S10 6.9465632569349e-07 arcsec\n",
"R22_S11 3.294776848423244e-08 arcsec\n",
"R22_S12 6.941710649870074e-07 arcsec\n",
"R22_S20 7.426626738283292e-07 arcsec\n",
"R22_S21 4.877863141380839e-07 arcsec\n",
"R22_S22 7.191050656540824e-07 arcsec\n",
"R23_S00 1.517322453316218e-06 arcsec\n",
"R23_S01 2.324434431112034e-06 arcsec\n",
"R23_S02 3.1471060601004372e-06 arcsec\n",
"R23_S10 1.3920246010322965e-06 arcsec\n",
"R23_S11 2.1065595077861174e-06 arcsec\n",
"R23_S12 2.895786501044027e-06 arcsec\n",
"R23_S20 1.3063172597429615e-06 arcsec\n",
"R23_S21 1.9739438551695517e-06 arcsec\n",
"R23_S22 2.721433105210485e-06 arcsec\n",
"R24_S00 4.175982560934465e-06 arcsec\n",
"R24_S01 5.747482642330133e-06 arcsec\n",
"R24_S02 9.603075148346651e-06 arcsec\n",
"R24_S10 3.941847479186602e-06 arcsec\n",
"R24_S11 5.573866340044021e-06 arcsec\n",
"R24_S12 9.81152269887281e-06 arcsec\n",
"R24_S20 3.7215638522257794e-06 arcsec\n",
"R24_S21 5.518217566275126e-06 arcsec\n",
"R24_S22 1.0854861082451119e-05 arcsec\n",
"R30_S00 1.0017236920462652e-05 arcsec\n",
"R30_S01 5.9422669473925256e-06 arcsec\n",
"R30_S02 4.308993954067273e-06 arcsec\n",
"R30_S10 1.1476235236682728e-05 arcsec\n",
"R30_S11 6.081246293572697e-06 arcsec\n",
"R30_S12 4.3230023095924425e-06 arcsec\n",
"R30_S20 5.945557569144434e-06 arcsec\n",
"R30_S21 6.830938773476881e-06 arcsec\n",
"R30_S22 4.225583768396624e-06 arcsec\n",
"R31_S00 3.170048812388974e-06 arcsec\n",
"R31_S01 2.188900907812993e-06 arcsec\n",
"R31_S02 1.4973043254192379e-06 arcsec\n",
"R31_S10 3.0144734113947926e-06 arcsec\n",
"R31_S11 2.2670820677722642e-06 arcsec\n",
"R31_S12 1.7577331173237774e-06 arcsec\n",
"R31_S20 3.101438016010477e-06 arcsec\n",
"R31_S21 2.501484905012948e-06 arcsec\n",
"R31_S22 2.1371761284735044e-06 arcsec\n",
"R32_S00 1.0514580415030835e-06 arcsec\n",
"R32_S01 9.759213835514153e-07 arcsec\n",
"R32_S02 1.022961861834511e-06 arcsec\n",
"R32_S10 1.5278763537360677e-06 arcsec\n",
"R32_S11 1.475563599167544e-06 arcsec\n",
"R32_S12 1.4116999893750321e-06 arcsec\n",
"R32_S20 2.1930427652684625e-06 arcsec\n",
"R32_S21 2.023010669809824e-06 arcsec\n",
"R32_S22 1.881710194580054e-06 arcsec\n",
"R33_S00 1.4539624644790492e-06 arcsec\n",
"R33_S01 2.019470831720156e-06 arcsec\n",
"R33_S02 2.724100596664895e-06 arcsec\n",
"R33_S10 1.724750596310033e-06 arcsec\n",
"R33_S11 2.1714554063574315e-06 arcsec\n",
"R33_S12 2.811880208478471e-06 arcsec\n",
"R33_S20 2.036897369118429e-06 arcsec\n",
"R33_S21 2.318973229878527e-06 arcsec\n",
"R33_S22 3.051145783347603e-06 arcsec\n",
"R34_S00 3.749451035186719e-06 arcsec\n",
"R34_S01 6.1130645824845245e-06 arcsec\n",
"R34_S02 1.3462606453280697e-05 arcsec\n",
"R34_S10 4.000617092722815e-06 arcsec\n",
"R34_S11 7.975109228562888e-06 arcsec\n",
"R34_S12 1.239076437597338e-05 arcsec\n",
"R34_S20 5.553458151902675e-06 arcsec\n",
"R34_S21 1.189067681809683e-05 arcsec\n",
"R34_S22 7.97652580786108e-06 arcsec\n",
"R41_S00 3.27887620201187e-06 arcsec\n",
"R41_S01 2.844243222045283e-06 arcsec\n",
"R41_S02 2.691016367475978e-06 arcsec\n",
"R41_S10 3.5872512072044937e-06 arcsec\n",
"R41_S11 3.4524245646682166e-06 arcsec\n",
"R41_S12 3.8054942971284477e-06 arcsec\n",
"R41_S20 4.067613752630853e-06 arcsec\n",
"R41_S21 4.7327789118198416e-06 arcsec\n",
"R41_S22 6.282363878144554e-06 arcsec\n",
"R42_S00 2.940192988058819e-06 arcsec\n",
"R42_S01 2.7486082414921117e-06 arcsec\n",
"R42_S02 2.5140707996985763e-06 arcsec\n",
"R42_S10 4.002620663265936e-06 arcsec\n",
"R42_S11 3.876953190771095e-06 arcsec\n",
"R42_S12 3.7531969732197456e-06 arcsec\n",
"R42_S20 6.321685493119289e-06 arcsec\n",
"R42_S21 6.846602190265953e-06 arcsec\n",
"R42_S22 7.741481345894385e-06 arcsec\n",
"R43_S00 2.500785225224768e-06 arcsec\n",
"R43_S01 2.939119014135222e-06 arcsec\n",
"R43_S02 5.119268546228302e-06 arcsec\n",
"R43_S10 4.31645631061019e-06 arcsec\n",
"R43_S11 6.415219265070235e-06 arcsec\n",
"R43_S12 1.1143159837594767e-05 arcsec\n",
"R43_S20 1.0250233264242218e-05 arcsec\n",
"R43_S21 1.4361027593873208e-05 arcsec\n",
"R43_S22 7.255121260807699e-06 arcsec\n"
]
}
],
"source": [
"# Repeat focal -> field -> focal for all detectors\n",
"for det in camera:\n",
" wcs = BatoidWCS(\n",
" boresight, rotSkyPos, rotTelPos,\n",
" batoid.Optic.fromYaml(\"LSST_r.yaml\"),\n",
" 620, # nm\n",
" camera, det\n",
" )\n",
" \n",
" fieldx = wcs._field_x + np.deg2rad(2000*0.2/3600)*np.linspace(-1, 1, 800)\n",
" fieldy = wcs._field_y + np.deg2rad(2000*0.2/3600)*np.linspace(-1, 1, 800)\n",
" fieldx, fieldy = np.meshgrid(fieldx, fieldy)\n",
" fieldx = fieldx.ravel()\n",
" fieldy = fieldy.ravel()\n",
"\n",
" focalx, focaly = wcs._fieldToFocalMany(fieldx, fieldy)\n",
" fieldx1, fieldy1 = wcs._focalToFieldMany(focalx, focaly)\n",
" \n",
" print(det.getName(), 3600*np.rad2deg(np.median(np.abs(fieldx1-fieldx))), \"arcsec\")"
]
},
{
"cell_type": "code",
"execution_count": 1049,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.2471421630143595e-05 pixel\n",
"4.088101013621781e-05 pixel\n",
"5.574985804668309e-06 arcsec\n",
"7.499402895974616e-06 arcsec\n"
]
}
],
"source": [
"# Try pixel -> radec -> pixel -> radec roundtrip\n",
"x = np.random.uniform(-50, 4050, size=100000)\n",
"y = np.random.uniform(-50, 4050, size=100000)\n",
"ra, dec = wcs.xyToradec(x, y, units='rad')\n",
"x1, y1 = wcs.radecToxy(ra, dec, units='rad')\n",
"ra1, dec1 = wcs.xyToradec(x1, y1, units='rad')\n",
"\n",
"print(np.median(np.abs(x-x1)), \"pixel\")\n",
"print(np.median(np.abs(y-y1)), \"pixel\")\n",
"print(206265*np.median(np.abs(ra-ra1)), \"arcsec\")\n",
"print(206265*np.median(np.abs(dec-dec1)), \"arcsec\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment