{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imaging observations of HL Tau"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this notebook is to show how to use a FITS image as input to create a `Source` object for imaging observations with METIS. For L-band we will work with an ALMA 233 GHz image of HL Tau, pretending that this is the structure of the object in the mid-infrared.  For the N-band, we will use a simulation of an AGN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from astropy.io import fits\n",
    "from astropy import units as u\n",
    "from astropy.wcs import WCS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "%matplotlib inline"
   ]
  },
  {
   "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.rc.__config__[\"!SIM.file.local_packages_path\"] = \"../../../\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate observations with METIS, the instrument packages for METIS, the ELT, and Armazones are required. The packages are downloaded from the server and installed into sub-directory `inst_pkgs` in the current working directory. If you have not got the packages yet, uncomment the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sim.download_package(['instruments/METIS', 'telescopes/ELT', 'locations/Armazones'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses two example files that can be downloaded from our server. If you already have these files, make sure that they are located in the current working directory. You can then skip the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.download_example_data([\"HL_Tau_prep_for_Scopesim.fits\", \"AGN_uc0890_image_l12_i090_p000.fits\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a source object from an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image that we use for the L-band is `HL_Tau_prep_for_Scopesim.fits`. This has been derived from `HLTau_B6cont_mscale_ap.image.fits`, which was retrieved from the ESO Science archive. The main step in preparing this file was a rescaling of the image values such that they directly represent flux in Jy per pixel (the flux scale was arbitrarily set such that the average over the first ring corresponds to 0.01 Jy/arcsec2). The image has a pixel scale of 5 mas, which was retained.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_hdul = fits.open(\"HL_Tau_prep_for_Scopesim.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The header contains the information necessary for Scopesim, a WCS and the `BUNIT` keyword that gives the units of the pixel values. Scopesim mainly needs the `CDELT` keywords, which contain the pixel scale of the input image. The reference coordinates `CRVAL` have been set to zero as other values may confuse Scopesim at this stage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_hdul[1].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several ways to instantiate a `Source` object in ScopeSim. In the present case, we provide an image (given as a FITS `HDU`). The header keyword `BUNIT` gives the units of the pixel values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = sim.Source(image_hdu=input_hdul[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## L-band observations of the source"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The source as defined above could be observed in any waveband -- we start in the L-band. The configuration is defined by instantiating a `UserCommands` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd_l = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"img_lm\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default filter is L':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Lp'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cmd_l[\"!OBS.filter_name\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A different filter could be selected by setting this parameter to a different value. Parameters can also be set by giving a properties dictionary, e.g:\n",
    "```python\n",
    "cmd_l = sim.UserCommands(use_instrument=\"METIS\", set_modes=[\"img_lm\"],\n",
    "                         properties={\"!OBS.filter_name\": \"Mp\", \"!OBS.exptime\": 100.})\n",
    "```\n",
    "We shall stick to the L' filter here, and set the exposure time directly in the `readout()` commands later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now build the `OpticalTrain` (which includes atmosphere, telescope and instrument) and `observe` the source. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis_l = sim.OpticalTrain(cmd_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis_l.observe(src, update=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `observe` command creates the ideal image that arrives at the detector focal plane and stores it as\n",
    "```python\n",
    "metis_l.image_planes[0]\n",
    "```\n",
    "This is an HDU object, the image is accessed by appending `.data`.\n",
    "\n",
    "   Photon noise and detector noise are created by the `readout` command, which can be applied several times. In this example we read out with two different exposure times. The exposure time is split automatically into DITs such that the detector is not saturated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_10 = metis_l.readout(exptime=10.)[0]\n",
    "print(\"----\")\n",
    "result_100 = metis_l.readout(exptime=100.)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`readout` creates a list of HDU lists (a list of multi-extension fits files). In this example the list contains a single HDU list, which we have already extracted by adding `[0]` above. The detector image is in the first extension of the HDU list. In the following cells we display (sections of) the detector images, with the input image for comparison, and then save the results to fits files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.subplot(131)\n",
    "plt.imshow(input_hdul[1].data[520:1160, 483:1123], origin='lower', vmin=0, vmax=5e-7,\n",
    "          extent=(-1.6, 1.6, -1.6, 1.6))\n",
    "plt.xlabel(\"[arcsec]\")\n",
    "plt.ylabel(\"[arcsec]\")\n",
    "plt.title(\"Input image\")\n",
    "plt.subplot(132)\n",
    "plt.imshow(result_10[1].data[768:1352, 732:1316], origin='lower', extent=(-1.6, 1.6, -1.6, 1.6))\n",
    "plt.xlabel(\"[arcsec]\")\n",
    "plt.ylabel(\"[arcsec]\")\n",
    "plt.title(r\"Scopesim simulation, $T_\\mathrm{exp} = 10\\,\\mathrm{s}$\")\n",
    "plt.subplot(133)\n",
    "plt.imshow(result_100[1].data[768:1352, 732:1316], origin='lower', extent=(-1.6, 1.6, -1.6, 1.6))\n",
    "plt.xlabel(\"[arcsec]\")\n",
    "plt.ylabel(\"[arcsec]\")\n",
    "plt.title(r\"Scopesim simulation, $T_\\mathrm{exp} = 100\\,\\mathrm{s}$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_10.writeto(\"hl_tau_img_l_10s.fits\", overwrite=True)\n",
    "result_100.writeto(\"hl_tau_img_l_100s.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## N-band observation of the source"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we load the image, here an image from a simulation of an AGN with 1% Eddington ratio, observed edge-on at 12 micron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_hdul = fits.open(\"AGN_uc0890_image_l12_i090_p000.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the image's header to see if it is complete for ScopeSim:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_hdul[0].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is complete, but the pixel scale is very small, since the image was created for VLTI observations. We modify the pixel scale and also adjust the total flux to represent a real object, in this case a powerful AGN with a total extent of about 1\" and a total flux of about 1 Jy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_hdul[0].header[\"CDELT1\"]=1.39e-6\n",
    "input_hdul[0].header[\"CDELT2\"]=1.39e-6 ## increase pixel scaling by factor of 10\n",
    "input_hdul[0].data*=50 ## sets total flux to ca. 1 Jy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us confirm that the total flux is now about 1 Jy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(input_hdul[0].data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now set up a new source object. Unfortunately, the value 'JY/PIXEL' of the header keyword `BUNIT` is not compliant with the FITS standard and is not understood by Scopesim. We could simply change it to the compliant form 'Jy'. Alternatively, the `flux` parameter can be used to stipulate that an image value of 1 corresponds to 1 Jy. Implicitly, this assigns a spectrum that has a constant value of $f_{\\nu}$ to each pixel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = sim.Source(image_hdu=input_hdul[0], flux=1*u.Jy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set up the instrument for N-band imaging, the total requested exposure time is set to 3600 s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd_n = sim.UserCommands(use_instrument='METIS', set_modes=['img_n'], \n",
    "                        properties={\"!OBS.exptime\": 3600})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default filter is N2. As before it could be changed by modifying the configuration keyword"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd_n['!OBS.filter_name']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The instrument itself (the `OpticalTrain`), including atmosphere and telescope, is built by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis_n = sim.OpticalTrain(cmd_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The list of \"effects\" in the `OpticalTrain` can be accessed by `metis_n.effects` (the method `show_in_notebook` adds some eye candy for display in this jupyter notebook). The effects that are actually applied are those with `included=True`, the others are ignored. Note that there are effects that serve the same purpose albeit in different ways. For example, both `armazones_atmo_skycalc_ter_curve` and `armazones_atmo_default_ter_curve` provide atmospheric absorption and emission. The first effect downloads the data from ESO's skycalc server; the second effect uses default data included in the instrument package. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis_n.effects.show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now \"observe\" the source. As before, this command creates the ideal image just in front of the detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis_n.observe(src, update=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The N band configuration of METIS includes a detector effect `ChopNodCombiner`. By default, this effect is switched off so that the readout corresponds to a staring observation, which should be sufficient for most science applications. \n",
    "When the effect is switched on it creates four detector images with chopping and nodding offsets and combines them into a chop-nod image with positive and negative beam images. This corresponds to the primary data product of the METIS data reduction pipeline. \n",
    "The default pattern is perpendicular chopping and nodding with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metis_n['chop_nod'].include = True\n",
    "print(\"Chopping:\", metis_n.cmds['!OBS.chop_offsets'])\n",
    "print(\"Nodding: \", metis_n.cmds['!OBS.nod_offsets'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exposure time given to the `readout()` method is the integration time on one position. With `exptime=3600`, we therefore get a total integration time of 14400 s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outhdul_perpendicular = metis_n.readout(exptime=3600)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndit = metis_n.cmds['!OBS.ndit']\n",
    "dit = metis_n.cmds['!OBS.dit']\n",
    "outhdul_perpendicular[0].header['DIT'] = dit\n",
    "outhdul_perpendicular[0].header['NDIT'] = ndit * 4\n",
    "outhdul_perpendicular[0].header['EXPTIME'] = 4 * dit * ndit\n",
    "outhdul_perpendicular.writeto('agn_er01_90deg_n_img_perpendicular.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "plt.imshow(outhdul_perpendicular[1].data[750:1750, 750:1750], origin='lower', extent=(750, 1750, 750, 1750))\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell combines the four beams into a single image. The chop/nod offsets were 3 arcsec, the N-band pixel scale is 6.79 mas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = outhdul_perpendicular[1].data\n",
    "shift = int(3 / 0.00679)\n",
    "xmin, xmax, ymin, ymax = 823, 1232, 853, 1253\n",
    "beam_1 = image[ymin:ymax, xmin:xmax]\n",
    "beam_2 = image[ymin:ymax, xmin+shift:xmax+shift]\n",
    "beam_3 = image[ymin+shift:ymax+shift, xmin:xmax]\n",
    "beam_4 = image[ymin+shift:ymax+shift, xmin+shift:xmax+shift]\n",
    "final_img = beam_1 - beam_2 - beam_3 + beam_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "hw_arcsec = (xmax - xmin)/2 * 0.00679\n",
    "plt.imshow(final_img, origin='lower', extent=(-hw_arcsec, hw_arcsec, -hw_arcsec, hw_arcsec))\n",
    "plt.xlabel(\"[arcsec]\")\n",
    "plt.ylabel(\"[arcsec]\")\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_wcs = WCS(naxis=2)\n",
    "final_wcs.wcs.ctype = ['LINEAR', 'LINEAR']\n",
    "final_wcs.wcs.crpix = [(final_img.shape[1] - 1)/2, (final_img.shape[0] - 1)/2]\n",
    "final_wcs.wcs.crval = [0., 0.]\n",
    "final_wcs.wcs.cdelt = [0.00679, 0.00679] ## this is the METIS detector plate scale in arcseconds\n",
    "final_wcs.wcs.cunit = ['arcsec', 'arcsec']\n",
    "final_hdr = final_wcs.to_header()\n",
    "final_hdr['EXPTIME'] = 4 * dit * ndit\n",
    "fits.writeto(\"agn_er01_90deg_n_img_reconstructed.fits\", data=final_img, header=final_hdr, overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}