diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py index 47238988..2baae5bd 100644 --- a/scopesim/effects/spectral_efficiency.py +++ b/scopesim/effects/spectral_efficiency.py @@ -11,7 +11,7 @@ from .effects import Effect from .ter_curves import TERCurve -from ..utils import figure_factory, get_logger, check_keys +from ..utils import figure_factory, get_logger, check_keys, from_currsys from .data_container import DataContainer from ..optics import echelle @@ -149,7 +149,7 @@ def __init__(self, **kwargs): self.efficiency_generator = self._generate_efficiency_curve_func() def _generate_efficiency_curve_func(self) -> Callable: - trace_params = self.table + trace_params = from_currsys(self.table, self.cmds) spectrographs = {} for row in trace_params: prefix = row["prefix"] # note trance ids are assumed to be prefix_{order} @@ -173,13 +173,17 @@ def _generate_efficiency_curve_func(self) -> Callable: xdisp_groove_length = u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq'] echelle_groove_length = u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq'] pix_per_res_elem = row['fwhm'] + prescription = {} + if 'collimator1_focal_length' in trace_params.colnames: + for key in ('collimator1_focal_length', 'collimator2_focal_length', 'entrance_slit_scale'): + prescription[key] = row[key] * u.Unit(trace_params.meta[key + '_unit']) spectrographs[prefix] = echelle.spectrograph_factory(min_wave, max_wave, focal_len, design_res, echelle_angle, min_order, max_order, echelle_groove_length, pix_per_res_elem, disp_npix, xdisp_npix, pix_size, xdisp_groove_length=xdisp_groove_length, xdisp_beta_center=xdisp_beta_center, - dispersion_focal_len=dispersion_focal_len) + dispersion_focal_len=dispersion_focal_len, **prescription) self._spectrographs = spectrographs def efficiency_curve(trace_id, wavelength): @@ -197,9 +201,12 @@ def apply_to(self, obj, **kwargs): """Interface between FieldOfView and SpectralEfficiency.""" trace_id = obj.trace_id - swcs = WCS(obj.hdu.header).spectral - with u.set_enabled_equivalencies(u.spectral()): - wave = swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])) << u.um + if obj.hdu is None: + wave = obj.cube_waveset + else: + swcs = WCS(obj.hdu.header).spectral + with u.set_enabled_equivalencies(u.spectral()): + wave = swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])) << u.um efficiency = self.efficiency_generator(trace_id, wave) params = {"description": trace_id} @@ -207,8 +214,7 @@ def apply_to(self, obj, **kwargs): effic_curve = TERCurve(array_dict={"wavelength": wave, "transmission": efficiency}, **params) self.efficiencies[trace_id] = effic_curve - obj.hdu.data *= effic_curve.throughput(wave).value[:, None, None] - return obj + return effic_curve.apply_to(obj) def plot(self): """Plot the grating efficiencies.""" diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 6597ebb2..b10cbb49 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -678,6 +678,19 @@ class EchelleSpectralTraceList(SpectralTraceList): nIR 2 0 40 24 970 2500 20000 64.2 225 225 4.7 0.7 10.0 10 0.015 4096 4096 45 175 10 x ---------------------------------------------------------------- + Optional ``detector_offset_x`` and ``detector_offset_y`` columns specify + the rotated analytical camera origin relative to detector center, positive + right/up. Both columns and their ``*_unit`` metadata are required together; + values must have length units. Explicit zeros place that origin at detector + center. Offsets are applied after rotation, before detector clipping, with + no automatic recentering. Omitting both columns retains automatic placement. + + With ``physical_slit_dispersion=True``, supply ``collimator1_focal_length``, + ``collimator2_focal_length`` and ``entrance_slit_scale`` columns with their + ``*_unit`` metadata. The entrance scale is length per on-sky arcsecond. + These determine the illumination angle and asymmetric pupil relay; + the along-slit plate scale is not used for across-slit dispersion. + The calculated traces are stored in the same HDUList format as required by SpectralTraceList, and supplied to the parent class through hdulist kwarg. @@ -690,10 +703,13 @@ def __init__(self, cmds=None, **kwargs): self.cmds = cmds trace_param_filename = kwargs.pop("filename") + self.physical_slit_dispersion = kwargs.pop("physical_slit_dispersion", False) + self.slit_samples = kwargs.pop("slit_samples", 3) save_generated_hdulist = kwargs.pop("save_generated_hdulist", False) generated_hdulist_filename = kwargs.pop( "generated_hdulist_filename", "analytical_echelle_traces.fits") trace_params = DataContainer(filename=trace_param_filename) + trace_params.table = from_currsys(trace_params.table, self.cmds) hdulist = self._generate_trace_hdulist(trace_params) if save_generated_hdulist: output_path = Path(generated_hdulist_filename).expanduser() @@ -712,6 +728,7 @@ def _generate_trace_hdulist(self, trace_params): hdul[0].header["EDATA"] = 2 trace_ids, ap_ids, im_ids = [], [], [] + self._spectrographs = {} for row in trace_params.table: prefix = row["prefix"] min_order = row['m0'] - row['n'] + 1 @@ -740,13 +757,19 @@ def _generate_trace_hdulist(self, trace_params): xdisp_groove_length = u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq'] echelle_groove_length = u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq'] pix_per_res_elem = row['fwhm'] + prescription = {} + if 'collimator1_focal_length' in trace_params.table.colnames: + for key in ('collimator1_focal_length', 'collimator2_focal_length', 'entrance_slit_scale'): + prescription[key] = row[key] * u.Unit(trace_params.meta[key + '_unit']) + nominal_width = row['nominal_slit_width'] * u.Unit(trace_params.meta['nominal_slit_width_unit']) + pix_per_res_elem = (nominal_width * prescription['entrance_slit_scale'] * focal_len / prescription['collimator2_focal_length'] / pix_size).to_value(u.dimensionless_unscaled) ss = echelle.spectrograph_factory(min_wave, max_wave, focal_len, design_res, echelle_angle, min_order, max_order, echelle_groove_length, pix_per_res_elem, disp_npix, xdisp_npix, pix_size, xdisp_groove_length=xdisp_groove_length, xdisp_beta_center=xdisp_beta_center, - dispersion_focal_len=dispersion_focal_len) + dispersion_focal_len=dispersion_focal_len, **prescription) _warn_if_echelle_design_res_inconsistent( prefix, design_res, @@ -757,7 +780,7 @@ def _generate_trace_hdulist(self, trace_params): ) slit_edge = (row['slitlength'] / 2) * u.Unit(trace_params.meta["slitlength_unit"]) - slit_pos = np.linspace(-slit_edge, slit_edge, num=3) + slit_pos = np.linspace(-slit_edge, slit_edge, num=self.slit_samples) if "plate_scale" in trace_params.table.colnames: plate_scale = row["plate_scale"] * u.Unit( trace_params.meta["plate_scale_unit"]) @@ -776,74 +799,52 @@ def _generate_trace_hdulist(self, trace_params): ).to_value(u.deg) cang = np.cos(np.deg2rad(detector_angle)) sang = np.sin(np.deg2rad(detector_angle)) - detector_x_padding = max( - 0.0, - ( - disp_npix * abs(cang) - + xdisp_npix * abs(sang) - - disp_npix - ) / 2, - ) - slit_x_padding = ( - abs(np.tan(np.deg2rad(detector_angle))) - * float(np.nanmax(np.abs(slit_offset_pix))) - ) - raw_x_min = 0.5 - detector_x_padding - slit_x_padding - raw_x_max = disp_npix - 0.5 + detector_x_padding + slit_x_padding - - def raw_detector_pixels(wave, order): - x = ss.wavelength_to_x_pixel(wave, order) - y = ss.wavelength_to_y_pixel(wave) - x = u.Quantity(x, copy=False).to_value( - u.dimensionless_unscaled) - y = u.Quantity(y, copy=False).to_value( - u.dimensionless_unscaled) - xpix = np.broadcast_to(x, (slit_offset_pix.size, wave.size)) - ypix = y[None, :] + slit_offset_pix[:, None] - return xpix, ypix - - rotated_x_min = rotated_y_min = np.inf - rotated_x_max = rotated_y_max = -np.inf - # The analytical echelle layout is relative, not yet detector - # placed. Rotate the footprint first, then place that rotated - # footprint on the detector. Detector padding is not part of the - # optical trace geometry; downstream detector/FOV code is - # responsible for clipping padded display or extraction regions. - for i, order in enumerate(ss.orders): - raw_x = np.linspace( - raw_x_min, raw_x_max, num=max(disp_npix, 2)) - wave = ss.x_pixel_to_wavelength(raw_x, order) - xpix, ypix = raw_detector_pixels(wave, order) - xrot = cang * xpix - sang * ypix - yrot = sang * xpix + cang * ypix - rotated_x_min = min(rotated_x_min, float(np.nanmin(xrot))) - rotated_x_max = max(rotated_x_max, float(np.nanmax(xrot))) - rotated_y_min = min(rotated_y_min, float(np.nanmin(yrot))) - rotated_y_max = max(rotated_y_max, float(np.nanmax(yrot))) - rotated_x_center = rotated_x_min + (rotated_x_max - rotated_x_min) / 2 - rotated_y_center = rotated_y_min + (rotated_y_max - rotated_y_min) / 2 - - def rotated_detector_pixels(wave, order): - xpix, ypix = raw_detector_pixels(wave, order) - xrot = cang * xpix - sang * ypix - yrot = sang * xpix + cang * ypix - return ( - xrot - rotated_x_center + disp_npix / 2, - yrot - rotated_y_center + xdisp_npix / 2, - ) + ss._detector_cos, ss._detector_sin = cang, sang + ss._slit_offset_pix = slit_offset_pix + self._spectrographs[prefix] = ss + if 'detector_offset_x' in trace_params.table.colnames or 'detector_offset_y' in trace_params.table.colnames: + # Multiplying an ASCII table's masked cell by a Unit yields + # zero, so reject missing values before attaching units. + if np.ma.is_masked(row['detector_offset_x']) or np.ma.is_masked(row['detector_offset_y']): + raise ValueError('Both detector offsets must have values.') + offset_x = row['detector_offset_x'] * u.Unit(trace_params.meta['detector_offset_x_unit']) + offset_y = row['detector_offset_y'] * u.Unit(trace_params.meta['detector_offset_y_unit']) + # Offsets locate the rotated camera origin relative to detector + # center. The raw mapping includes the unrotated pixel origin. + rotated_x_center = cang*disp_npix/2 - sang*xdisp_npix/2 - (offset_x/pix_size).to_value(u.dimensionless_unscaled) + rotated_y_center = sang*disp_npix/2 + cang*xdisp_npix/2 - (offset_y/pix_size).to_value(u.dimensionless_unscaled) + + # Project the translated detector rectangle back onto model x; + # slit offsets are in model y and cancel in this projection. + raw_x_center = cang*rotated_x_center + sang*rotated_y_center + raw_x_halfspan = (abs(cang)*disp_npix + abs(sang)*xdisp_npix)/2 + candidate_raw_x = np.linspace(raw_x_center - raw_x_halfspan, raw_x_center + raw_x_halfspan, disp_npix + 1) + else: + # Preserve the original automatic placement and sampling for + # parameter files that do not specify a detector position. + detector_x_padding = max(0., (disp_npix*abs(cang) + xdisp_npix*abs(sang) - disp_npix)/2) + slit_x_padding = abs(np.tan(np.deg2rad(detector_angle))) * np.max(np.abs(slit_offset_pix)) + raw_x_min = .5 - detector_x_padding - slit_x_padding + raw_x_max = disp_npix - .5 + detector_x_padding + slit_x_padding + candidate_raw_x = np.linspace(raw_x_min, raw_x_max, num=max(disp_npix, 2)) + rotated_x_min = rotated_y_min = np.inf + rotated_x_max = rotated_y_max = -np.inf + for order in ss.orders: + wave = ss.x_pixel_to_wavelength(candidate_raw_x, order) + xrot, yrot = ss._rotated_detector_pixels(wave, order) + rotated_x_min = min(rotated_x_min, xrot.min()) + rotated_x_max = max(rotated_x_max, xrot.max()) + rotated_y_min = min(rotated_y_min, yrot.min()) + rotated_y_max = max(rotated_y_max, yrot.max()) + rotated_x_center = rotated_x_min + (rotated_x_max - rotated_x_min)/2 + rotated_y_center = rotated_y_min + (rotated_y_max - rotated_y_min)/2 + ss._detector_center = (rotated_x_center, rotated_y_center) def on_detector(wave, order): - xpix, ypix = rotated_detector_pixels(wave, order) - return ( - (xpix >= 0) - & (xpix <= disp_npix) - & (ypix >= 0) - & (ypix <= xdisp_npix) - ) + xpix, ypix = ss._detector_pixels(wave, order) + return (xpix >= 0) & (xpix <= disp_npix) & (ypix >= 0) & (ypix <= xdisp_npix) - for i, order in enumerate(ss.orders): - candidate_raw_x = np.linspace( - raw_x_min, raw_x_max, num=max(disp_npix, 2)) + for order in ss.orders: candidate_wave = ss.x_pixel_to_wavelength(candidate_raw_x, order) valid_wave = np.any( on_detector(candidate_wave, order), axis=0) @@ -872,7 +873,7 @@ def on_detector(wave, order): continue s = np.tile(slit_pos, wave.size).reshape(wave.size, slit_pos.size).T.ravel() w = np.tile(wave, slit_offset_pix.size) - xpix, ypix = rotated_detector_pixels(wave, order) + xpix, ypix = ss._detector_pixels(wave, order) xval = (xpix.ravel() - disp_npix / 2) * pix_size.to(u.mm) yval = (ypix.ravel() - xdisp_npix / 2) * pix_size.to(u.mm) @@ -881,7 +882,12 @@ def on_detector(wave, order): 'x': xval, 'y': yval}) + if self.physical_slit_dispersion: + order_table['dlam_deta'] = ss.slit_dispersion(w, order) + trace_hdu = fits.BinTableHDU(order_table) + if self.physical_slit_dispersion: + trace_hdu.header['SLITCOL'] = 'dlam_deta' trace_hdu.header['DISPDIR'] = row['dispdir'] trace_hdu.header["EXTNAME"] = f'{prefix}_{order:d}' trace_hdu.header["DESIGNR"] = ( diff --git a/scopesim/optics/echelle.py b/scopesim/optics/echelle.py index 4c7614e6..709b2fd2 100644 --- a/scopesim/optics/echelle.py +++ b/scopesim/optics/echelle.py @@ -11,7 +11,10 @@ def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, echelle_groove_length: float|u.Quantity, pix_per_res_elem: float, disp_npix: int, xdisp_npix: int, pix_size: float|u.Quantity, xdisp_groove_length: float|u.Quantity = 0.0, xdisp_beta_center: float|u.Quantity = 0.0, - dispersion_focal_len: float|u.Quantity = None): + dispersion_focal_len: float|u.Quantity = None, + collimator1_focal_length: u.Quantity = None, + collimator2_focal_length: u.Quantity = None, + entrance_slit_scale: u.Quantity = None): """ Parameters @@ -47,6 +50,11 @@ def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, dispersion_focal_len: float|u.Quantity Effective focal length for the echelle/main-dispersion detector coordinate. If float, assumed to be in mm. Defaults to focal_len. + collimator1_focal_length, collimator2_focal_length: u.Quantity + Physical collimator lengths. When supplied together, derive the + dispersion focal length from focal_len times their ratio. + entrance_slit_scale: u.Quantity + Entrance-slit length per on-sky angle, required with the collimators. Returns ------- @@ -98,6 +106,9 @@ def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, pixels_per_res_elem=pix_per_res_elem, focal_length=focal_len, dispersion_focal_length=dispersion_focal_len, + collimator1_focal_length=collimator1_focal_length, + collimator2_focal_length=collimator2_focal_length, + entrance_slit_scale=entrance_slit_scale, grating=echelle_grating, detector=Detector(disp_npix, xdisp_npix, pix_size), cross_disperser=cross_disperser) @@ -191,16 +202,21 @@ def estimate_xdisp_angle_with_groove_len( if not isinstance(d, u.Quantity): d = d * u.mm - # assume we want a dispersion at the midpoint that matches the range to the detector length - t = detector_length / focal_length - k = (l_max - l_min) / d - x = k * (1 + np.sqrt(1 + t ** 2)) / 2 / t - beta_center_guess = np.arccos(np.sqrt(x)) - beta_half_range = np.arcsin(k / 2 / np.sqrt(x)) - beta_max = beta_center_guess + beta_half_range - beta_min = beta_center_guess - beta_half_range - alpha = np.arcsin(l_max / d - np.sin(beta_max)) - return alpha.to(u.rad), ((beta_min + beta_max) / 2).to(u.rad) + # The camera maps beta_center +/- half_angle to +/- detector_length/2. + # Subtracting the two grating equations eliminates the incidence angle: + # (l_max - l_min)/d = 2*cos(beta_center)*sin(half_angle). + if l_max <= l_min or detector_length <= 0 * u.mm or focal_length <= 0 * u.mm or d <= 0 * u.mm: + raise ValueError("Wavelength range, detector length, focal length and groove spacing must be positive.") + half_angle = np.arctan(detector_length / (2 * focal_length)) + cos_center = (l_max - l_min) / (2 * d * np.sin(half_angle)) + if not 0 < cos_center <= 1: + raise ValueError("The wavelength range cannot fit this detector with the specified groove spacing.") + beta_center = np.arccos(cos_center) + sin_alpha = l_min / d - np.sin(beta_center - half_angle) + if beta_center + half_angle >= 90 * u.deg or not -1 <= sin_alpha <= 1: + raise ValueError("The requested detector span requires a non-propagating grating angle.") + alpha = np.arcsin(sin_alpha) + return alpha.to(u.rad), beta_center.to(u.rad) @staticmethod def estimate_xdisp_freq_and_angle(l_min: float | u.Quantity, @@ -364,6 +380,9 @@ def __init__( detector: Detector, cross_disperser: GratingSetup = None, dispersion_focal_length: u.Quantity = None, + collimator1_focal_length: u.Quantity = None, + collimator2_focal_length: u.Quantity = None, + entrance_slit_scale: u.Quantity = None, ): """ @@ -398,6 +417,12 @@ def __init__( if dispersion_focal_length is None else dispersion_focal_length.to(u.mm) ) + if collimator1_focal_length is not None: + self.collimator1_focal_length = collimator1_focal_length.to(u.mm) + self.collimator2_focal_length = collimator2_focal_length.to(u.mm) + self.entrance_slit_scale = entrance_slit_scale.to(u.mm / u.arcsec) + # Angular dispersion is amplified by the asymmetric pupil relay. + self.dispersion_focal_length = self.focal_length * self.collimator1_focal_length / self.collimator2_focal_length self.pixel_scale = np.arctan( self.detector.pixel_size / self.dispersion_focal_length) self.beta_central_pixel = self.grating.beta_center @@ -432,6 +457,22 @@ def set_beta_center(self, beta, littrow: bool = False): if littrow: self.grating.alpha = beta + def slit_dispersion(self, wavelength, order): + """Apparent wavelength shift per across-slit arcsecond at fixed incident wavelength. + + m*lambda = d*(sin(alpha) + sin(beta)); increasing entrance-slit eta + increases alpha and reduces the wavelength inferred at nominal alpha. + The first collimator sets this angle; the pupil relay is already in + the camera's dispersion focal length and must not enter a second time. + """ + derivative = -self.grating.d * np.cos(self.grating.alpha) * self.entrance_slit_scale / (order * self.collimator1_focal_length) + return np.ones(np.shape(wavelength)) * derivative.to(u.um / u.arcsec) + + def projected_slit_width(self, wavelength, order, slit_width): + """Unrotated echelle-coordinate slit width, including local anamorphism.""" + half_shift = self.slit_dispersion(wavelength, order) * slit_width / 2 + return abs(self.wavelength_to_x_pixel(wavelength + half_shift, order) - self.wavelength_to_x_pixel(wavelength - half_shift, order)) * u.pixel + @property def orders(self): try: @@ -558,17 +599,15 @@ def order_mask(self, wave, fsr_edge: bool = False): def edge_wave(self, fsr=True): """ :param fsr: True to return the FSR edge wavelengths, False to return the detector edges - :return: wavelengths at the edges of the detector (or FSR) for each order, always include central pixel wavelength + :return: the two endpoint wavelengths for each order """ - pix = np.array([0, self.detector.n_pix_x // 2, self.detector.n_pix_x-1]) + .5 - fiducial_waves = self.x_pixel_to_wavelength(pix, self.orders[:, None]) + pixels = np.array([0, self.detector.n_pix_x//2, self.detector.n_pix_x - 1]) + .5 + fiducial_waves = self.x_pixel_to_wavelength(pixels, self.orders[:, None]) if not fsr: return fiducial_waves[:, [0, -1]] - + # Retain the original central-pixel wavelength convention for the FSR. central_fsr = fiducial_waves[:, 1] / self.orders - fsr_edges = (u.Quantity([-central_fsr / 2, central_fsr / 2]) + fiducial_waves[:, 1]).T - - return fsr_edges + return (u.Quantity([-central_fsr/2, central_fsr/2]) + fiducial_waves[:, 1]).T def central_wave(self, order): """ @@ -616,6 +655,22 @@ def x_pixel_wavelengths(self, pixel_phase=.5): def wavelength_to_y_pixel(self, wave): return self.y_pixel_for_beta(self.cross_disperser.beta(wave, 1)) + def _rotated_detector_pixels(self, wave, order): + x = self.wavelength_to_x_pixel(wave, order).to_value(u.dimensionless_unscaled) + y = self.wavelength_to_y_pixel(wave).to_value(u.dimensionless_unscaled) + xpix = np.broadcast_to(x, (self._slit_offset_pix.size, wave.size)) + ypix = y[None, :] + self._slit_offset_pix[:, None] + return self._detector_cos*xpix - self._detector_sin*ypix, self._detector_sin*xpix + self._detector_cos*ypix + + def _detector_pixels(self, wave, order): + """The trace generator's placed coordinates for all configured slit nodes. + + Unlike the clipped trace table, this mapping can also inspect the full + FSR outside the detector. It uses the same rotation and placement. + """ + x, y = self._rotated_detector_pixels(wave, order) + return x - self._detector_center[0] + self.detector.n_pix_x/2, y - self._detector_center[1] + self.detector.n_pix_y/2 + def y_pixel_to_wavelength(self, pixel, m): """ :param pixel: pixel diff --git a/scopesim/tests/tests_optics/test_echelle.py b/scopesim/tests/tests_optics/test_echelle.py new file mode 100644 index 00000000..2370f03d --- /dev/null +++ b/scopesim/tests/tests_optics/test_echelle.py @@ -0,0 +1,102 @@ +"""Planar grating geometry checked through the camera's forward projection.""" + +import numpy as np +import pytest +from astropy import units as u + +from scopesim.optics.echelle import GratingSetup, spectrograph_factory + + +# Wavelength endpoints [nm], cross-disperser density [mm^-1], detector height [pix]. +CHANNELS = [(310, 430, 2000, 4096), (430, 650, 1100, 4096), (650, 980, 625, 4096), + (960, 1340, 350, 2048), (1470, 1830, 300, 2048), (1990, 2450, 270, 2048)] + + +def test_asymmetric_pupil_slit_dispersion_and_anamorphism(): + instrument = spectrograph_factory(960, 1340, 225, 18000, 65.6, 58, 80, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=1/350, xdisp_beta_center=10.9, + collimator1_focal_length=1500*u.mm, collimator2_focal_length=1250*u.mm, entrance_slit_scale=.5306*u.mm/u.arcsec) + assert instrument.dispersion_focal_length == 270*u.mm + center = instrument.central_wave(65) + derivative = instrument.slit_dispersion(center, 65) + assert derivative < 0*u.um/u.arcsec + for slit, expected in ((.7, 17800), (.33, 37800)): + resolution = (center / abs(derivative * slit*u.arcsec)).to_value(u.dimensionless_unscaled) + assert resolution == pytest.approx(expected, rel=.01) + width = instrument.projected_slit_width(center, 65, .7*u.arcsec) + assert width.to_value(u.pixel) == pytest.approx(4.45704, rel=1e-5) + wavelengths = center * np.array([.99, 1., 1.01]) + widths = instrument.projected_slit_width(wavelengths, 65, .7*u.arcsec) + assert np.all(np.diff(widths) > 0*u.pixel) + half_angle = (.35*u.arcsec * instrument.entrance_slit_scale / instrument.collimator1_focal_length).value * u.rad + beta_low = np.arcsin(65*wavelengths/instrument.grating.d - np.sin(instrument.grating.alpha + half_angle)) + beta_high = np.arcsin(65*wavelengths/instrument.grating.d - np.sin(instrument.grating.alpha - half_angle)) + edge_width = abs(instrument.x_pixel_for_beta(beta_high) - instrument.x_pixel_for_beta(beta_low)) + # Independent finite slit edges retain the same anamorphism as dispersion. + np.testing.assert_allclose(widths.to_value(u.pixel), edge_width.value, rtol=2e-7) + step = center*1e-7 + dispersion = (instrument.wavelength_to_x_pixel(wavelengths+step, 65)-instrument.wavelength_to_x_pixel(wavelengths-step, 65))/(2*step) + measured_delta = widths.to_value(u.pixel)/dispersion + np.testing.assert_allclose(measured_delta.to_value(u.nm), abs(derivative*.7*u.arcsec).to_value(u.nm), rtol=2e-7) + eta = .0825*u.arcsec + dalpha = (eta * instrument.entrance_slit_scale / instrument.collimator1_focal_length).to_value(u.dimensionless_unscaled) * u.rad + exact = instrument.grating.d / 65 * (np.sin(instrument.grating.alpha) - np.sin(instrument.grating.alpha + dalpha)) + assert (derivative * eta).to_value(u.nm) == pytest.approx(exact.to_value(u.nm), rel=1e-4) + + +def test_incomplete_physical_prescription_fails(): + with pytest.raises(AttributeError): + spectrograph_factory(960, 1340, 225, 18000, 65.6, 58, 80, 1/24, 4.46, 2048, 2048, .015, + xdisp_groove_length=1/350, collimator1_focal_length=1500*u.mm) + + +@pytest.mark.parametrize("low,high,density,height", CHANNELS) +def test_fixed_groove_endpoints_fill_camera(low, high, density, height): + instrument = spectrograph_factory(low, high, 225, 18000, 65.6, 40, 50, 1/24, 4.46, + 2048, height, .015, xdisp_groove_length=1/density, dispersion_focal_len=270) + wavelength = np.array([low, high]) * u.nm + pixels = instrument.wavelength_to_y_pixel(wavelength) + # Algebraic inverse/forward evaluation; 1e-7 pixel allows floating-point cancellation. + np.testing.assert_allclose(pixels.value, [0, height], atol=1e-7, rtol=0) + recovered = instrument.y_pixel_to_wavelength(pixels, 1) + np.testing.assert_allclose(recovered.to_value(u.nm), wavelength.value, atol=1e-9, rtol=0) + + +def test_units_and_signed_incidence(): + numeric = GratingSetup.estimate_xdisp_angle_with_groove_len(1470, 1830, 30.72, 225, 1/300) + physical = GratingSetup.estimate_xdisp_angle_with_groove_len(1.47*u.um, 1.83*u.um, 3.072*u.cm, .225*u.m, (1/300)*u.mm) + assert physical[0] < 0 * u.rad + for plain, quantity in zip(numeric, physical): + assert plain.unit == u.rad + assert quantity.unit == u.rad + assert plain.to_value(u.rad) == pytest.approx(quantity.to_value(u.rad), abs=1e-12) + + +@pytest.mark.parametrize("inputs", [(310, 430, 1, 225, .0005), (3000, 3010, 61.44, 225, .0001), + (430, 310, 61.44, 225, .0005), (310, 430, 0, 225, .0005)]) +def test_impossible_geometry_fails(inputs): + with pytest.raises(ValueError): + GratingSetup.estimate_xdisp_angle_with_groove_len(*inputs) + + +def test_explicit_angle_keeps_configured_geometry(): + instrument = spectrograph_factory(960, 1340, 225, 18000, 65.6, 57, 80, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=1/350, xdisp_beta_center=12) + assert instrument.cross_disperser.alpha.to_value(u.deg) == pytest.approx(12) + assert instrument.cross_disperser.beta_center.to_value(u.deg) == pytest.approx(12) + + +def test_free_frequency_endpoints_fill_camera(): + instrument = spectrograph_factory(960, 1340, 225, 18000, 65.6, 57, 80, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=np.nan) + pixels = instrument.wavelength_to_y_pixel(np.array([960, 1340]) * u.nm) + np.testing.assert_allclose(pixels.value, [0, 2048], atol=1e-7, rtol=0) + + +def test_fsr_retains_the_central_pixel_wavelength_convention(): + instrument = spectrograph_factory(1990, 2450, 225, 18000, 65.6, 31, 38, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=1/270, dispersion_focal_len=270) + edges = instrument.edge_wave(fsr=True) + center = instrument.x_pixel_to_wavelength(instrument.detector.n_pix_x//2 + .5, instrument.orders) + np.testing.assert_allclose(edges.mean(axis=1).to_value(u.nm), center.to_value(u.nm), rtol=1e-12) + np.testing.assert_allclose(np.diff(edges, axis=1).ravel().to_value(u.nm), (center/instrument.orders).to_value(u.nm), rtol=1e-12)