Coverage for python / lsst / cbp / coordUtils.py: 19%
97 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-23 08:26 +0000
1# This file is part of cbp.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21"""Coordinate conversion functions"""
23__all__ = ["fieldAngleToVector", "vectorToFieldAngle", "pupilPositionToVector",
24 "computeShiftedPlanePos", "convertVectorFromBaseToPupil", "convertVectorFromPupilToBase",
25 "computeAzAltFromBasePupil", "getFlippedPos", "rotate2d"]
27import math
29import numpy as np
31from lsst.sphgeom import Vector3d
32from lsst.geom import SpherePoint, radians
34# Error data from computeAzAltFromBasePupil.
35# Error is determined by calling convertVectorFromPupilToBase
36# using the az/alt computed by computeAzAltFromBasePupil
37# and measuring the separation between that base vector and the input.
38_RecordError = False # record errors?
39_ErrorLimitArcsec = None # Only record errors larger than this.
40_ErrorList = None # A list of tuples: (errorArcsec, vectorBase, vectorPupil).
42ZeroSpherePoint = SpherePoint(0, 0, radians)
45def startRecordingErrors(errorLimit=1e-11*radians):
46 """Start recording numeric errors in computeAzAltFromBasePupil
47 and reset the error list.
49 Parameters
50 ----------
51 errorLimit : `lsst.geom.Angle`
52 Only errors larger than this limit are recorded.
53 """
54 global _RecordError, _ErrorLimitArcsec, _ErrorList
55 _RecordError = True
56 _ErrorLimitArcsec = errorLimit.asArcseconds()
57 _ErrorList = []
60def stopRecordingErrors():
61 """Stop recording numeric errors in computeAzAltFromBasePupil.
62 """
63 global _RecordError
64 _RecordError = False
67def getRecordedErrors():
68 """Get recorded numeric errors in computeAzAltFromBasePupil,
69 sorted by increasing error.
71 Returns
72 -------
73 errorList : `list` of tuples
74 Recorded error as a list of tuples containing:
75 - |error| in radians
76 - vectorBase
77 - vectorPupil
78 """
79 if _ErrorList is None:
80 raise RuntimeError("Errors were never recorded")
81 return sorted(_ErrorList, key=lambda elt: elt[0])
84def getFlippedPos(xyPos, flipX):
85 """Get a 2-dimensional position with the x axis properly flipped.
87 Parameters
88 ----------
89 xyPos : pair of `float`
90 Position to rotate.
91 flipX : `bool`
92 True if the x axis should be flipped (negated).
94 Returns
95 -------
96 xyResult : pair of `float`
97 ``xyPos`` with the x axis flipped or not, according to ``flipX``
98 """
99 return (-xyPos[0], xyPos[1]) if flipX else xyPos
102def fieldAngleToVector(xyrad, flipX):
103 """Convert a pupil field angle to a pupil unit vector.
105 Parameters
106 ----------
107 xyrad : `tuple` of 2 `float`
108 x,y :ref:`pupil field angle <lsst.cbp.pupil_field_angle>`
109 (radians).
111 Returns
112 -------
113 vector : `numpy.array` of 3 `float`
114 A unit vector in the :ref:`pupil frame <lsst.cbp.pupil_frame>`.
115 """
116 assert len(xyrad) == 2
117 xyradRightHanded = getFlippedPos(xyrad, flipX=flipX)
118 amount = math.hypot(*xyradRightHanded)*radians
119 bearing = math.atan2(xyradRightHanded[1], xyradRightHanded[0])*radians
120 return np.array(ZeroSpherePoint.offset(bearing=bearing, amount=amount).getVector())
123def vectorToFieldAngle(vec, flipX):
124 """Convert a vector to a pupil field angle.
126 Parameters
127 ----------
128 vec : sequence of 3 `float`
129 3-dimensional vector in the :ref:`pupil frame
130 <lsst.cbp.pupil_frame>`; the magnitude is ignored,
131 but must be large enough to compute an accurate unit vector.
132 flipX : bool
133 Set True if the x axis of the focal plane is flipped
134 with respect to the pupil.
136 Returns
137 -------
138 fieldAngle : `tuple` of 2 `float`
139 x,y :ref:`pupil field angle <lsst.cbp.pupil_field_angle>`
140 (radians).
141 """
142 sp = SpherePoint(Vector3d(*vec))
143 amountRad = ZeroSpherePoint.separation(sp).asRadians()
144 bearingRad = ZeroSpherePoint.bearingTo(sp).asRadians()
145 xyrad = (amountRad*math.cos(bearingRad),
146 amountRad*math.sin(bearingRad))
147 return getFlippedPos(xyrad, flipX=flipX)
150def pupilPositionToVector(xyPos, flipX):
151 """Convert a pupil plane position to a 3D vector.
153 Parameters
154 ----------
155 xyPos : sequence of 2 `float`
156 :ref:`pupil plane position <lsst.cbp.pupil_position>` (mm).
157 flipX : `bool`
158 True if the x axis of the position is flipped (negated).
160 Returns
161 -------
162 vector : sequence of 3 `float`
163 3-dimensional vector in the
164 :ref:`pupil frame <lsst.cbp.pupil_frame>`:
165 x = 0, y = plane position x, z = plane position y.
166 """
167 xyPosRightHanded = getFlippedPos(xyPos, flipX=flipX)
168 return np.array([0, xyPosRightHanded[0], xyPosRightHanded[1]])
171def computeShiftedPlanePos(planePos, fieldAngle, shift):
172 """Compute the plane position of a vector on a plane
173 shifted along the optical axis.
175 Parameters
176 ----------
177 planePos : pair of `float`
178 Plane position at which the vector intersects the plane (x, y).
179 fieldAngle : pair of `float`
180 Field angle of vector (x, y radians).
181 shift : `float`
182 Amount by which the new plane is shifted along the optical axis.
183 If ``shift`` and both components of ``fieldAngle``
184 are positive then both axes of the shifted plane position
185 will be larger (more positive) than ``planePos``.
187 Returns
188 -------
189 shiftedPlanePos : tuple of `float`
190 Plane position at which vector intersects the shifted plane (x, y).
192 Notes
193 -----
194 `flipX` is not an input because for this computation it does not matter
195 if the x axis is flipped: fieldAngle and planePos are either both
196 flipped or not, and that cancels out the effect.
197 """
198 # unit vector y = plane x, unit vector z = plane y
199 # (ignoring flipped x axis, which cancels out)
200 unitVector = fieldAngleToVector(fieldAngle, False)
201 dxy = [shift*unitVector[i]/unitVector[0] for i in (1, 2)]
202 return tuple(planePos[i] + dxy[i] for i in range(2))
205def convertVectorFromBaseToPupil(vectorBase, pupilAzAlt):
206 """Given a vector in base coordinates and the pupil pointing,
207 compute the vector in pupil coordinates.
209 Parameters
210 ----------
211 vectorBase : sequence of 3 `float`
212 3-dimensional vector in the :ref:`base frame
213 <lsst.cbp.base_frame>`.
214 pupilAzAlt : `lsst.geom.SpherePoint`
215 Pointing of the pupil frame as :ref:`internal azimuth, altitude
216 <lsst.cbp.internal_angles>`.
218 Returns
219 -------
220 vectorPupil : `np.array` of 3 `float`
221 3-dimensional vector in the :ref:`pupil frame
222 <lsst.cbp.pupil_frame>`.
224 Notes
225 -----
226 This could be implemented as the following Euler angle
227 rotation matrix, which is:
228 - first rotate about the z axis by azimuth
229 - then rotate about the rotated -y axis by altitude
230 - there is no third rotation
232 c1*c2 -s1 -c1*s2
233 c2*s1 c1 s1s2
234 s1 0 c2
236 where angle 1 = azimuth, angle 2 = altitude,
237 sx = sine(angle x) and cx = cosine(angle x).
239 Knowing this matrix is helpful, e.g. for math inside
240 computeAzAltFromBasePupil.
241 """
242 vectorMag = np.linalg.norm(vectorBase)
243 vectorSpBase = SpherePoint(Vector3d(*vectorBase))
245 telBase = SpherePoint(0, 0, radians)
247 # rotate vector around base z axis by -daz
248 daz = pupilAzAlt[0] - telBase[0]
249 zaxis = SpherePoint(Vector3d(0, 0, 1))
250 vectorSpRot1 = vectorSpBase.rotated(axis=zaxis, amount=-daz)
252 # rotate vector around pupil -y axis by -dalt
253 dalt = pupilAzAlt[1] - telBase[1]
254 negYAxis = SpherePoint(Vector3d(0, -1, 0))
255 vectorSpPupil = vectorSpRot1.rotated(axis=negYAxis, amount=-dalt)
256 return np.array(vectorSpPupil.getVector()) * vectorMag
259def convertVectorFromPupilToBase(vectorPupil, pupilAzAlt):
260 """Given a vector in pupil coordinates and the pupil pointing,
261 compute the vector in base coords.
263 Parameters
264 ----------
265 vectorPupil : sequence of 3 `float`
266 3-dimesional vector in the :ref:`pupil frame
267 <lsst.cbp.pupil_frame>`.
268 pupilAzAlt : `lsst.geom.SpherePoint`
269 Pointing of the pupil frame as :ref:`internal azimuth, altitude
270 <lsst.cbp.internal_angles>`.
272 Returns
273 -------
274 vectorBase : `np.array` of 3 `float`
275 3-dimensional vector in the :ref:`base frame
276 <lsst.cbp.base_frame>`.
277 """
278 vectorMag = np.linalg.norm(vectorPupil)
279 vectorSpPupil = SpherePoint(Vector3d(*vectorPupil))
281 telBase = SpherePoint(0, 0, radians)
283 # rotate vector around pupil -y axis by dalt
284 dalt = pupilAzAlt[1] - telBase[1]
285 negYAxis = SpherePoint(Vector3d(0, -1, 0))
286 vectorSpRot1 = vectorSpPupil.rotated(axis=negYAxis, amount=dalt)
288 # rotate that around base z axis by daz
289 daz = pupilAzAlt[0] - telBase[0]
290 zaxis = SpherePoint(Vector3d(0, 0, 1))
291 vectorSpBase = vectorSpRot1.rotated(axis=zaxis, amount=daz)
292 return np.array(vectorSpBase.getVector()) * vectorMag
295def computeAzAltFromBasePupil(vectorBase, vectorPupil):
296 """Compute az/alt from a vector in the base frame
297 and the same vector in the pupil frame.
299 Parameters
300 ----------
301 vectorBase : `iterable` of three `float`
302 3-dimensional vector in the :ref:`base frame
303 <lsst.cbp.base_frame>`.
304 vectorPupil : `iterable` of `float`
305 The same vector in the :ref:`pupil frame <lsst.cbp.pupil_frame>`.
306 This vector should be within 45 degrees or so of the optical axis
307 for accurate results.
309 Returns
310 -------
311 pupilAzAlt : `lsst.geom.SpherePoint`
312 Pointing of the pupil frame as :ref:`internal azimuth, altitude
313 <lsst.cbp.internal_angles>`.
315 Raises
316 ------
317 ValueError
318 If vectorPupil x <= 0
320 Notes
321 -----
322 The magnitude of each vector is ignored, except that a reasonable
323 magnitude is required in order to compute an accurate unit vector.
324 """
325 if vectorPupil[0] <= 0:
326 raise ValueError("vectorPupil x must be > 0: {}".format(vectorPupil))
328 # Compute telescope altitude using:
329 #
330 # base z = sin(alt) pupil x + cos(alt) pupil z
331 #
332 # One way to derive this is from the last row of an Euler rotation
333 # matrix listed in the comments for convertVectorFromBaseToPupil.
334 spBase = SpherePoint(Vector3d(*vectorBase))
335 spPupil = SpherePoint(Vector3d(*vectorPupil))
336 xb, yb, zb = spBase.getVector()
337 xp, yp, zp = spPupil.getVector()
338 factor = 1 / math.fsum((xp**2, zp**2))
339 addend1 = xp * zb
340 addend2 = zp * math.sqrt(math.fsum((xp**2, zp**2, -zb**2)))
341 if zp == 0:
342 sinAlt = zb/xp
343 else:
344 sinAlt = factor*(addend1 - addend2)
345 alt = math.asin(sinAlt)*radians
347 # Consider the spherical triangle connecting the telescope pointing
348 # (pupil frame x axis), the vector, and zenith (the base frame z axis).
349 # The length of all sides is known:
350 # - sideA is the side connecting the vector to the pupil frame x axis
351 # (since 0, 0 is a unit vector pointing along pupil frame x);
352 # - sideB is the side connecting telescope pointing to the zenith
353 # - sideC is the side connecting the vector to the zenith
354 #
355 # Solve for angleA, the angle between the sides at the zenith;
356 # that angle is the difference in azimuth between the telescope pointing
357 # and the azimuth of the base vector.
358 sideA = SpherePoint(0, 0, radians).separation(spPupil).asRadians()
359 sideB = math.pi/2 - alt.asRadians()
360 sideC = math.pi/2 - spBase[1].asRadians()
362 # sideA can be small or zero so use a half angle formula
363 # sides B and C will always be well away from 0 and 180 degrees
364 semiPerimeter = 0.5*math.fsum((sideA, sideB, sideC))
365 sinHalfAngleA = math.sqrt(math.sin(semiPerimeter - sideB) * math.sin(semiPerimeter - sideC)
366 / (math.sin(sideB) * math.sin(sideC)))
367 daz = 2*math.asin(sinHalfAngleA)*radians
368 if spPupil[0].wrapCtr() > 0:
369 daz = -daz
370 az = spBase[0] + daz
371 if _RecordError: # to study sources of numerical imprecision
372 sp = SpherePoint(az, alt)
373 vectorBaseRT = convertVectorFromPupilToBase(vectorPupil, sp)
374 errorArcsec = SpherePoint(Vector3d(*vectorBaseRT)).separation(
375 SpherePoint(Vector3d(*vectorBase))).asArcseconds()
376 if errorArcsec > _ErrorLimitArcsec:
377 _ErrorList.append((errorArcsec, vectorBase, vectorPupil))
378 return SpherePoint(az, alt)
381def rotate2d(pos, angle):
382 """Rotate a 2-dimensional position by a given angle.
384 Parameters
385 ----------
386 pos : pair of `float`
387 Position to rotate.
388 angle : `lsst.geom.Angle`
389 Amount of rotation.
391 Returns
392 -------
393 rotPos : pair of `float`
394 Rotated position.
396 Examples
397 --------
398 ``rotate2d((1, 2), 90*lsst.geom.degrees, False)`` returns `(-2, 1)`.
399 """
400 angRad = angle.asRadians()
401 sinAng = math.sin(angRad)
402 cosAng = math.cos(angRad)
403 return (
404 cosAng*pos[0] - sinAng*pos[1],
405 sinAng*pos[0] + cosAng*pos[1]
406 )