notebooks/miri_photon_noise_and_contrast.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates the use of the Pandeia engine with the MIRI coronagraphs. Specifically, it provides examples of:\n",
    "* Scene construction and instrument configuration\n",
    "* Engine calculations (using the bundled precomputed PSF library and with on-the-fly PSFs with WebbPSF)\n",
    "* Basic data reduction with PSF subtraction\n",
    "* Iterating over instrument configurations to optimize observations\n",
    "* Contrast calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from pandeia.engine.calc_utils import build_default_calc\n",
    "\n",
    "import jwst_pancake as pancake\n",
    "\n",
    "from copy import deepcopy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constructing a Scene\n",
    "\n",
    "We'll start by defining the source and instrument properties for our desired observation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_mV = 4.7\n",
    "ref_mV = 3.\n",
    "\n",
    "target_Sp = 'a5v'\n",
    "ref_Sp = 'a3v'\n",
    "\n",
    "subarray = 'mask1065'\n",
    "filter_c = 'f1065c'\n",
    "mask_c = 'fqpm1065'\n",
    "\n",
    "ngroup = 10\n",
    "nint = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll load in a MIRI template and configure the instrument for our observation. (In this case, it's already configured for a F1065C observation, but we'll demonstrate setting the calculation keys regardless.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a default MIRI coronagraphy calculation\n",
    "config = build_default_calc('jwst', 'miri', 'coronagraphy')\n",
    "\n",
    "# Set the coronagraph and filter\n",
    "config['configuration']['detector']['subarray'] = subarray\n",
    "config['configuration']['detector']['ngroup'] = ngroup\n",
    "config['configuration']['instrument']['aperture'] = mask_c\n",
    "config['configuration']['instrument']['filter'] = filter_c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This template contains a scene with a single star. We'll set the star properties and then duplicate it to create a planetary companion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pull out the target (the first entry in the 'scene' list)\n",
    "targetstar = config['scene'][0]\n",
    "targetstar['id'] = 1\n",
    "targetstar['spectrum'] = {'spectrum_parameters': ['normalization', 'sed']}\n",
    "targetstar['spectrum']['normalization'] = {'type': 'photsys', 'bandpass': 'johnson,j', 'norm_fluxunit': 'abmag'}\n",
    "targetstar['spectrum']['normalization']['norm_flux'] = target_mV\n",
    "targetstar['spectrum']['sed'] = {'sed_type': 'phoenix', 'key': target_Sp}\n",
    "\n",
    "# Copy the target star and turn it into a planet\n",
    "planetA = deepcopy(targetstar)\n",
    "planetA['id'] = 2 #each source must have a unique ID, starting at 1\n",
    "\n",
    "# A different way to normalize source flux\n",
    "planetA['spectrum']['normalization']['bandpass'] = 'miri,imaging,f1500w'\n",
    "planetA['spectrum']['normalization']['norm_flux'] = 15.3\n",
    "planetA['spectrum']['normalization']['type'] = 'jwst'\n",
    "planetA['spectrum']['sed']['sed_type'] = 'blackbody'\n",
    "planetA['spectrum']['sed']['temp'] = 900.\n",
    "del planetA['spectrum']['sed']['key'] #unnecessary now\n",
    "\n",
    "# Source offset\n",
    "planetA['position']['x_offset'] = 2.3 #arcsec\n",
    "planetA['position']['y_offset'] = 3.1\n",
    "\n",
    "# Update calculation file with the new planet\n",
    "config['scene'].append(planetA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rotate the scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pancake.scene.rotate_scene(config['scene'],35.,center=[0.,0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll add in a global offset to the entire scene to capture the effect of target acquisition error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errx, erry = pancake.scene.get_ta_error()\n",
    "pancake.scene.offset_scene(config['scene'], errx, erry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now for the reference scene for PSF subtraction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "refstar = config['strategy']['psf_subtraction_source']\n",
    "refstar['spectrum'] = deepcopy(targetstar['spectrum'])\n",
    "refstar['spectrum']['normalization']['norm_flux'] = ref_mV\n",
    "refstar['spectrum']['sed']['key'] = ref_Sp\n",
    "refstar['id'] = 3\n",
    "\n",
    "# And add target acquisition error\n",
    "errx_ref, erry_ref = pancake.scene.get_ta_error()\n",
    "pancake.scene.offset_scene([refstar], errx_ref, erry_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we'll plot the two scenes we've constructed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,5))\n",
    "plt.subplot(121,polar=True)\n",
    "pancake.scene.plot_scene(config['scene'],'Target Scene w/ Companion',newfig=False)\n",
    "ax = plt.gca()\n",
    "ax.set_rlim(0,5.)\n",
    "plt.subplot(122,polar=True)\n",
    "pancake.scene.plot_scene([refstar],'Reference Scene',newfig=False)\n",
    "ax = plt.gca()\n",
    "ax.set_rlim(0,5.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the Pandeia Engine\n",
    "\n",
    "Now we pass our calculation files to the pandeia engine to create the slope images (and a number of other products).\n",
    "\n",
    "The most direct way to do so is to pass the calculation file to ```engine.perform_calculation.``` However, because we want to calculate both the target and the reference source, using the convenience functions ```engine.calculate_target``` and ```engine.calculate_reference``` will take the appropriate parts of the coronagraphy configuration and calculate them.\n",
    "\n",
    "### Wave Sampling\n",
    "\n",
    "An aside on performance and accuracy: The ```engine.options.wave_sampling``` parameter provides a hook into the wavelength sampling of the 3D (x,y,wavelength) cube. By default, Pandeia adopts some large value for the wavelength sampling (typically 150+); however, this is the primary time sink in the calculation. Setting ```engine.options.wave_sampling = 10``` while developing your simulation provides dramatic time savings while getting within ~5% of the \"true\" value. By ```engine.options.wave_sampling = 40```, one can expect agreement to within roughly 1%.\n",
    "\n",
    "### On-the-fly PSF Calculations\n",
    "\n",
    "The Pandeia engine relies on a library of precomputed PSFs that are sparsely sampled across the coronagraphic field of view. For the MIRI coronagraphs, this sparse sampling will often be insufficient for accurately capturing PSF variations arising from small offsets.\n",
    "\n",
    "Pandeia-Coronagraphy gives the option (```engine.options.on_the_fly_PSFs```) to circumvent the use of this precomputed library and force recomputing each PSF on the fly in WebbPSF. We'll toggle it to True, at the cost of some speed.\n",
    "\n",
    "### Saving and loading Options\n",
    "\n",
    "The options object offers several convenience functions to simplify its use. The ```options.save_options``` function stores a single set of options in the object for later retrieval (this allows you to adjust the options several times to see the result whilst still preserving the original configuration), and ```options.restore_options``` restores whatever is kept in that single internal backup. If you need to save several sets of configurations, or store configurations for later use, the ```options.current_options``` property contains a dictionary describing the current configuration state, and its values can be saved within a script. In addition, by setting ```options.current_options``` to a dictionary of your choice (whether created manually or obtained from the current_options property), you can set the calculation parameters as you choose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options = pancake.engine.options\n",
    "options.save_options()\n",
    "options.wave_sampling = 21\n",
    "options.on_the_fly_PSFs = True\n",
    "\n",
    "results = pancake.engine.calculate_all(config)\n",
    "targ_results = results['target']\n",
    "ref_results = results['reference']\n",
    "options.restore_options()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the detector images (in e$^-$/s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_slope = targ_results['2d']['detector']\n",
    "reference_slope = ref_results['2d']['detector']\n",
    "\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(121)\n",
    "plt.imshow(target_slope)\n",
    "plt.title('Target Slope Image')\n",
    "plt.colorbar().set_label('e$^{-}$/s')\n",
    "plt.subplot(122)\n",
    "plt.imshow(reference_slope)\n",
    "plt.title('Reference Slope Image')\n",
    "plt.colorbar().set_label('e$^{-}$/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And check for saturation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_sat = targ_results['2d']['saturation']\n",
    "reference_sat = ref_results['2d']['saturation']\n",
    "\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(121)\n",
    "plt.imshow(target_sat)\n",
    "plt.title('Target Saturation Image')\n",
    "plt.subplot(122)\n",
    "plt.imshow(reference_sat)\n",
    "plt.title('Reference Saturation Image')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Post-Processing\n",
    "\n",
    "Subtract the registered and scaled reference PSF from the target image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Handle any NaNs in the image (caused by saturation)\n",
    "clean_reference = reference_slope.copy()\n",
    "clean_reference[np.isnan(reference_slope)] = np.nanmax(reference_slope)\n",
    "\n",
    "# Register the reference to the target and renormalize\n",
    "registered_ref = pancake.analysis.register_to_target(clean_reference,target_slope,rescale_reference=True)\n",
    "centered_target = target_slope - np.nanmean(target_slope)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(centered_target - registered_ref)\n",
    "plt.colorbar().set_label('e$^{-}$/s')\n",
    "plt.title('After Reference-Subtraction')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimizing Observation Parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngroup_list = range(2,8)\n",
    "\n",
    "raw_images = []\n",
    "refsub_images = []\n",
    "for ngroup in ngroup_list:\n",
    "    \n",
    "    #Exposure parameters\n",
    "    numberofgroups = ngroup\n",
    "    numberofints = 1\n",
    "    current_config = deepcopy(config)\n",
    "    current_config['configuration']['detector']['ngroup'] = numberofgroups\n",
    "    current_config['configuration']['detector']['nint'] = numberofints\n",
    "\n",
    "    #Pandeia calculation\n",
    "    # Note that accurately capturing spatial variation is unnecessary for this calculation\n",
    "    # As such, it isn't necessary to turn on on-the-fly PSFs\n",
    "    current_result = pancake.engine.calculate_all(current_config)\n",
    "    occ_slope = current_result['target']['2d']['detector']\n",
    "    ref_slope = current_result['reference']['2d']['detector']\n",
    "\n",
    "    #PSF subtraction assuming photon noise, the normalization is done properly \n",
    "    centered_occ = occ_slope - np.nanmean(occ_slope)\n",
    "    centered_ref = ref_slope - np.nanmean(ref_slope)\n",
    "    reg_ref = pancake.analysis.register_to_target(centered_ref, centered_occ)\n",
    "    ref_sub = centered_occ - reg_ref\n",
    "    \n",
    "    raw_images.append(centered_occ)\n",
    "    refsub_images.append(ref_sub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for raw, reduced, groups in zip(raw_images, refsub_images, ngroup_list):\n",
    "    plt.figure(figsize=(10,4))\n",
    "    plt.subplot(121)\n",
    "    plt.imshow(raw)\n",
    "    plt.title('Raw Data ({} groups)'.format(groups))\n",
    "    plt.colorbar().set_label('e$^{-}$/s')\n",
    "    plt.subplot(122)\n",
    "    plt.imshow(reduced)\n",
    "    plt.title('Reduced Data ({} groups)'.format(groups))\n",
    "    plt.colorbar().set_label('e$^{-}$/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save Calculation File\n",
    "\n",
    "Save out your scene and instrument parameters for quick loading with a future call to ```engine.load_calculation```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pancake.engine.save_calculation(config, 'mygreatcalculation.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save Pandeia Images\n",
    "\n",
    "```engine.save_to_fits``` is provided as a convenience function for quickly saving out arrays or cubes to a FITS file. This doesn't preserve any header values. See http://docs.astropy.org/en/stable/io/fits/ for a more complete treatment of reading and writings FITS files in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save out 2d slop images\n",
    "pancake.engine.save_to_fits(targ_results['2d']['detector'], 'target_slope.fits')\n",
    "pancake.engine.save_to_fits(ref_results['2d']['detector'], 'reference_slope.fits')\n",
    "\n",
    "# Save out cube\n",
    "pancake.engine.save_to_fits(raw_images,'raw_cube.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Limiting Contrast Calculation\n",
    "\n",
    "\n",
    "This is where we calculate the contrast as a function of exposure time. Note that in this example the limiting contrast is driven by:\n",
    "\n",
    "1. detector noise, in the Pandeia model\n",
    "2. photon noise on the speckles\n",
    "3. TA error between target and reference. \n",
    "\n",
    "There is no error due to wavefront thermal drifts or dynamical vibrations.\n",
    "\n",
    "For the occulted source, we'll copy the target calculation from before and pop the planet out of the scene list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "occulted = deepcopy(config)\n",
    "occulted['scene'].pop(-1)\n",
    "\n",
    "#and plot\n",
    "pancake.scene.plot_scene(occulted['scene'],'Occulted Star')\n",
    "ax = plt.gca()\n",
    "ax.set_rlim(0,1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unocculted source for the contrast normalization is just the target source moved from behind the coronagraphic mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#copy the occulted calculation file\n",
    "unocculted = deepcopy(occulted)\n",
    "#apply an offset\n",
    "unocculted['scene'][0]['position']['x_offset'] = 0.8 # arcsec\n",
    "unocculted['scene'][0]['position']['y_offset'] = 0.8 # arcsec\n",
    "\n",
    "#and plot\n",
    "pancake.scene.plot_scene(unocculted['scene'],'Unocculted Star')\n",
    "ax = plt.gca()\n",
    "ax.set_rlim(0,1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " We'll use the same reference as in the previous calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import fftconvolve\n",
    "\n",
    "def quick_contrast(occulted_image, unocculted_image, n_annuli=20):\n",
    "    '''\n",
    "    A quick and dirty contrast calculation just for demonstration\n",
    "    purposes.\n",
    "    '''\n",
    "    # Convolve the unocculted image with an aperture and pick out the max \n",
    "    # as the normalization constant\n",
    "    kernel = np.array([[0,0,1,0,0], #simple aperture\n",
    "                   [0,1,1,1,0],\n",
    "                   [1,1,1,1,1],\n",
    "                   [0,1,1,1,0],\n",
    "                   [0,0,1,0,0]]).astype(float)\n",
    "    unocc_aperture = fftconvolve(unocculted_image,kernel,mode='valid')\n",
    "    norm = np.max(unocc_aperture)\n",
    "\n",
    "    # Convolve reference-subtract and raw frames with the aperture as well\n",
    "    occ_aperture = fftconvolve(occulted_image,kernel,mode='valid')\n",
    "\n",
    "    # Compute radial distance from center (in pixels)\n",
    "    indices = np.indices(unocc_aperture.shape)\n",
    "    center = np.array(unocc_aperture.shape) / 2.\n",
    "    radial = np.sqrt( (indices[0] - center[0])**2 + (indices[1] - center[1])**2 )\n",
    "    # Compute 20 annuli (uniform in radius)\n",
    "    radial_bins = np.linspace(0,np.max(radial),num=n_annuli)\n",
    "    annuli_inds = np.digitize(radial,radial_bins)\n",
    "\n",
    "    # Take the variance of raw and reference-subtracted images in each annulus and normalize by unocculted max\n",
    "    contrast = np.array([np.std(occ_aperture[annuli_inds == a]) for a in np.unique(annuli_inds)]) / norm    \n",
    "    return radial_bins, contrast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll iterate over exposure parameters and calculate the raw and reference-subtracted contrast at each iteration.\n",
    "\n",
    "NB: In the loop below, you'll notice we've replaced the call ```engine.calculate_batch([occulted,unocculted,reference])``` with three separate calls to ```engine.perform_calculation```. This is motivated by an [LRU cache](https://en.wikipedia.org/wiki/Cache_replacement_policies#Least_Recently_Used_.28LRU.29) implemented to cache the results of calls to WebbPSF to avoid duplicating unecessary calculations. However, this caching implementation is incompatible with multiprocessing, so to take advantage of it, we've serialized the calls to the Pandeia engine. This translates into a substantial times savings for sequential calculations that modify things like exposure time but not source properties or filters (which would require new calculations in WebbPSF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options.save_options()\n",
    "options.on_the_fly_PSFs = True\n",
    "options.wave_sampling = 11\n",
    "n_annuli = 30\n",
    "ngroup_list = range(2,11)\n",
    "\n",
    "raw_contrast_list = []\n",
    "refsub_contrast_list = []\n",
    "for i, ngroup in enumerate(ngroup_list):\n",
    "    print(\"Starting {} of {}\".format(i+1, len(ngroup_list)))\n",
    "    # Exposure parameters\n",
    "    numberofgroups = ngroup\n",
    "    numberofints = 1\n",
    "    occulted['configuration']['detector']['ngroup'] = numberofgroups\n",
    "    occulted['configuration']['detector']['nint'] = numberofints\n",
    "    unocculted['configuration']['detector']['ngroup'] = 200\n",
    "    unocculted['configuration']['detector']['nint'] = 10\n",
    "\n",
    "    # Pandeia calculation\n",
    "    # Here we avoid calculate_batch to take advantage\n",
    "    # of LRU caching, which is incompatible with\n",
    "    # multiprocessing.\n",
    "    result = pancake.engine.calculate_all(occulted)\n",
    "    options.set_saturation(False)\n",
    "    result_unocc = pancake.engine.calculate_target(unocculted)\n",
    "    options.set_saturation(True)\n",
    "    occ_slope = result['target']['2d']['detector']\n",
    "    unocc_slope = result_unocc['2d']['detector']\n",
    "    ref_slope = result['reference']['2d']['detector']\n",
    "\n",
    "    # PSF subtraction assuming photon noise, the normalization is done properly \n",
    "    centered_occ = occ_slope - np.nanmean(occ_slope)\n",
    "    centered_ref = ref_slope - np.nanmean(ref_slope)\n",
    "    reg_ref = pancake.analysis.register_to_target(centered_ref,centered_occ)\n",
    "    ref_sub = centered_occ - reg_ref\n",
    "\n",
    "    radial_bins, raw_contrast = quick_contrast(occ_slope,unocc_slope,n_annuli=n_annuli)\n",
    "    radial_bins, refsub_contrast = quick_contrast(ref_sub,unocc_slope,n_annuli=n_annuli)\n",
    "\n",
    "    raw_contrast_list.append(raw_contrast)\n",
    "    refsub_contrast_list.append(refsub_contrast)\n",
    "    print(\"Finishing {} of {}\".format(i+1, len(ngroup_list)))\n",
    "\n",
    "options.restore_options()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pix_scale = 0.11 # see https://jwst-docs.stsci.edu/display/JTI/MIRI+Coronagraphic+Imaging\n",
    "mask_radii = 0.33\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "for i,contrast in enumerate(raw_contrast_list):\n",
    "    plt.semilogy(radial_bins * pix_scale,contrast,label='Nints=1, Ngroups={}'.format(ngroup_list[i]))\n",
    "    plt.fill_between([0,mask_radii],1e-6,20,alpha=0.2,lw=0,color=[0.5,0.5,0.5])\n",
    "    plt.xlabel('Radial Separation (arcsec)')\n",
    "    plt.ylabel('Contrast')\n",
    "    plt.xlim(0,5.6)\n",
    "    plt.ylim(1e-5,1e-2)\n",
    "    plt.title('Raw Target, {} V={}, Reference {}, V={}'.format(target_Sp,target_mV,ref_Sp,ref_mV),fontsize=15)\n",
    "    plt.legend(loc='upper right')\n",
    "    \n",
    "plt.figure(figsize=(10,5))\n",
    "for i,contrast in enumerate(refsub_contrast_list):\n",
    "    plt.semilogy(radial_bins * pix_scale,contrast,label='Nints=1, Ngroups={}'.format(ngroup_list[i]))\n",
    "    plt.fill_between([0,mask_radii],1e-6,20,alpha=0.2,lw=0,color=[0.5,0.5,0.5])\n",
    "    plt.xlabel('Radial Separation (arcsec)')\n",
    "    plt.ylabel('Contrast')\n",
    "    plt.xlim(0,5.6)\n",
    "    plt.ylim(1e-5,1e-1)\n",
    "    plt.title('Reduced Target, {} V={}, Reference {}, V={}'.format(target_Sp,target_mV,ref_Sp,ref_mV),fontsize=15)\n",
    "    plt.legend(loc='upper right')\n",
    "    "
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}