Simulating and Inverting a CTIS Observation

position_scene = na.Cartesian2dVectorLinearSpace(
    start=-10 * u.arcsec,
    stop=10 * u.arcsec,
    axis=na.Cartesian2dVectorArray("scene_x", "scene_y"),
    num=na.Cartesian2dVectorArray(64 + 1, 64 + 1),
)

position_sensor = na.Cartesian2dVectorArray(
    x=na.arange(0, 128 + 1, axis="sensor_x") * u.pix,
    y=na.arange(0, 64 + 1, axis="sensor_y") * u.pix,
)
        

position_scene.x.shape
{'scene_x': 65}
position_scene.shape
{'scene_x': 65, 'scene_y': 65}
position_sensor.shape
{'sensor_x': 129, 'sensor_y': 65}

The scene has 64 × 64 cells across 20″, a 0.3125″ pitch; the plate scale in step 8 is 0.4″ per pixel.

coordinates_scene = na.DopplerPositionalVectorArray.from_velocity(
    velocity=velocity,
    wavelength_rest=wavelength_rest,
    position=position_scene,
)

coordinates_sensor = na.DopplerPositionalVectorArray.from_velocity(
    velocity=velocity,
    wavelength_rest=wavelength_rest,
    position=position_sensor,
)
        

coordinates_scene.shape
{'wavelength': 21, 'scene_x': 65, 'scene_y': 65}
coordinates_sensor.shape
{'wavelength': 21, 'sensor_x': 129, 'sensor_y': 65}

colorbar = na.plt.rgbmesh(
    C=scene,
    axis_wavelength="wavelength",
    ax=ax,
    vmin=0,
    vmax=scene.outputs.max(),
)
na.plt.pcolormesh(C=colorbar, axis_rgb="wavelength", ax=cax)
        

• Collapses the "wavelength" axis into color with colorsynth: each pixel's spectrum is stretched onto the visible range and weighted by the CIE 1931 color-matching functions.
• Hue shows the shape of the spectrum, here the Doppler shift; brightness shows the radiance between vmin and vmax.
• Returns a 2D colorbar, color against radiance and wavelength, drawn with na.plt.pcolormesh(..., axis_rgb=...).

spectrum = scene.outputs.mean(("scene_x", "scene_y"))

ax2 = ax.twiny()
na.plt.stairs(velocity, spectrum, ax=ax)
na.plt.stairs(scene.inputs.wavelength, spectrum, ax=ax2)
        

• Reduce over the spatial axes by name.
• stairs takes the 21 vertices and the 20 bin values directly.
• Two components, about ±125 km/s from rest, on the 1% background.

instrument = ctis.instruments.IdealInstrument(
    area_effective=1 * u.cm ** 2,
    timedelta_exposure=20 * u.s,
    plate_scale=.4 * u.arcsec / u.pix,
    dispersion=dispersion,
    angle=angle,
    wavelength_ref=wavelength_rest,
    position_ref=na.Cartesian2dVectorArray(64, 32) * u.pix,
    coordinates_scene=coordinates_scene,
    coordinates_sensor=coordinates_sensor,
    channel="dispersion angle = " + angle.to_string_array("%03d"),
    axis_channel="channel",
    axis_wavelength="wavelength",
    axis_scene_xy=("scene_x", "scene_y"),
    axis_sensor_xy=("sensor_x", "sensor_y"),
)
        

solution = inversion.solution
spectrum_inverted = solution.outputs.mean(
    ("scene_x", "scene_y")
)

na.plt.stairs(scene.inputs.wavelength, spectrum,
              ax=ax, label="original")
na.plt.stairs(scene.inputs.wavelength,
              spectrum_inverted,
              ax=ax, label="reconstructed")
        

• The peaks come back 17% and 14% low.
• The dip between the components comes back 65% high.
• The missing flux goes into the line wings.

na.plt.plot(inversion.iteration,
            inversion.mean_chi_squared,
            ax=ax[0],
            axis=inversion.inverter.axis_iteration,
            label=instrument.channel)
na.plt.plot(inversion.iteration,
            inversion.correlation_residual,
            ax=ax[1],
            axis=inversion.inverter.axis_iteration,
            label=instrument.channel)
        

• $\langle \chi^2 \rangle$ per channel uses σ of the predicted image; the stopping rule uses its mean over channels.
• The residual's correlation with the predicted image tracks structure left in the residual.
• $\langle \chi^2 \rangle$ ends at 0.33: below 1, so the later iterations fit noise.