{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook presents METIS LSS L-band simulation of three models of a young stellar object (YSO)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scopesim as sim\n", "sim.bug_report()\n", "\n", "# Edit this path if you have a custom install directory, otherwise comment it out. [For ReadTheDocs only]\n", "sim.link_irdb(\"../../../\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "from astropy.io import fits\n", "from astropy import units as u\n", "\n", "from matplotlib import pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you haven't got the instrument packages yet, uncomment the following cell, which will install the packages into `./inst_pkgs`, a subdirectory of your current working directory. This is the default location where scopesim looks for instrument packages." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# sim.download_packages([\"METIS\", \"ELT\", \"Armazones\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have downloaded the packages but to a different location, you can set\n", "```python\n", "sim.set_inst_pkgs_path(\"/path/to/inst/pkgs\")\n", "```\n", "We recommend, however, to create a working directory for each simulation project and to use the default installation of packages into a subdirectory. Keeping simulation results and configuration files together that way makes it easy to reconstruct later the exact conditions under which a simulation was run." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Preparation of source cubes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The input data are cubes of three different models of the same YSO, HD100546. We keep the names of FITS files, the `Source` objects and the results of the ScopeSim simulations in dictionaries, indexed by short names for the models." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fitsfiles = {\n", " \"cav\": \"models_Lband_HD100546_cav_f100PAH.cube_3.0mas.fits\",\n", " \"emptycav\": \"models_Lband_HD100546_empytcav.cube_3.0mas.fits\",\n", " \"gap\": \"models_Lband_HD100546_gap100.cube_3.0mas.fits\",\n", "}\n", "print(\"Model names:\", list(fitsfiles))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The FITS files can be downloaded from the ScopeSim server. If you already have them, make sure that the files are in the correct location (e.g. current working directory, see also the note below). The next code cell will replace the file names with absolute paths to the download cache location. If you already have the files in the current working directory, simply skip that line and ScopeSim will look for them there." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Note: ScopeSim v0.9.0 or later will now download the example data to a hidden cache directory by default. To change this, pass the optional save_dir argument to sim.download_example_data() with the desired download location, e.g. save_dir=\"./\" for the current working directory (the default download location in previous versions).\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fitsfiles = dict(zip(fitsfiles, sim.download_example_data(*fitsfiles.values())))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The file headers are not yet in the form that ScopeSim understands and we make two minor modifications: \n", "- Set CRVAL to 0, because Scopesim cannot look elsewhere\n", "- Set BUNIT keyword (the files have the keyword UNITS, which is non-standard)\n", "- The cubes contain some negative values. We replace these with 0.\n", "- We introduce a factor `scale_delt` to increase the pixel size, which makes features more visible. If you want to simulate the original source pixel scale, set `scale_delt` to 1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sources = {}\n", "scale_cdelt = 5.\n", "for model, fitsfile in fitsfiles.items():\n", " with fits.open(fitsfile) as hdul:\n", " hdul[0].header['CRVAL1'] = 0.\n", " hdul[0].header['CRVAL2'] = 0.\n", " hdul[0].header['CDELT1'] *= scale_cdelt\n", " hdul[0].header['CDELT2'] *= scale_cdelt\n", " hdul[0].header['BUNIT'] = hdul[0].header['UNITS']\n", " hdul[0].data.clip(min=0, out=hdul[0].data)\n", " sources[model] = sim.Source(cube=hdul)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get an impression of what the data look like we display the cubes summed along the wavelength direction." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Determine plot limits in arcsec from header keywords (these are in degrees)\n", "hdr = sources['cav'].cube_fields[0].header\n", "i_lim = np.array([0, hdr['NAXIS1']])\n", "x_lim = hdr['CRVAL1'] + hdr['CDELT1'] * (i_lim + 1 - hdr['CRPIX1']) * 3600\n", "j_lim = np.array([0, hdr['NAXIS2']])\n", "y_lim = hdr['CRVAL2'] + hdr['CDELT2'] * (j_lim + 1 - hdr['CRPIX2']) * 3600\n", "\n", "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 4))\n", "for i, (model, src) in enumerate(sources.items()):\n", " im = axes[i].imshow(src.cube_fields[0].data.sum(axis=0) + 1e-14, # add small positive value to avoid 0 in LogNorm\n", " origin=\"lower\", norm=\"log\", vmin=1e-4, vmax=10,\n", " extent=(x_lim[0], x_lim[-1], y_lim[0], y_lim[-1]))\n", " axes[i].set_xlabel(\"arcsec\")\n", " axes[i].set_ylabel(\"arcsec\")\n", " axes[i].set_title(model);\n", "fig.subplots_adjust(right=0.9)\n", "cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75])\n", "fig.colorbar(im, cax=cbar_ax);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In wavelength, the cubes are sampled on a linear wavelength grid from 3.1 to 4 $\\mu\\mathrm{m}$, with a step size of $0.2\\,\\mu\\mathrm{m}$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(sources['cav'].cube_fields[0].waveset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the source object does not cover the entire spatial and wavelength range permitted by METIS: The spatial extent is about 2.5 arcsec, compared to the METIS slit length of 8 arcsec. The wavelength range is 3.1 to 4 $\\mu$m, whereas METIS permits 2.9 to 4.2 $\\mu$m. This will be visible in the simulated raw and rectified spectra below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation with Scopesim" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cubes are observed in the long-slit spectroscopic mode in the L band. As usual, there are four steps: `UserCommands` -> `OpticalTrain` -> `observe` -> `readout` to arrive at a detector image. The optical train can be reused for observation of different sources when `update=True` is set in `observe()`. \n", "\n", "In addition to observing the three model cubes we will make a blank sky observation for background subtraction, and we want to ensure that all observations are made with the same detector settings so that the background can be subtracted without any rescaling. The detector settings are (1) the detector mode, which can be \"fast\" or \"slow\", with different gain values and noise characteristics; and (2) the breakdown of the requested exposure time into `NDIT` subexposures of integration time `DIT`. The parameters can be chosen automatically by the effects `metis['detector_readout_parameters']` and `metis['auto_exposure']`, optimising signal-to-noise ratio while avoiding detector saturation. We shall determine these parameters through an observation of the first model and then apply them to all observations, including the sky observation for which the effects would choose different parameters otherwise." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "exptime = 3600 * u.s\n", "cmd = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"lss_l\"],\n", " properties={\"!OBS.exptime\": exptime.value, \"!OBS.dit\": None, \"!OBS.ndit\": None})" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "metis = sim.OpticalTrain(cmd)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get exposure parameters\n", "metis.observe(sources['cav'])\n", "tmpresult = metis.readout(detector_readout_mode='auto', dit=None, ndit=None)[0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "detmode = tmpresult[0].header[\"HIERARCH ESO DET1 MODE\"]\n", "dit = tmpresult[0].header[\"HIERARCH ESO DET1 DIT\"]\n", "ndit = tmpresult[0].header[\"HIERARCH ESO DET1 NDIT\"]\n", "\n", "print(\"Detector mode:\", detmode)\n", "print(\"DIT:\", dit)\n", "print(\"NDIT:\", ndit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now set these parameters in the calls to `metis.readout()` after observing all models." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results = {}\n", "for model, src in sources.items():\n", " print(f'Observing model \"{model}\"...')\n", " metis.observe(src, update=True)\n", " results[model] = metis.readout(detector_readout_mode=detmode, dit=dit, ndit=ndit)[0]\n", " print(\"-----\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5.2))\n", "for axis, (model, result) in zip(axes, results.items()):\n", " axis.imshow(result[1].data, origin=\"lower\", norm=\"log\")\n", " axis.set_xlabel(\"pixel\")\n", " axis.set_ylabel(\"pixel\")\n", " axis.set_title(model)\n", "fig.suptitle(\"Raw spectra\", fontsize=20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following is a simulation of a blank-sky observation, which will be used for background subtraction. We use the same readout parameters as for the models. Given that the source is fairly compact, one could estimate the background from the science observation as well." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sky = sim.source.source_templates.empty_sky()\n", "metis.observe(sky, update=True)\n", "bgresult = metis.readout(detector_readout_mode=detmode, dit=dit, ndit=ndit)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We subtract the background from the data before rectifying the spectra." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for result in results.values():\n", " result[1].data -= bgresult[1].data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now use the `rectify_traces` method to resample the spectra to a linear grid in wavelength and spatial position. The method is attached to the `SpectraTraceList` effect, which holds the information on the geometrical mapping onto the detector. Note that this is an optimistic approach as it uses the same transformations used before to simulate the detector data. When reducing real data, the transformation parameters would have to be estimated from dedicated calibration measurements and would therefore have a degree of uncertainty." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rectified = {}\n", "tracelist = metis['spectral_traces']\n", "for model, result in results.items():\n", " rectified[model] = tracelist.rectify_traces(result, -4.0, 4.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before showing the rectified spectra, we convert the pixel numbers to wavelength and spatial position along the slit. For this purpose we use the WCS keywords in the headers of the rectified HDUs. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.wcs import WCS\n", "hdr = rectified['cav'][1].header\n", "wcs = WCS(hdr)\n", "naxis1, naxis2 = hdr['NAXIS1'], hdr['NAXIS2']\n", "det_xi = wcs.sub((2,)).pixel_to_world(np.arange(naxis2)).to(u.arcsec)\n", "det_wave = wcs.spectral.pixel_to_world(np.arange(naxis1)).to(u.um) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(15, 15))\n", "for axis, (model, result) in zip(axes, rectified.items()):\n", " img = axis.imshow(\n", " result[1].data, origin=\"lower\", norm=\"log\", aspect=\"auto\",\n", " extent=(det_wave[0].value, det_wave[-1].value, det_xi[0].value, det_xi[-1].value))\n", " axis.set_ylabel(f\"position along slit [{det_xi.unit.to_string('latex')}]\")\n", " axis.set_xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\")\n", " axis.set_title(model);\n", " fig.colorbar(img, ax=axis)\n", "fig.suptitle(\"Rectified spectra\", y=.99, fontsize=24)\n", "fig.tight_layout();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As stated above the input source object does not fill the entire spatial and wavelength ranges covered by METIS. In a real observation the entire frame shown here would be filled." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results have to be saved to disk explicitely so they can be analysed with external tools. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "# Change to True if you want to save the rectified data\n", "save_rectified = False\n", "if save_rectified:\n", " for model, result in rectified.items():\n", " outfile = Path(fitsfiles[model]).stem + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n", " result.writeto(outfile, overwrite=True)\n", " print(fitsfiles[model], \"--->\", outfile)\n", " bgresult.writeto(\"models_Lband_HD100546-background_simulated_LSS_L.fits\", overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Rotating the field\n", "In Scopesim, the METIS slit is always aligned with the x-axis. To simulate different slit orientations, the input cube has to be rotated. We define a function to do this (for more explanation of the steps see the notebook `LSS_AGN-01_preparation.ipynb`)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.wcs import WCS\n", "from scipy.interpolate import RectBivariateSpline" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def rotate_cube(incube, angle):\n", " \"\"\"Rotate input cube by angle\n", " \n", " Parameters\n", " ----------\n", " incube : ImageHDU\n", " angle : float\n", " angle in degrees, counterclockwise from positive x-axis\n", " \"\"\"\n", " rangle = angle * np.pi / 180 # degrees to radians\n", " \n", " wcs_orig = WCS(incube.header).sub(2)\n", " wcs_orig.wcs.ctype = [\"LINEAR\", \"LINEAR\"] # avoids discontinuity around RA=0 degrees\n", " wcs_rot = WCS(incube.header).sub(2)\n", " wcs_rot.wcs.ctype = [\"LINEAR\", \"LINEAR\"]\n", " \n", " wcs_rot.wcs.pc = [[np.cos(rangle), np.sin(rangle)],\n", " [-np.sin(rangle), np.cos(rangle)]]\n", " i_orig = np.arange(incube.header['NAXIS1'])\n", " j_orig = np.arange(incube.header['NAXIS2'])\n", " x_orig, y_orig = wcs_orig.all_pix2world(i_orig, j_orig, 0)\n", " x_orig[x_orig >180] -= 360. # RA continuous across 0\n", " i_rot, j_rot = np.meshgrid(i_orig, j_orig)\n", " x_rot, y_rot = wcs_rot.all_pix2world(i_rot, j_rot, 0)\n", " \n", " for iplane, plane in enumerate(incube.data):\n", " interp = RectBivariateSpline(y_orig, x_orig, plane, kx=1, ky=1)\n", " incube.data[iplane,] = interp(y_rot, x_rot, grid=False)\n", " incube.header['ANGLE'] = angle, \"slit rotation [deg]\"\n", " return incube" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function is applied to the data cubes before creating the `Source` object. We choose an angle of 45 degrees. The simulation then proceeds as above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "angle = 45\n", "rotsources = {}\n", "for model, fitsfile in fitsfiles.items():\n", " with fits.open(fitsfile) as hdul: \n", " hdul[0].header['CRVAL1'] = 0.\n", " hdul[0].header['CRVAL2'] = 0.\n", " hdul[0].header['CDELT1'] *= scale_cdelt\n", " hdul[0].header['CDELT2'] *= scale_cdelt\n", " hdul[0].header['BUNIT'] = hdul[0].header['UNITS']\n", " hdul[0].data.clip(min=0, out=hdul[0].data)\n", " hdul[0] = rotate_cube(hdul[0], angle)\n", " rotsources[model] = sim.Source(cube=hdul)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rotresults = {}\n", "for model, src in rotsources.items():\n", " print(f'Observing model \"{model}\"...')\n", " metis.observe(src, update=True)\n", " rotresults[model] = metis.readout(detector_readout_mode=detmode, dit=dit, ndit=ndit)[0]\n", " print(\"-----\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rotrectified = {}\n", "for model, result in rotresults.items():\n", " result[1].data -= bgresult[1].data\n", " rotrectified[model] = tracelist.rectify_traces(result, -4.0, 4.0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(15, 15))\n", "for axis, (model, result) in zip(axes, rotrectified.items()):\n", " img = axis.imshow(\n", " result[1].data, origin=\"lower\", norm=\"log\", aspect=\"auto\",\n", " extent=(det_wave[0].value, det_wave[-1].value, det_xi[0].value, det_xi[-1].value))\n", " axis.set_ylabel(f\"position along slit [{det_xi.unit.to_string('latex')}]\")\n", " axis.set_xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\")\n", " axis.set_title(model);\n", " fig.colorbar(img, ax=axis)\n", "fig.suptitle(\"Rectified rotated spectra\", y=.99, fontsize=24)\n", "fig.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Change to True if you want to save the rectified data\n", "save_rotrectified = False\n", "if save_rotrectified:\n", " for model, result in rotresults.items():\n", " outfile = Path(fitsfiles[model]).stem + f\"-rot_{angle}-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n", " result.writeto(outfile, overwrite=True)\n", " print(fitsfiles[model], \"--->\", outfile)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compare the rotated spectrum for the `emptycav` model to the spectrum without rotation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6, 12))\n", "\n", "ax1.plot(det_xi, rectified['emptycav'][1].data[:, 1500], label=\"no rotation\")\n", "ax1.plot(det_xi, rotrectified['emptycav'][1].data[:, 1500], label=f\"rotation {angle} deg\")\n", "ax1.set_xlim(x_lim[0], x_lim[-1])\n", "ax1.set_ylim(1e-1, 3e4)\n", "ax1.set_xlabel(\"arcsec\")\n", "ax1.semilogy()\n", "ax1.legend()\n", "\n", "fig.subplots_adjust(left=0.1)\n", "ax2.imshow(rotsources['emptycav'].cube_fields[0].data.sum(axis=0) + 1e-14,\n", " norm=\"log\", vmin=1e-4, vmax=10,\n", " extent=(x_lim[0], x_lim[-1], y_lim[0], y_lim[-1]))\n", "ax2.set_xlabel(\"arcsec\")\n", "ax2.set_ylabel(\"arcsec\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Flux calibration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The pixel values in the detector images give electrons accumulated over the exposure time. To get back to physical units, e.g. Jy, one has to perform a flux calibration. As in real observations, we do this here with the observation of a standard star. For simplicity, we use a star with a spectrum that is constant at $f_{\\nu} = 1\\,\\mathrm{Jy}$. We observe it with the identical setup as the science targets above, except that the exposure time is reduced to 1 second." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "std = sim.source.source_templates.star(flux=1 * u.Jy)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "metis.observe(std, update=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "std_exptime = 1 * u.s\n", "std_result = metis.readout(exptime=std_exptime.value, detector_readout_mode=detmode, dit=None, ndit=None)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We subtract the background (scaled to the exposure time of the standard exposure) and rectify." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "std_dit = std_result[0].header[\"HIERARCH ESO DET1 DIT\"]\n", "std_result[1].data -= bgresult[1].data / dit * std_dit\n", "std_rectified = tracelist.rectify_traces(std_result, -4, 4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(15, 7))\n", "plt.imshow(std_rectified[1].data, origin=\"lower\", norm=\"log\", aspect=\"auto\",\n", " extent=(det_wave[0].value, det_wave[-1].value, det_xi[0].value, det_xi[-1].value))\n", "plt.ylabel(f\"position along slit [{det_xi.unit.to_string('latex')}]\")\n", "plt.xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\")\n", "plt.title(\"Standard star\")\n", "plt.colorbar();" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xmin, xmax = 675, 825\n", "xaxis = np.arange(xmin, xmax)\n", "y_1, y_2 = 200, 1750\n", "lam_1, lam_2 = det_wave[y_1], det_wave[y_2]\n", "\n", "plt.figure(figsize=(10, 5))\n", "plt.plot(xaxis, std_rectified[1].data[xmin:xmax, y_1], label=fr\"$\\lambda = ${lam_1.to_string(precision=3, format='latex')}\")\n", "plt.plot(xaxis, std_rectified[1].data[xmin:xmax, y_2], label=fr\"$\\lambda = ${lam_2.to_string(precision=3, format='latex')}\")\n", "plt.plot(xaxis, std_rectified[1].data[xmin:xmax, :].mean(axis=1), label='average')\n", "plt.vlines((732 - 10, 732 + 10), 1, 1e5, colors='black', linestyles='dashed', label='extraction aperture')\n", "plt.xlabel(\"pixel\")\n", "plt.ylabel(\"ADU\")\n", "plt.semilogy()\n", "plt.ylim(10, 4e4)\n", "plt.legend()\n", "plt.title(\"Spatial cuts through standard spectrum\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The total flux density of 1 Jy from the star is spread over the point spread function, which extends far from the central position. It is cut by the finite aperture of the spectroscopic slit, and will be cut further by the finite extent of the extraction aperture over which we will sum the two-dimensional spectrum. Most of the flux is contained in the core of the PSF; from the figure, we can take an aperture size in the spatial direction of 20 pixels - outside this aperture, the spectrum is dominated by noise. \n", "The slit that we used is" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(metis['slit_wheel'].current_slit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with a width of 38.1 mas, corresponding to 7 pixels. As this notebook focuses on the simulator rather than the data reduction, we neglect slit and aperture losses and assume that the signal integrated over 20 pixels in the spatial direction corresponds to the input flux density of 1 Jy. The error is on the order of 2 per cent. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hwidth = 10 # pixels, half width of extraction aperture\n", "# exptime = std_rectified[0].header[\"EXPTIME\"] * u.s\n", "std_1d = std_rectified[1].data[731-hwidth:731+hwidth].sum(axis=0) / std_dit" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(det_wave, std_1d)\n", "plt.xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\")\n", "plt.ylabel(\"ADU\")\n", "plt.title(fr\"Extracted standard spectrum ($t_\\mathrm{{exp}} = ${std_exptime.to_string(format='latex')})\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We construct the flux conversion curve that takes the observed spectrum to the calibrated spectrum in Jy. Note that there are strong excursions at the edge of the spectroscopic filter and at deep atmospheric absorption features. The object flux at these wavelengths is not recoverable." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "flux_conv = 1 * u.Jy / std_1d" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(det_wave, flux_conv)\n", "plt.ylim(0, 1e-5)\n", "plt.xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\")\n", "plt.ylabel(flux_conv.unit.to_string(\"latex_inline\"))\n", "plt.title(\"Flux conversion\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We apply the flux calibration to the two-dimensional spectra from above, taking into account the different exposure time, and display them. Note how the PAH feature at around 3.3 um now becomes apparent in the disk." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "emptycav_Jy = flux_conv * rectified['emptycav'][1].data * u.electron / exptime\n", "emptycav_Jy_rot = flux_conv * rotrectified['emptycav'][1].data * u.electron / exptime" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(13, 12))\n", "ax1.imshow(emptycav_Jy.value, origin=\"lower\", norm=\"log\", vmin=1e-5, vmax=6e-5, aspect=\"auto\",\n", " extent=(det_wave[0].value, det_wave[-1].value, det_xi[0].value, det_xi[-1].value))\n", "ax1.set_ylabel(f\"position along slit [{det_xi.unit.to_string('latex')}]\")\n", "ax1.set_xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\")\n", "ax1.set_title(f\"{model}, no rotation\")\n", "\n", "img = ax2.imshow(emptycav_Jy_rot.value, origin=\"lower\", norm=\"log\", vmin=1e-5, vmax=6e-5, aspect=\"auto\", \n", " extent=(det_wave[0].value, det_wave[-1].value, det_xi[0].value, det_xi[-1].value))\n", "ax2.set_ylabel(f\"position along slit [{det_xi.unit.to_string('latex')}]\")\n", "ax2.set_xlabel(f\"Wavelength [{det_wave.unit.to_string('latex')}]\") \n", "ax2.set_title(f\"{model}, rotated by {angle} degrees\")\n", "\n", "fig.subplots_adjust(right=0.8)\n", "cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n", "fig.colorbar(img, cax=cbar_ax, label=emptycav_Jy_rot.unit);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we perform background subtraction and flux calibration on all simulations. The data are now normalised to an exposure time of 1 second and therefore have units 'electrons per second'." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Change to True if you want to save the calibrated data\n", "save_calibrated = False\n", "for model, result in rectified.items():\n", " result[1].data = flux_conv.value * result[1].data / exptime\n", " result[1].header['BUNIT'] = \"e/s\"\n", " if save_calibrated:\n", " outfile = (Path(fitsfiles[model]).stem + \"-simulated_LSS_L-bgsub_fluxcal\"\n", " + Path(fitsfiles[model]).suffix)\n", " result.writeto(outfile, overwrite=True)\n", " print(\"--->\", outfile)\n", "\n", "for model, result in rotrectified.items():\n", " result[1].data = flux_conv.value * result[1].data / exptime\n", " result[1].header['BUNIT'] = \"e/s\"\n", " if save_calibrated:\n", " outfile = (Path(fitsfiles[model]).stem + f\"-rot_{angle}-simulated_LSS_L-bgsub_fluxcal\"\n", " + Path(fitsfiles[model]).suffix)\n", " result.writeto(outfile, overwrite=True)\n", " print(\"--->\", outfile)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Comparison to input spectrum" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "result_extract = rectified['emptycav'][1].data[731-hwidth:731+hwidth].sum(axis=0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "src_extract = sources['emptycav'].cube_fields[0].data[:, 88:91, 87:92].sum(axis=(1, 2))\n", "src_wave = sources['emptycav'].cube_fields[0].waveset" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "skycalc = metis['skycalc_atmosphere']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(12, 6))\n", "plt.plot(det_wave, result_extract, label=\"flux-calibrated simulated spectrum\")\n", "# The scaling for the source spectrum has been estimated by eye\n", "plt.plot(src_wave, src_extract / 31000, label=\"input spectrum\")\n", "plt.plot(det_wave, 0.0002 * skycalc.throughput(det_wave), label=\"sky transmission (* 0.0002)\")\n", "plt.xlabel(f\"Wavelength [{src_wave.unit.to_string('latex')}]\")\n", "plt.ylabel(\"Relative flux\")\n", "plt.legend()\n", "plt.xlim(3.1, 4.0);" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.9" } }, "nbformat": 4, "nbformat_minor": 4 }