The purpose of this notebook is to create a spectral cube that can be used as a Source for METIS IFU (a.k.a. LMS) simulations with Scopesim. The goal is not to create a realistic astrophysical source nor to probe any signal-to-noise limits, but to demonstrate the formal requirements on a FITS file that can be used by Scopesim, and maybe hint at some of the capabilities of the METIS IFU.

import numpy as np
from matplotlib import pyplot as plt
from astropy.io import fits
from astropy import units as u
from astropy.wcs import WCS

# A function to set background colour of individual cells
from IPython.display import HTML, display

def set_background(color):    
    script = (
        "var cell = this.closest('.jp-CodeCell');"
        "var editor = cell.querySelector('.jp-Editor');"
        "editor.style.background='{}';"
        "this.parentNode.removeChild(this)"
    ).format(color)
    
    display(HTML('<img src onerror="{}" style="display:none">'.format(script)))

The idea is to use an (ALMA) image of HL Tau and turn it into a spectral cube including Doppler shifts due to disk rotation. The image is taken for our server.

import scopesim as sim
path = sim.download_example_data("HL_Tau_prep_for_Scopesim.fits")

hdul = fits.open(path[0])
hdul.info()

Any fits file that is to be used as source for Scopesim needs to have BUNIT keyword, which specifies the units of the data values, as well as the WCS keywords, from which Scopesim takes the pixel scale.

print("Data unit (BUNIT):", hdul[1].header["BUNIT"]) 

pixscale = np.abs(hdul[1].header['CDELT1'] * u.Unit(hdul[1].header['CUNIT1']))
print(f"Pixel scale: {pixscale.to(u.mas)}")

HL Tau is at a distance of 140 pc, but we will place it at 50 pc, which gives a scale of 20 mas per AU:

au_mas = 1000 * u.mas / 50  # 1 AU is 1 arcsec at 1 pc
print(f"1 AU at 50 pc: {au_mas}")

We shall model the rotation of the disk as Keplerian with a reference rotational velocity of 30 km/s at 1 AU, which gives a velocity shift of $$\frac{\Delta\lambda}{\lambda} = \frac{v}{c} = 10^{-4} \left(\frac{R_\mathrm{ell}}{20,\mathrm{mas}}\right)^{-1},$$ where $R_\mathrm{ell}$ is the projected radius in the disk onto the plane of the sky. The disk appears elliptical with $e = b/a\approx 0.75$ at a position angle of $\phi\approx 45^\circ$, so $$R_\mathrm{ell} = \sqrt{ x’^2 + (y’/e)^2}$$ with $$x’ = x \cos\phi - y \sin\phi$$ $$y’ = x\sin\phi + y \cos\phi$$

To implement this, we first cut the image to a size that is a little larger than the field of view of the image slicer of the METIS IFU, which is roughly $0.9,\mathrm{arcsec}\times 0.6,\mathrm{arcsec}$. The following numbers have been determined by inspection of the original image such that the disk is centred in the new image.

npix = int(0.9 / pixscale.to(u.arcsec).value)   # number of pixels needed to cover slice length (0.9 arcsec)

# cut image around disk centre
x1 = int(802.5 - npix / 2)
x2 = int(802.5 + npix / 2) + 1
y1 = int(839.5 - npix / 2)
y2 = int(839.5 + npix / 2) + 1
img_cut = hdul[1].data[y1:y2, x1:x2]

The WCS has to be adapted to the new image size by moving the reference point CRPIXi. For simplicity we change the projection from the original sine to a simple linear projection. With that we can proceed to updating the HDU list with the new data section and new WCS.

wcs = WCS(hdul[1].header)
wcs.wcs.crpix = [89.5, 89.5]
wcs.wcs.crval = [0., 0.]
wcs.wcs.ctype = ["LINEAR", "LINEAR"]

Note: The following cell is a temporary workaround to avoid a bug in Scopesim v0.11.2 and earlier. The bug will be fixed in the next release and cdelt can then be negative.

set_background("lightblue")
wcs.wcs.cdelt[0] = np.abs(wcs.wcs.cdelt[0])
wcs.wcs.cdelt[1] = np.abs(wcs.wcs.cdelt[1])

hdul[1].data = img_cut
hdul[1].header.update(wcs.to_header())

hdul.writeto("HL_Tau_img_cut.fits", overwrite=True)

Next, we set the velocity field derived above. We start by creating numpy arrays for the coordinates $x$ and $y$ and then transform them to $x’$ and $y’$, $R_\mathrm{ell}$ and finally to $\Delta \lambda/\lambda$ (i.e. dloglambda).

# Derive two images for the world coordinates from that
ny, nx = img_cut.shape
i, j = np.meshgrid(np.arange(nx), np.arange(ny))
x, y = wcs.all_pix2world(i, j, 0)
x = (x * u.Unit(wcs.wcs.cunit[0]) << u.mas).value 
y = (y * u.Unit(wcs.wcs.cunit[1]) << u.mas).value

phi = -45 * np.pi / 180
xx = x * np.cos(phi) - y * np.sin(phi)
yy = x * np.sin(phi) + y * np.cos(phi)

inc = 45 * np.pi / 180
ell = np.cos(inc)
R_ell = np.sqrt(xx**2 + (yy/ell)**2)

plt.imshow(hdul[1].data, origin='lower', norm='log')
plt.contour(R_ell)
plt.title(r"Contours of $R_{\mathrm{ell}}$ overlaid on HL Tau");

v_los = 30 * xx * au_mas.value / R_ell**2 * np.sin(inc)
dloglambda = v_los / 2.998e6

from matplotlib.colors import SymLogNorm
from matplotlib.cm import coolwarm
plt.imshow(v_los, norm=SymLogNorm(linthresh=13), cmap=coolwarm, origin='lower')
plt.colorbar(label="line-of-sight velocity [km/s]")
plt.title("Line-of-sight velocity field");

Spectral dimension

We choose a central wavelength of 3.5$,\mu$m. The LMS covers a wavelength range from 3.472 to 3.526$,\mu$m (these values can be obtained from metis['lms_spectral_traces'].meta['wave_min'] and metis['lms_spectral_traces'].meta['wave_max'] when metis is set up as an OpticalTrain at the given central wavelength). A data cube used as Source for a Scopesim simulation should cover this wavelength range (a smaller range should work but fills the detector images only partially) at a wavelength sampling that is comparable to or better than the dispersion (in microns per pixel) on the LMS detectors.An indication is given by the scopesim parameter !SIM.spectral.spectral_bin_width, which is set to $10^{-5},\mu$m for the LMS (and should not be changed by the user). The following cell creates the wavelength vector and expands the image WCS that we created above a third, spectral dimension that describes this vector.

wave_min, wave_max = 3.472, 3.526   # microns
dwave = 1e-5                        
wavelength = np.arange(wave_min, wave_max, dwave)
nwave = len(wavelength)

wcs_cube = wcs.sub([1, 2, 0])   # extend image wcs to cube wcs
wcs_cube.wcs.ctype[2] = "WAVE"
wcs_cube.wcs.crpix[2] = 1
wcs_cube.wcs.crval[2] = wave_min
wcs_cube.wcs.cdelt[2] = dwave
wcs_cube.wcs.cunit[2] = "um"

For the flux vector we choose a series of narrow emission lines that are equally spaced in wavelength, on a slowly varying continuum. Note that emission lines need to be Nyquist sampled in the input spectrum! Failing this, scopesim may miss them. As we want to Doppler shift the flux vector we define it on a slightly longer wavelength vector and compute an interpolation function.

lwavelength = np.arange(0.9*wave_min, 1.1*wave_max, dwave)
flux = np.zeros_like(lwavelength)
sigma = 1.5 * dwave
for lamc in [3.48, 3.49, 3.5, 3.51, 3.52]:
    flux += 0.2 * np.exp(-(lwavelength - lamc)**2/(2 * sigma**2))
flux += 0.4 - ((lwavelength - lamc)/0.3)**2

from scipy.interpolate import interp1d
iflux = interp1d(lwavelength, flux)

The following cell is rather compact and does several things:

• np.outer() creates a cube from the wavelength vector and the dloglambda image. The values in this cube are $\lambda/(1+v/c)$, i.e. a shifted wavelength vector at each image position.

• These “wavelengths” are given as an argument to the flux interpolation function, thus filling the cube with values $f(\lambda/(1+v/c))$. This has the effect of shifting the emission lines to observed wavelength $\lambda_0 (1+v/c)$.

• Finally, the cube is multiplied by the image of HL Tau to scale the spectra. This then is our final data cube.

cube = img_cut[None, :, :] * iflux( np.outer(wavelength, (1./(1+dloglambda))).reshape(nwave, ny, nx) )

_, axes = plt.subplots(1, 2, figsize=(9, 4))
axes[0].imshow(cube[802,] - cube[800,], origin='lower')
axes[0].text(10, 170, "difference layers 802 and 800", color='white')
axes[1].plot(cube[:, 90, 91], label="Pixel (91, 90)")
axes[1].plot(cube[:, 90, 89], label="Pixel (89, 90)")
plt.legend()
axes[1].set_xlim(700, 920);

We finally adjust the data values somewhat arbitrarily, by looking at the minimum and maximum values in the cube:

print(cube.min(), cube.max())

Real circumstellar disks have total fluxes of around 100 mJy (optimistically), so we need to scale our cube by about a factor of 50. We’ll make it 100 so SNR won’t be an issue later on.

cube *= 100

total_spec = cube.sum(axis=(1, 2))

plt.plot(wavelength, total_spec)
plt.xlabel("Wavelength [um]")
plt.ylabel("Jy")
plt.title("Integrated spectrum");

We finally put the cube into a fits file. The essential metadata are the WCS, which describes the spatial as well as the spectral dimensions, and the BUNIT keyword. For the latter we stick with Jy, although this is somewhat incongruous as the values are spatially integrated over one pixel, whereas spectrally they are taken per Hz (i.e. a spectral density). Many units are possible as long as they make it clear whether (photon or energy) flux is taken per spatial pixel or solid angle and per spectral bin or wavelength or frequency interval.

hdr = wcs_cube.to_header()
hdr['BUNIT'] = "Jy"

hdu_cube = fits.PrimaryHDU(header=hdr, data=cube.astype(np.float32))

hdu_cube.writeto("our_first_cube.fits", overwrite=True)

In the next notebook, METIS_IFU-02-Simulation.ipynb, we will use this cube as the Source object for METIS-LMS observations.