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