Wavelength Mapping

This notebook contains illustrations of how various wavelengths are mapped to different pixels within the detector.

import numpy as np
import astropy.units as u
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
from sunpy.coordinates import Helioprojective, get_horizons_coord
from astropy.wcs.utils import pixel_to_pixel

from mocksipipeline.detector.response import SpectrogramChannel
from overlappy.wcs import overlappogram_fits_wcs, pcij_matrix

line_ids = [
    #('Fe XVIII',14.21*u.angstrom),  # also targeted by MaGIXS
    ('Fe XVII', 15.01*u.angstrom),  # also targeted by MaGIXS
    ('Fe XVII', 16.78*u.AA),
    ('Fe XVII', 17.05*u.AA),
    ('O VII', 21.60*u.angstrom),  # also targeted by MaGIXS
    ('O VII', 21.81*u.angstrom),
    ('O VII', 22.10*u.AA),
    ('O VIII', 18.97*u.angstrom),  # also targeted by MaGIXS
    ('Fe XXV', 1.86*u.AA),
    ('Ca XIX', 3.21*u.AA),
    ('Si XIII', 6.74*u.AA),
    ('Mg XI', 9.32*u.AA),
    ('Fe XVII', 11.25*u.AA),
    ('Fe XX', 12.83*u.AA),
    ('Ne IX', 13.45*u.AA),
    ('Fe XIX', 13.53*u.AA),
    ('C VI', 33.73*u.AA),
    ('C V', 40.27*u.AA),
    ('Si XII', 44.16*u.AA),
    ('Si XI', 49.18*u.AA),
]

def add_line_ids_to_axis(ax, line_ids):
    for ion,line in line_ids:
        ax.axvline(x=line.to_value('Angstrom'), ls=':', color='k',)
    ax2 = ax.secondary_xaxis('top')
    ax2.set_xticks(u.Quantity([l for _,l in line_ids]).to_value('Angstrom'),
                   labels=[f'{ion}, {line.to_string(format="latex_inline")}' for ion,line in line_ids],
                   rotation=90,
                   horizontalalignment='center',
                   #verticalalignment='center'
                  )

First, we are going to define a WCS for the 0th order and the 1st order

chan_0 = SpectrogramChannel(0, [])
chan_1 = SpectrogramChannel(1, [])
chan_3 = SpectrogramChannel(3, [])
chan_m1 = SpectrogramChannel(-1, [])
chan_m3 = SpectrogramChannel(-3, [])

observer = get_horizons_coord('SDO', '2020-11-09 18:00:00')
hpc_frame = Helioprojective(observer=observer, obstime=observer.obstime)

INFO: Obtained JPL HORIZONS location for Solar Dynamics Observatory (spac [sunpy.coordinates.ephemeris]

roll_angle = -90 * u.deg
dispersion_angle = 0 * u.deg

wcs0 = overlappogram_fits_wcs(
    chan_0.detector_shape,
    chan_0.wavelength,
    (chan_0.resolution[0], chan_0.resolution[1], chan_0.spectral_resolution),
    reference_pixel=chan_0.reference_pixel,
    reference_coord=(0*u.arcsec, 0*u.arcsec,0*u.angstrom),
    pc_matrix=pcij_matrix(roll_angle, dispersion_angle, order=chan_0.spectral_order, dispersion_axis=0),
    observer=observer,
)
wcs1 = overlappogram_fits_wcs(
    chan_1.detector_shape,
    chan_1.wavelength,
    (chan_1.resolution[0], chan_1.resolution[1], chan_1.spectral_resolution),
    reference_pixel=chan_1.reference_pixel,
    reference_coord=(0*u.arcsec, 0*u.arcsec,0*u.angstrom),
    pc_matrix=pcij_matrix(roll_angle, dispersion_angle, order=chan_1.spectral_order, dispersion_axis=0),
    observer=observer,
)
wcs3 = overlappogram_fits_wcs(
    chan_3.detector_shape,
    chan_3.wavelength,
    (chan_3.resolution[0], chan_3.resolution[1], chan_3.spectral_resolution),
    reference_pixel=chan_3.reference_pixel,
    reference_coord=(0*u.arcsec, 0*u.arcsec,0*u.angstrom),
    pc_matrix=pcij_matrix(roll_angle, dispersion_angle, order=chan_3.spectral_order, dispersion_axis=0),
    observer=observer,
)
wcsm1 = overlappogram_fits_wcs(
    chan_m1.detector_shape,
    chan_m1.wavelength,
    (chan_m1.resolution[0], chan_m1.resolution[1], chan_m1.spectral_resolution),
    reference_pixel=chan_m1.reference_pixel,
    reference_coord=(0*u.arcsec, 0*u.arcsec, 0*u.angstrom),
    pc_matrix=pcij_matrix(roll_angle, dispersion_angle, order=chan_m1.spectral_order, dispersion_axis=0),
    observer=observer,
)
wcsm3 = overlappogram_fits_wcs(
    chan_m3.detector_shape,
    chan_m3.wavelength,
    (chan_m3.resolution[0], chan_m3.resolution[1], chan_m3.spectral_resolution),
    reference_pixel=chan_m3.reference_pixel,
    reference_coord=(0*u.arcsec, 0*u.arcsec, 0*u.angstrom),
    pc_matrix=pcij_matrix(roll_angle, dispersion_angle, order=chan_m3.spectral_order, dispersion_axis=0),
    observer=observer,
)

#ar_coord = SkyCoord(Tx=[100]*u.arcsec, Ty=[-450]*u.arcsec, frame=hpc_frame)
ar_coord = SkyCoord(Tx=[-1000]*u.arcsec, Ty=[0]*u.arcsec, frame=hpc_frame)
#ar_coord = SkyCoord(Tx=[0]*u.arcsec, Ty=[-1200]*u.arcsec, frame=hpc_frame)
ar_coord_x_pix, ar_coord_y_pix, _ = ar_coord_pixel = wcs0.world_to_pixel(ar_coord, 0*u.angstrom)

wave_indices = range(1300)
center_x_o0, _, _ = pixel_to_pixel(wcs0, wcs0, ar_coord_x_pix[0], ar_coord_y_pix[0], wave_indices)
center_x_o1, _, _ = pixel_to_pixel(wcs0, wcs1, ar_coord_x_pix[0], ar_coord_y_pix[0], wave_indices)
center_x_o3, _, _ = pixel_to_pixel(wcs0, wcs3, ar_coord_x_pix[0], ar_coord_y_pix[0], wave_indices)
center_x_om1, _, _ = pixel_to_pixel(wcs0, wcsm1, ar_coord_x_pix[0], ar_coord_y_pix[0], wave_indices)
center_x_om3, _, _ = pixel_to_pixel(wcs0, wcsm3, ar_coord_x_pix[0], ar_coord_y_pix[0], wave_indices)

_,wavelength = wcs0.pixel_to_world(*(chan_0.reference_pixel[:2]-1*u.pix).value, wave_indices)

fig = plt.figure(figsize=(18,5))
ax = fig.add_subplot()
ax.plot(wavelength.to('Angstrom'),center_x_o0, label='order=0')
ax.plot(wavelength.to('Angstrom'),center_x_o1, label='order=1')
ax.plot(wavelength.to('Angstrom'),center_x_o3, label='order=3')
ax.plot(wavelength.to('Angstrom'),center_x_om1, label='order=-1')
ax.plot(wavelength.to('Angstrom'),center_x_om3, label='order=-3')
add_line_ids_to_axis(ax, line_ids)
for ion,line in line_ids:
    #if ion in ('Fe XVII', 'O VII', 'O VIII'):
    pix_pos,_,_ = wcs1.world_to_pixel(ar_coord, line)
    ax.hlines(pix_pos, 0, line.to_value('Angstrom'), ls=':', color='C1')
    pix_pos,_,_ = wcs3.world_to_pixel(ar_coord, line)
    ax.hlines(pix_pos, 0, line.to_value('Angstrom'), ls=':', color='C2')
    pix_pos,_,_ = wcsm1.world_to_pixel(ar_coord, line)
    ax.hlines(pix_pos, 0, line.to_value('Angstrom'), ls=':', color='C3')
    pix_pos,_,_ = wcsm3.world_to_pixel(ar_coord, line)
    ax.hlines(pix_pos, 0, line.to_value('Angstrom'), ls=':', color='C4')
ax.axhline(y=0, ls='-', color='k')
ax.axhline(y=2000, ls='-', color='k')
ax.axvline(x=0, ls='-', color='k')
ax.set_xlim(0,70)
ax.set_ylim(1500,2500)
ax.legend(loc=4)
ax.set_ylabel('Position of Disk Center on detector')
ax.set_xlabel('Wavelength [Å]')

Text(0.5, 0, 'Wavelength [Å]')