Coverage for python / lsst / summit / utils / simonyi / mountAnalysis.py: 0%
356 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-01 09:03 +0000
1# This file is part of summit_utils.
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/>.
22from __future__ import annotations
24__all__ = [
25 "calculateMountErrors",
26 "plotMountErrors",
27 "MountErrors",
28 "getLinearRates",
29 "getAltAzOverPeriod",
30 "calculateHexRms",
31]
33import copy
34import logging
35from dataclasses import dataclass
36from typing import TYPE_CHECKING
37from zoneinfo import ZoneInfo
39import astropy.units as u
40import matplotlib.dates as mdates
41import matplotlib.pyplot as plt
42import numpy as np
43from astropy.coordinates import AltAz, EarthLocation, SkyCoord
44from matplotlib.dates import num2date
45from matplotlib.ticker import FuncFormatter
47from lsst.summit.utils.dateTime import dayObsIntToString
48from lsst.summit.utils.tmaUtils import filterBadValues
49from lsst.utils.plotting.figures import make_figure
51from .mountData import getAzElRotHexDataForExposure
53if TYPE_CHECKING:
54 from astropy.time import Time
55 from lsst_efd_client import EfdClient
56 from matplotlib.figure import Figure
58 from lsst.daf.butler import DimensionRecord
60 from .mountData import MountData
63NON_TRACKING_IMAGE_TYPES = ["BIAS", "FLAT"]
65COMCAM_ANGLE_TO_EDGE_OF_FIELD_ARCSEC = 1800.0
66LSSTCAM_ANGLE_TO_EDGE_OF_FIELD_ARCSEC = 8500.0
68# These levels determine the colouring of the cells in the RubinTV.
69# Yellow for warning level, red for bad level
70MOUNT_IMAGE_WARNING_LEVEL = 0.05
71MOUNT_IMAGE_BAD_LEVEL = 0.10 # and red for this
73N_REPLACED_WARNING_LEVEL = 999999 # fill these values in once you've spoken to Craig and Brian
74N_REPLACED_BAD_LEVEL = 999999 # fill these values in once you've spoken to Craig and Brian
76SIMONYI_LOCATION = EarthLocation.of_site("Rubin:Simonyi")
77EARTH_ROTATION = 15.04106858 # degrees/hour
80@dataclass
81class MountErrors:
82 azRms: float
83 elRms: float
84 rotRms: float
85 camHexRms: float
86 m2HexRms: float
87 imageAzRms: float
88 imageElRms: float
89 imageRotRms: float
90 imageImpactRms: float
91 residualFiltering: bool
92 nReplacedAz: int
93 nReplacedEl: int
96def tickFormatter(value: float, tick_number: float) -> str:
97 # Convert the value to a string without subtracting large numbers
98 # tick_number is unused.
99 return f"{value:.2f}"
102def calculateMountErrors(
103 expRecord: DimensionRecord,
104 client: EfdClient,
105 maxDelta: float = 0.1,
106 doFilterResiduals: bool = False,
107 useMockPointingModelResidualsAboveAzEl: float = 10.0,
108 useMockPointingModelResidualsAboveRot: float = 15.0,
109) -> tuple[MountErrors, MountData] | tuple[None, None]:
110 """Queries the EFD over a given exposure and calculates the RMS errors
111 for the axes, optionally using a pointing model to calculate residuals.
113 Parameters
114 ----------
115 expRecord : `DimensionRecord`
116 The exposure record containing the necessary fields for calculations.
117 client : `EfdClient`
118 The EFD client to query for mount data.
119 maxDelta : `float`, optional
120 The maximum delta for filtering bad values, by default 0.1.
121 doFilterResiduals : `bool`, optional
122 Whether to filter residuals.
123 useMockPointingModelResidualsAboveAzEl : `float`, optional
124 The threshold above which to use the mock pointing model residuals, as
125 an RMS, in arcseconds, for the azimuth and elevation axes.
126 useMockPointingModelResidualsAboveRot : `float`, optional
127 The threshold above which to use the mock pointing model residuals, as
128 an RMS, in arcseconds, for the rotator.
130 Returns
131 -------
132 tuple[MountErrors, MountData] | tuple[None, None]
133 A tuple containing the mount errors and mount data, or (None, None) if
134 the exposure type is non-tracking.
135 """
136 logger = logging.getLogger(__name__)
138 imgType = expRecord.observation_type.upper()
139 if imgType in NON_TRACKING_IMAGE_TYPES:
140 logger.info(f"Skipping mount torques for non-tracking image type {imgType} for {expRecord.id}")
141 return None, None
143 mountData = getAzElRotHexDataForExposure(client, expRecord)
145 elevation = 90 - expRecord.zenith_angle
147 azError = mountData.azimuthData["azError"].to_numpy()
148 elError = mountData.elevationData["elError"].to_numpy()
149 rotError = mountData.rotationData["rotError"].to_numpy()
150 nReplacedAz = 0
151 nReplacedEl = 0
152 if doFilterResiduals:
153 # Filtering out bad values
154 nReplacedAz = filterBadValues(azError, maxDelta)
155 nReplacedEl = filterBadValues(elError, maxDelta)
156 mountData.azimuthData["azError"] = azError
157 mountData.elevationData["elError"] = elError
159 # Calculate the linear demand model
160 if len(mountData.azimuthData) == len(mountData.elevationData):
161 azModelValues, elModelValues = getAltAzOverPeriod(expRecord, nPoints=len(mountData.azimuthData))
162 else:
163 azModelValues, _ = getAltAzOverPeriod(expRecord, nPoints=len(mountData.azimuthData))
164 _, elModelValues = getAltAzOverPeriod(expRecord, nPoints=len(mountData.elevationData))
166 _, _, rotRate = getLinearRates(expRecord)
168 azimuthData = mountData.azimuthData
169 azValues = np.asarray(azimuthData["actualPosition"])
170 azMedian = np.median(azValues)
171 azModelMedian = np.median(azModelValues)
172 # subtract the overall offset
173 azModelValues -= azModelMedian - azMedian
174 azimuthData["linearModel"] = azModelValues
175 azLinearError = (azValues - azModelValues) * 3600
176 azLinearRms = np.sqrt(np.mean(azLinearError * azLinearError))
177 if azLinearRms > useMockPointingModelResidualsAboveAzEl:
178 logger.warning(
179 f"Azimuth pointing model RMS error {azLinearRms:.3f} arcsec is above threshold of "
180 f"{useMockPointingModelResidualsAboveAzEl:.3f} arcsec, calculating errors vs astropy."
181 )
182 # If linear error is large, replace demand errors with linear error
183 azimuthData["azError"] = azLinearError
185 elevationData = mountData.elevationData
186 elValues = np.asarray(elevationData["actualPosition"])
187 elMedian = np.median(elValues)
188 elModelMedian = np.median(elModelValues)
189 # subtract the overall offset
190 elModelValues -= elModelMedian - elMedian
191 elevationData["linearModel"] = elModelValues
192 elLinearError = (elValues - elModelValues) * 3600
193 elLinearRms = np.sqrt(np.mean(elLinearError * elLinearError))
194 if elLinearRms > useMockPointingModelResidualsAboveAzEl:
195 logger.warning(
196 f"Elevation pointing model RMS error {elLinearRms:.3f} arcsec is above threshold of "
197 f"{useMockPointingModelResidualsAboveAzEl:.3f} arcsec, calculating errors vs astropy."
198 )
199 # If linear error is large, replace demand errors with linear error
200 elevationData["elError"] = elLinearError
202 rotationData = mountData.rotationData
203 rotValues = np.asarray(rotationData["actualPosition"])
204 rotValTimes = np.asarray(rotationData["timestamp"])
205 rotModelValues = np.zeros_like(rotValues)
206 rotMedian = np.median(rotValues)
207 rotTimesMedian = np.median(rotValTimes)
208 rotModelValues = rotMedian + rotRate * (rotValTimes - rotTimesMedian)
209 rotationData["linearModel"] = rotModelValues
210 rotLinearError = (rotValues - rotModelValues) * 3600
211 rotLinearRms = np.sqrt(np.mean(rotLinearError * rotLinearError))
212 if rotLinearRms > useMockPointingModelResidualsAboveRot:
213 logger.warning(
214 f"Rotation pointing model RMS error {rotLinearRms:.3f} arcsec is above threshold of "
215 f"{useMockPointingModelResidualsAboveAzEl:.3f} arcsec, calculating errors vs astropy."
216 )
217 # If linear error is large, replace demand errors with linear error
218 rotationData["rotError"] = rotLinearError
220 azError = mountData.azimuthData["azError"].to_numpy()
221 elError = mountData.elevationData["elError"].to_numpy()
222 rotError = mountData.rotationData["rotError"].to_numpy()
224 azRms = np.sqrt(np.mean(azError * azError))
225 elRms = np.sqrt(np.mean(elError * elError))
226 rotRms = np.sqrt(np.mean(rotError * rotError))
228 # Calculate Image impact RMS
229 imageAzRms = azRms * np.cos(elevation * np.pi / 180.0)
230 imageElRms = elRms
231 imageRotRms = rotRms * LSSTCAM_ANGLE_TO_EDGE_OF_FIELD_ARCSEC * np.pi / 180.0 / 3600.0
232 camHexRms, m2HexRms = calculateHexRms(mountData)
233 # TODO should the hex RMS values be added in quadrature?
234 imageImpactRms = np.sqrt(imageAzRms**2 + imageElRms**2 + imageRotRms**2 + camHexRms**2 + m2HexRms**2)
236 mountErrors = MountErrors(
237 azRms=float(azRms),
238 elRms=float(elRms),
239 rotRms=float(rotRms),
240 camHexRms=float(camHexRms),
241 m2HexRms=float(m2HexRms),
242 imageAzRms=float(imageAzRms),
243 imageElRms=float(imageElRms),
244 imageRotRms=float(imageRotRms),
245 imageImpactRms=float(imageImpactRms),
246 residualFiltering=doFilterResiduals,
247 nReplacedAz=nReplacedAz,
248 nReplacedEl=nReplacedEl,
249 )
251 return (mountErrors, mountData)
254def plotMountErrors(
255 mountData: MountData,
256 mountErrors: MountErrors,
257 figure=None,
258 saveFilename: str = "",
259) -> Figure:
260 mountData = copy.deepcopy(mountData) # Ensure we don't modify the original data
261 mountErrors = copy.deepcopy(mountErrors)
263 imageImpactRms = mountErrors.imageImpactRms
264 expRecord = mountData.expRecord
265 if expRecord is not None:
266 dayObsString = dayObsIntToString(expRecord.day_obs)
267 dataIdString = f"{expRecord.instrument} {dayObsString} - seqNum {expRecord.seq_num}"
268 title = f"{dataIdString} - Exposure time = {expRecord.exposure_time:.1f}s"
269 else:
270 title = "Mount Errors" # if the data is of unknown provenance
272 if figure is None:
273 figure = make_figure(figsize=(12, 8))
274 else:
275 figure.clear()
276 ax = figure.gca()
277 ax.clear()
279 utc = ZoneInfo("UTC")
280 chile_tz = ZoneInfo("America/Santiago")
282 # Function to convert UTC to Chilean time
283 def offset_time_aware(utc_time):
284 # Ensure the time is timezone-aware in UTC
285 if utc_time.tzinfo is None:
286 utc_time = utc.localize(utc_time)
287 return utc_time.astimezone(chile_tz)
289 [[ax1, ax4], [ax2, ax5], [ax3, ax6]] = figure.subplots(
290 3,
291 2,
292 sharex="col",
293 sharey=False,
294 gridspec_kw={
295 "wspace": 0.35,
296 "hspace": 0,
297 "height_ratios": [2.5, 1, 1],
298 "width_ratios": [1.4, 1],
299 "left": 0.07,
300 "right": 0.60,
301 },
302 )
303 [ax7, ax8] = figure.subplots(
304 2,
305 1,
306 sharex="col",
307 sharey=False,
308 gridspec_kw={
309 "wspace": 0.35,
310 "hspace": 0,
311 "height_ratios": [1, 1],
312 "width_ratios": [1],
313 "left": 0.73,
314 "right": 0.94,
315 },
316 )
317 # [ax1, ax4] = [azimuth, rotator]
318 # [ax2, ax5] = [azError, rotError]
319 # [ax3, ax6] = [azTorque, rotTorque]
320 # [ax7, ax8] = [camHex, m2Hex]
322 # Use the native color cycle for the lines. Because they're on
323 # different axes they don't cycle by themselves
324 axs = [ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8]
325 lineColors = [p["color"] for p in plt.rcParams["axes.prop_cycle"]]
326 nColors = len(lineColors)
327 colorCounter = 0
329 ax1.plot(
330 mountData.azimuthData["actualPosition"],
331 label="Azimuth position",
332 c=lineColors[colorCounter % nColors],
333 )
334 colorCounter += 1
335 ax1.plot(
336 mountData.azimuthData["linearModel"],
337 label="Azimuth linear model",
338 ls="--",
339 c=lineColors[colorCounter % nColors],
340 )
341 colorCounter += 1
342 ax1.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
343 ax1.set_ylabel("Azimuth (degrees)")
345 ax1_twin = ax1.twinx()
346 ax1_twin.plot(
347 mountData.elevationData["actualPosition"],
348 label="Elevation position",
349 c=lineColors[colorCounter % nColors],
350 )
351 colorCounter += 1
352 ax1_twin.plot(
353 mountData.elevationData["linearModel"],
354 label="Elevation linear model",
355 ls="--",
356 c=lineColors[colorCounter % nColors],
357 )
358 colorCounter += 1
360 ax2.plot(
361 mountData.azimuthData["azError"],
362 label="Azimuth tracking error",
363 c=lineColors[colorCounter % nColors],
364 )
365 colorCounter += 1
366 ax2.plot(
367 mountData.elevationData["elError"],
368 label="Elevation tracking error",
369 c=lineColors[colorCounter % nColors],
370 )
371 colorCounter += 1
372 ax2.axhline(0.01, ls="-.", color="black")
373 ax2.axhline(-0.01, ls="-.", color="black")
374 ax2.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
375 ax2.set_ylabel("Tracking error (arcsec)")
376 ax2.set_xticks([]) # remove x tick labels on the hidden upper x-axis
377 ax2.set_ylim(-0.05, 0.05)
378 ax2.set_yticks([-0.04, -0.02, 0.0, 0.02, 0.04])
379 ax2.legend(loc="lower center")
380 ax2.text(0.1, 0.9, f"Image impact = {imageImpactRms:.3f} arcsec (w/ rot&hex).", transform=ax2.transAxes)
381 if mountErrors.residualFiltering:
382 ax2.text(
383 0.1,
384 0.8,
385 (
386 f"{mountErrors.nReplacedAz} bad az values and "
387 f"{mountErrors.nReplacedEl} bad el values were replaced"
388 ),
389 transform=ax2.transAxes,
390 )
391 ax3_twin = ax3.twinx()
392 ax3.plot(
393 mountData.azimuthData["actualTorque"],
394 label="Azimuth torque",
395 c=lineColors[colorCounter % nColors],
396 )
397 colorCounter += 1
398 ax3_twin.plot(
399 mountData.elevationData["actualTorque"],
400 label="Elevation torque",
401 c=lineColors[colorCounter % nColors],
402 )
403 colorCounter += 1
404 ax3.set_ylabel("Azimuth torque (Nm)")
405 ax3.set_xlabel("Time (UTC)") # yes, it really is UTC, matplotlib converts this automatically!
407 # put the ticks at an angle, and right align with the tick marks
408 ax3.set_xticks(ax3.get_xticks()) # needed to supress a user warning
409 xlabels = ax3.get_xticks()
410 ax3.set_xticklabels(xlabels)
411 ax3.tick_params(axis="x", rotation=45)
412 ax3.xaxis.set_major_locator(mdates.AutoDateLocator())
413 ax3.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
415 ax4.plot(
416 mountData.rotationData["actualPosition"],
417 label="Rotator position",
418 c=lineColors[colorCounter % nColors],
419 )
420 colorCounter += 1
421 ax4.plot(
422 mountData.rotationData["linearModel"],
423 label="Rotator linearModel",
424 ls="--",
425 c=lineColors[colorCounter % nColors],
426 )
427 colorCounter += 1
428 ax4.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
429 ax4.yaxis.tick_right()
430 ax4.set_ylabel("Rotator angle (degrees)")
431 ax4.yaxis.set_label_position("right")
432 ax5.plot(
433 mountData.rotationData["rotError"],
434 c=lineColors[colorCounter % nColors],
435 )
437 colorCounter += 1
438 ax5.axhline(0.1, ls="-.", color="black")
439 ax5.axhline(-0.1, ls="-.", color="black")
440 ax5.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
441 ax5.set_ylabel("Tracking error (arcsec)")
442 ax5.tick_params(labelbottom=False) # Hide x-axis tick labels without removing ticks
443 ax5.set_ylim(-0.5, 0.5)
444 ax5.set_yticks([-0.4, -0.2, 0.0, 0.2, 0.4])
445 ax5.yaxis.tick_right()
446 ax5.yaxis.set_label_position("right")
448 ax6.plot(mountData.rotationTorques["torque0"], label="Torque0", c=lineColors[colorCounter % nColors])
449 colorCounter += 1
450 ax6.plot(mountData.rotationTorques["torque1"], label="Torque1", c=lineColors[colorCounter % nColors])
451 ax6.set_xlabel("Time (UTC)") # yes, it really is UTC, matplotlib converts this automatically!
452 # put the ticks at an angle, and right align with the tick marks
453 ax6.set_xticks(ax6.get_xticks()) # needed to supress a user warning
454 xlabels = ax6.get_xticks()
455 ax6.set_xticklabels(xlabels)
456 ax6.tick_params(axis="x", rotation=45)
457 ax6.xaxis.set_major_locator(mdates.AutoDateLocator())
458 ax6.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
459 ax6.yaxis.tick_right()
460 ax6.yaxis.set_label_position("right")
461 ax6.legend()
463 hexNames = ["X", "Y", "Z", "U", "V", "W"]
464 for i in [0, 1, 2]:
465 camhex = mountData.camhexData[f"position{i}"]
466 camhex -= np.median(camhex)
467 ax7.plot(
468 camhex,
469 label=hexNames[i],
470 c=lineColors[colorCounter % nColors],
471 )
472 colorCounter += 1
473 ax7.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
474 ax7.set_ylabel("CamHex XYZ(micron) (minus median)")
475 ax7.legend()
476 ax7_twin = ax7.twinx()
477 for i in [3, 4]:
478 camhex = mountData.camhexData[f"position{i}"]
479 camhex *= 3600.0 # convert to arcseconds
480 camhex -= np.median(camhex)
481 ax7_twin.plot(
482 camhex,
483 label=hexNames[i],
484 c=lineColors[colorCounter % nColors],
485 )
486 colorCounter += 1
487 ax7_twin.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
488 ax7_twin.set_ylabel("CamHex UV(arcsec) (minus median)")
489 # ax7_twin.yaxis.tick_right()
490 # ax7.yaxis.set_label_position("right")
491 ax7_twin.legend()
493 for i in [0, 1, 2]:
494 m2hex = mountData.m2hexData[f"position{i}"]
495 m2hex -= np.median(m2hex)
496 ax8.plot(
497 m2hex,
498 label=hexNames[i],
499 c=lineColors[colorCounter % nColors],
500 )
501 colorCounter += 1
502 ax8.legend()
503 ax8.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
504 ax8.set_ylabel("M2Hex XYZ(micron) (minus median)")
505 ax8.set_xlabel("Time (UTC)") # yes, it really is UTC, matplotlib converts this automatically!
506 # put the ticks at an angle, and right align with the tick marks
507 ax8.set_xticks(ax8.get_xticks()) # needed to supress a user warning
508 xlabels = ax8.get_xticks()
509 ax8.set_xticklabels(xlabels)
510 ax8.tick_params(axis="x", rotation=45)
511 ax8.xaxis.set_major_locator(mdates.AutoDateLocator())
512 ax8.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
513 ax8_twin = ax8.twinx()
514 for i in [3, 4]:
515 m2hex = mountData.m2hexData[f"position{i}"]
516 m2hex *= 3600.0 # Convert to arcseconds
517 m2hex -= np.median(m2hex)
518 ax8_twin.plot(
519 m2hex,
520 label=hexNames[i],
521 c=lineColors[colorCounter % nColors],
522 )
523 colorCounter += 1
524 ax8_twin.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
525 ax8_twin.set_ylabel("M2Hex UV(arcsec) (minus median)")
526 # ax7_twin.yaxis.tick_right()
527 # ax7.yaxis.set_label_position("right")
528 ax8_twin.legend()
530 ax1_twin.yaxis.set_major_formatter(FuncFormatter(tickFormatter))
531 ax1_twin.set_ylabel("Elevation (degrees)")
532 ax1.tick_params(labelbottom=False) # Hide x-axis tick labels without removing ticks
533 # combine the legends and put inside the plot
534 handles1a, labels1a = ax1.get_legend_handles_labels()
535 handles1b, labels1b = ax1_twin.get_legend_handles_labels()
536 handles2a, labels2a = ax3.get_legend_handles_labels()
537 handles2b, labels2b = ax3_twin.get_legend_handles_labels()
538 handles = handles1a + handles1b + handles2a + handles2b
539 labels = labels1a + labels1b + labels2a + labels2b
540 # ax2 is "in front" of ax1 because it has the vlines plotted on it, and
541 # vlines are on ax2 so that they appear at the bottom of the legend, so
542 # make sure to plot the legend on ax2, otherwise the vlines will go on
543 # top of the otherwise-opaque legend.
544 ax1_twin.legend(handles, labels, facecolor="white", framealpha=1)
546 ax1.set_title("Azimuth and Elevation")
547 ax4.set_title("Rotator")
548 ax7.set_title("Hexapods")
549 ax4.legend()
550 figure.subplots_adjust(top=0.85) # Adjust the top margin to make room for the suptitle
551 figure.suptitle(title, fontsize=14, y=1.04) # Adjust y to move the title up
553 # Create the upper axis for Chilean time
554 ax1_twiny = ax1.twiny()
555 ax1_twiny.set_xlim(ax1.get_xlim()) # Set the limits of the upper axis to match the lower axis
556 utcTicks = ax1.get_xticks() # Use the same ticks as the lower UTC axis
557 utcTickLabels = [num2date(tick, tz=utc) for tick in utcTicks]
558 chileTickLabels = [offset_time_aware(label) for label in utcTickLabels]
559 # Set the same tick positions but with Chilean time labels
560 ax1_twiny.set_xticks(utcTicks)
561 ax1_twiny.set_xticklabels([tick.strftime("%H:%M:%S") for tick in chileTickLabels])
562 ax1_twiny.tick_params(axis="x", rotation=45)
563 ax1_twiny.set_xlabel("Time (Chilean Time)")
565 ax4_twiny = ax4.twiny()
566 ax4_twiny.set_xlim(ax4.get_xlim()) # Set the limits of the upper axis to match the lower axis
567 utcTicks = ax4.get_xticks() # Use the same ticks as the lower UTC axis
568 utcTickLabels = [num2date(tick, tz=utc) for tick in utcTicks]
569 chileTickLabels = [offset_time_aware(label) for label in utcTickLabels]
570 # Set the same tick positions but with Chilean time labels
571 ax4_twiny.set_xticks(utcTicks)
572 ax4_twiny.set_xticklabels([tick.strftime("%H:%M:%S") for tick in chileTickLabels])
573 ax4_twiny.tick_params(axis="x", rotation=45)
574 ax4_twiny.set_xlabel("Time (Chilean Time)")
576 ax7_twiny = ax7.twiny()
577 ax7_twiny.set_xlim(ax7.get_xlim()) # Set the limits of the upper axis to match the lower axis
578 utcTicks = ax7.get_xticks() # Use the same ticks as the lower UTC axis
579 utcTickLabels = [num2date(tick, tz=utc) for tick in utcTicks]
580 chileTickLabels = [offset_time_aware(label) for label in utcTickLabels]
581 # Set the same tick positions but with Chilean time labels
582 ax7_twiny.set_xticks(utcTicks)
583 ax7_twiny.set_xticklabels([tick.strftime("%H:%M:%S") for tick in chileTickLabels])
584 ax7_twiny.tick_params(axis="x", rotation=45)
585 ax7_twiny.set_xlabel("Time (Chilean Time)")
587 # Add exposure start and end:
588 for ax in axs:
589 if expRecord is not None:
590 # assert expRecord is not None, "expRecord is None"
591 ax.axvline(expRecord.timespan.begin.utc.datetime, ls="--", color="green")
592 ax.axvline(expRecord.timespan.end.utc.datetime, ls="--", color="red")
594 if saveFilename:
595 figure.savefig(saveFilename, bbox_inches="tight")
597 return figure
600def getLinearRates(expRecord: DimensionRecord) -> tuple[float, float, float]:
601 """Calculate the linear rates of motion for az, el, and rotation during an
602 exposure.
604 The rates are calculated based on the tracking RA and Dec, azimuth, zenith
605 angle, and the exposure timespan. The rates are returned in degrees per
606 second.
608 Parameters
609 ----------
610 expRecord : `DimensionRecord`
611 The exposure record containing the necessary fields for calculations.
613 Returns
614 -------
615 azRate, elRate, rotRate: `tuple`[`float`, `float`, `float`]
616 The azimuth rate, elevation rate, and rotator rate in degrees per
617 second.
618 """
619 begin: Time = expRecord.timespan.begin
620 end: Time = expRecord.timespan.end
621 dT: float = (expRecord.timespan.end - expRecord.timespan.begin).value * 86400.0
622 rotRate = (
623 -EARTH_ROTATION
624 * np.cos(SIMONYI_LOCATION.lat.rad)
625 * np.cos(expRecord.azimuth * u.deg)
626 / np.cos((90.0 - expRecord.zenith_angle) * u.deg)
627 / 3600.0
628 )
629 skyLocation = SkyCoord(expRecord.tracking_ra * u.deg, expRecord.tracking_dec * u.deg)
630 altAz1 = AltAz(obstime=begin, location=SIMONYI_LOCATION)
631 altAz2 = AltAz(obstime=end, location=SIMONYI_LOCATION)
632 obsAltAz1 = skyLocation.transform_to(altAz1)
633 obsAltAz2 = skyLocation.transform_to(altAz2)
634 elRate = float((obsAltAz2.alt.deg - obsAltAz1.alt.deg) / dT)
635 azRate = float((obsAltAz2.az.deg - obsAltAz1.az.deg) / dT)
637 # All rates are in degrees / second
638 return azRate, elRate, float(rotRate.value)
641def getAltAzOverPeriod(
642 expRecord: DimensionRecord,
643 nPoints: int,
644) -> tuple[np.ndarray, np.ndarray]:
645 """Get the AltAz coordinates over a period.
647 Parameters
648 ----------
649 begin : `Time`
650 The beginning of the period.
651 end : `Time`
652 The end of the period.
653 target : `SkyCoord`
654 The sky coordinates to track.
655 nPoints : `int`, optional
656 The number of points to sample, by default 100.
658 Returns
659 -------
660 tuple[np.ndarray, np.ndarray]
661 The azimuth and elevation coordinates in degrees.
662 """
663 begin = expRecord.timespan.begin
664 end = expRecord.timespan.end
665 times = begin + (end - begin) * np.linspace(0, 1, nPoints)
666 target = SkyCoord(expRecord.tracking_ra * u.deg, expRecord.tracking_dec * u.deg)
667 altAzFrame = AltAz(obstime=times, location=SIMONYI_LOCATION)
668 targetAltAz = target.transform_to(altAzFrame)
669 az = targetAltAz.az
670 if abs(az[0].degree) < 90.0:
671 az_wrapped = az.wrap_at(180.0 * u.deg)
672 else:
673 az_wrapped = az.wrap_at(0.0 * u.deg)
674 return az_wrapped.degree, targetAltAz.alt.degree
677def calculateHexRms(mountData: MountData) -> tuple[float, float]:
678 """Calculate the image impact of hexapod motions.
680 Parameters
681 ----------
682 mountData : MountData
683 The EFD data associated with the exposure
685 Returns
686 -------
687 tuple[float, float]
688 The image motions associated with the CamHex and M2Hex motions.
689 """
691 # The below image motion coefficients were calculated
692 # with a Batoid simulation by Josh Meyers
693 camHexXY = 1.00 # microns(image) / micron(hexapod)
694 camHexUV = 4.92 # microns(image) / arcsecond(hexapod)
695 m2HexXY = 1.13 # microns(image) / micron(hexapod)
696 m2HexUV = 37.26 # microns(image) / arcsecond(hexapod)
698 # Convert these to image impact in arcseconds
699 # The 10.0 is microns / pixel
700 pixelScale = 0.2 # arcseconds / pixel - find this elsewhere?
701 camHexXY = camHexXY / 10.0 * pixelScale # arcseconds(image) / micron(hexapod)
702 camHexUV = camHexUV / 10.0 * pixelScale # arcseconds(image) / arcsecond(hexapod)
703 camCoefs = [camHexXY, camHexXY, 0, camHexUV, camHexUV, 0]
704 m2HexXY = m2HexXY / 10.0 * pixelScale # arcseconds(image) / micron(hexapod)
705 m2HexUV = m2HexUV / 10.0 * pixelScale # arcseconds(image) / arcsecond(hexapod)
706 m2Coefs = [m2HexXY, m2HexXY, 0, m2HexUV, m2HexUV, 0]
708 camHexMs = 0.0
709 for i in [0, 1, 3, 4]:
710 camhex = copy.deepcopy(mountData.camhexData[f"error{i}"])
711 camhex *= camCoefs[i]
712 camHexMs += np.mean(camhex * camhex)
713 camHexRms = np.sqrt(camHexMs) # in arcseconds image impact
715 m2HexMs = 0.0
716 for i in [0, 1, 3, 4]:
717 m2hex = copy.deepcopy(mountData.m2hexData[f"error{i}"])
718 m2hex *= m2Coefs[i]
719 m2HexMs += np.mean(m2hex * m2hex)
720 m2HexRms = np.sqrt(m2HexMs) # in arcseconds image impact
721 return (float(camHexRms), float(m2HexRms))