NIRC2_Preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#NIRC2 preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We present here functions dedicated to Keck NIRC2 preprocessing. It consists in a basic preprocessing (masterflat), determination of the center of the VORTEX from sky or sci images, registration of the images and the creation of cube.\n",
    "\n",
    "\n",
    "More specifically, this notebook talks about:\n",
    "0. [Import](#import)\n",
    "1. [Basic preprocessing](#prepro)\n",
    "  0. [Open a fits file](#open)  \n",
    "  1. [Deal with header](#head)\n",
    "  2. [Master flat](#mflat)\n",
    "  3. [Preprocess files](#prefiles)\n",
    "  4. [Remove bad pixels](#badpix)\n",
    "  5. [Create a cube from fits images](#cube)\n",
    "2. [Find the VORTEX center from a sky image](#center)\n",
    "   0. [Introduction](#intro)\n",
    "   1. [Step-by-step procedure](#proc)\n",
    "       1. [Initialization](#ini)\n",
    "       2. [Define a model (optional)](#usermodel)\n",
    "       3. [Minimization](#min)\n",
    "       4. [Representation](#res)\n",
    "   4. [Routine](#routine)\n",
    "   5. [Evolution of the VORTEX position](#evo)\n",
    "3. [Find the VORTEX center directly from sci images](#centerdust)   \n",
    "4. [Registration](#regi)\n",
    "5. [Parallactic angles](#paral)\n",
    "    0. [PA for a set of images](#paral_set)\n",
    "    1. [PA for a single image](#paral_single)\n",
    "    2. [Calculating the PA directly from $\\mathrm{DEC}$, $\\phi$ and $H$](#paral_calc)\n",
    "6. [Crop the cube](#crop)\n",
    "7. [ADI PCA and RDI PCA: brief overview](#pca)\n",
    "\n",
    "**Tips**\n",
    "+ To obtain the complete docstring for a specific function, open a new cell and run: function_name?.\n",
    "+ To take a look at the source code of a specific function, open a new cell and run: function_name??.\n",
    "+ Cells which include both parameter(s) and code are splitted in two parts, respectively #Parameters and #Let's go. So, you just have to modify the #Parameters part."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Import <a id='import'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import all the functions from the module **NIRC2_Preprocessing**. Furthermore, we will make use of few VIP and numpy package functions. Whatever part of the notebook you plan to execute, you better first run this cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------\n",
      " oooooo     oooo ooooo ooooooooo.   \n",
      "  `888.     .8'  `888' `888   `Y88. \n",
      "   `888.   .8'    888   888   .d88' \n",
      "    `888. .8'     888   888ooo88P'  \n",
      "     `888.8'      888   888         \n",
      "      `888'       888   888         \n",
      "       `8'       o888o o888o        \n",
      "------------------------------------\n",
      "  Vortex Image Processing pipeline  \n",
      "------------------------------------\n"
     ]
    }
   ],
   "source": [
    "from NIRC2_Preprocessing import *\n",
    "import numpy as np\n",
    "from vip.fits import display_array_ds9, write_fits\n",
    "#from vip.calib.cosmetics_sdi import bad_pixel_removal\n",
    "from vip.calib.badpixremoval import bp_clump_removal, bp_annuli_removal, bp_isolated_removal\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Basic preprocessing <a id='prepro'></a> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Open a fits file <a id='open'></a>\n",
    "\n",
    "**Summary**: Open a fits file and display it with DS9. Optionally, the header is extracted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily open a fits image, even those which raise a *missing END card* error, by using the function **open_fits( )**. If the parameter **header** is set to *False* (default value), only 1 output is returned: image = open_fits(path); otherwise, 2 outputs are returned, respectively the image and the header. More information about how to handle **header** can be found [here](#head)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/n0022.fits'\n",
    "\n",
    "# Let's go\n",
    "image, header = open_fits(path, header=True, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "open_fits??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable **image** contains the fits image as a numpy array and can be displayed with DS9 using the VIP function **display_array_ds9( )**. As explained in the VIP docstring, DS9 should be already installed in your system along with XPA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print (image.shape)\n",
    "display_array_ds9(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can also open a list of fits files with DS9. To this aim, we list all the fits files included in a repository using the function **listing( )**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/'\n",
    "\n",
    "# Let's go\n",
    "file_list = listing(path, ext='fits')\n",
    "cube = [open_fits(f, header=False, verbose=False) for f in file_list]\n",
    "display_array_ds9(cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Deal with header <a id='head'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable **header** is a pyfits.header.Header object (which behaves like a dictionary object) which contains all the fits header cards. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/sci/flatted/N2.20091101.23009_flatted.fits'\n",
    "cards = ['OBSERVER','TELESCOP','CAMNAME','FILTER','OBJECT',\n",
    "         'RA','DEC','DATE-OBS','UTC','FILENAME',\n",
    "         'ITIME','COADDS','AIRMASS',]\n",
    "\n",
    "# Let's go\n",
    "_, header = open_fits(path, header=True, verbose=False)\n",
    "for card in cards:\n",
    "    print ('{} = {}'.format(card,header[card]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To display the full header with comments, just run the next cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function **extract_headers( )** allows to extract all the common header items from a list of fits images. It can be useful if, for instance, you want to check whether all the files of interest were taken with the same filter, ...  \n",
    "This function returns a dictionary where the keys are the header cards and the values are lists of the corresponding header values. Let us note that the last list entry for each key is a string which gives the corresponding header card definition (as it is in the original header)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path_flat = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/'\n",
    "cards = ['FILTER','ITIME','COADDS']\n",
    "short = True # if True, only the filename is displayed hereafter. Otherwise, the whole path is displayed.\n",
    "\n",
    "# Let's go\n",
    "file_list = listing(path_flat, ext='fits')\n",
    "h = extract_headers(file_list)\n",
    "for card in cards:\n",
    "    print ('{}: {}'.format(card,h[card][-1]))\n",
    "    for k,f in enumerate(file_list):\n",
    "        print ('{} = {}'.format(f[-f[::-1].find('/')*short:], h[card][k]))\n",
    "    print ('')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function **find_header_card( )** allows to search if there exists any header card which:\n",
    "+ contains a given substring\n",
    "+ starts with a given substring\n",
    "+ ends with a given substring\n",
    "\n",
    "It may be useful if you don't know the exact header card but you suspect that it should contain a given substring. For instance, if you want to obtain the total exposure time without knowing the corresponding header card, you should do:\n",
    "\n",
    "\\>>> result, info = find_header_card(header,'TIME',criterion='find', info=True)  \n",
    "\\>>> for key, value in result.items(): print '{} = {} /{}'.format(key, value, info[key])  \n",
    "WXTIME = 01:51:45.00 /UT of logged weather keyword values  \n",
    "ELAPTIME = 30.0 /ITIME * COADDS (added by KOA)  \n",
    "ITIME = 15.0 /Integration time per coadd  \n",
    "GUIDTIME = null /UT of logged GUIDFWHM value  \n",
    "\n",
    "Now, you know that the corresponding header card is *ELAPTIME*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/sci/N2.20091101.23009.fits'\n",
    "card = 'TIME'\n",
    "\n",
    "# Let's go\n",
    "_, header = open_fits(path, header=True, verbose=False)\n",
    "result, info = find_header_card(header,card,criterion='find', info=True)\n",
    "for key, value in result.items(): print ('{} = {} /{}'.format(key, value, info[key]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Master flat <a id='mflat'></a>\n",
    "\n",
    "**Summary**: create a master flat from a set of flat images. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the repository which contains all flats and list them. We list all the fits files using the function **listing( )**. If the parameter **selection** is set to *True*, each fits image will be opened with DS9 and you will be asked to keep it or discard it. If *False*, all fits images in the folder will be used to create the master flat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0002.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0003.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0004.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0005.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0006.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0007.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0008.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0009.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0010.fits\n",
      "/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/n0011.fits\n"
     ]
    }
   ],
   "source": [
    "# Parameters\n",
    "path_flat = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/flats/'\n",
    "selection = False\n",
    "\n",
    "# Let's go\n",
    "file_list = listing(path_flat, selection=selection, ext='fits')\n",
    "print ('')\n",
    "for filename in file_list:\n",
    "    print (filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the flat images listed in **file_list**, the master flat is created with the function **masterFlat( )**.\n",
    "\n",
    "+ If **header** is *True*, a dict which contains all individual header informations associated to all the flat images is returned.\n",
    "+ If **norm** is *True*, the master flat is normalized.\n",
    "+ If **display** is *True*, the master flat is automatically displayed with DS9.\n",
    "+ If **save** is *True*, the master flat is saved into the 'mflat.fits' file in the same folder as the individual flats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mflat, headers = masterFlat(file_list, header=True, norm=True, display=True, save=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable **headers** is similar to the one returned by **extract_headers( )** (see [here](#head)).\n",
    "Therefore, one can retrieve informations about all the flats as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "cards = ['DATE-OBS','EXPSTART']\n",
    "\n",
    "# Let's go\n",
    "for card in cards:\n",
    "    print ('{}: {}'.format(card,headers[card][-1]))\n",
    "    for k,f in enumerate(file_list):        \n",
    "        print ('{} = {}'.format(f, headers[card][k]))\n",
    "    print ('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Preprocess files <a id='prefiles'></a>\n",
    "**Summary**: divide a set of fits image by the master flat and save them into a folder. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define two paths which respectively contains: \n",
    "+ all the files to process\n",
    "+ the master flat \n",
    "\n",
    "If you want to process one file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "path_mflat = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/flat/mflat.fits' # Master flat\n",
    "file_list = ['/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/sci/N2.20091101.23009.fits']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... but if you want to process a list of files, then we list them using the function **listing( )**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n0022.fits\n",
      "n0023.fits\n",
      "n0024.fits\n",
      "n0025.fits\n",
      "n0026.fits\n",
      "n0027.fits\n",
      "n0028.fits\n",
      "n0029.fits\n",
      "n0030.fits\n",
      "n0031.fits\n",
      "n0032.fits\n",
      "n0033.fits\n",
      "n0034.fits\n",
      "n0035.fits\n",
      "n0036.fits\n",
      "n0037.fits\n",
      "n0038.fits\n",
      "n0039.fits\n",
      "n0040.fits\n",
      "n0041.fits\n",
      "n0042.fits\n",
      "n0043.fits\n",
      "n0044.fits\n",
      "n0045.fits\n",
      "n0046.fits\n",
      "n0047.fits\n",
      "n0048.fits\n",
      "n0049.fits\n",
      "n0050.fits\n",
      "n0051.fits\n",
      "n0052.fits\n",
      "n0053.fits\n",
      "n0054.fits\n",
      "n0055.fits\n",
      "n0056.fits\n",
      "n0057.fits\n",
      "n0058.fits\n",
      "n0059.fits\n",
      "n0060.fits\n",
      "n0061.fits\n",
      "n0062.fits\n",
      "n0063.fits\n",
      "n0064.fits\n",
      "n0065.fits\n",
      "n0066.fits\n",
      "n0067.fits\n",
      "n0068.fits\n",
      "n0069.fits\n",
      "n0070.fits\n",
      "n0071.fits\n",
      "n0072.fits\n",
      "n0073.fits\n",
      "n0074.fits\n",
      "n0075.fits\n",
      "n0076.fits\n",
      "n0077.fits\n",
      "n0078.fits\n",
      "n0079.fits\n",
      "n0080.fits\n",
      "n0081.fits\n",
      "n0082.fits\n",
      "n0083.fits\n",
      "n0084.fits\n",
      "n0085.fits\n",
      "n0086.fits\n",
      "n0087.fits\n",
      "n0088.fits\n",
      "n0089.fits\n",
      "n0090.fits\n",
      "n0091.fits\n",
      "n0092.fits\n",
      "n0093.fits\n",
      "n0094.fits\n",
      "n0095.fits\n",
      "n0096.fits\n",
      "n0097.fits\n",
      "n0098.fits\n",
      "n0099.fits\n"
     ]
    }
   ],
   "source": [
    "# Parameters\n",
    "path_mflat = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/calibration/mflat.fits' # Master flat\n",
    "path_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/'\n",
    "short = True # if True, only the filename is displayed hereafter. Otherwise, the whole path is displayed.\n",
    "\n",
    "# Let's go\n",
    "file_list = listing(path_files, selection=False, ext='fits')\n",
    "for f in file_list:\n",
    "    print (f[-f[::-1].find('/')*short:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, by the use of the function **applyFlat( )**, we divide all images by the master flat and (optionally) save them into a specific folder. \n",
    "\n",
    "+ If **header** is *True*, a dict which contains all individual header informations associated to all the flat images is returned.\n",
    "+ If **display** is *True*, the preprocessed images are automatically displayed with DS9.\n",
    "+ If **save** is *True*, the preprocessed images are saved into the 'path_files/flatted/' repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------------------\n",
      "Starting time: 2015-10-19 10:39:19\n",
      "-------------------------------------------------------------------\n",
      "PREPROCESSING IMAGES\n",
      "\n",
      "Save = True\n",
      "\n",
      "Fits files successfully created\n",
      "\n",
      "-------------------------------------------------------------------\n",
      "Images succesfully preprocessed\n",
      "Running time:  0:00:05.013454\n",
      "-------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "preprocessed, headers = applyFlat(file_list, path_mflat, header=True, display=False, save=True, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable **preprocessed** is a array-like which contains the preprocessed images for a possible future use. We display a raw image and its associated preprocess one with DS9.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "display_array_ds9(preprocessed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Remove bad pixels <a id='badpix'></a>\n",
    "\n",
    "This section is dedicated to the image bad pixels removal. We use the vip function **bp_clump_removal( )**. For more informations, please consult the vip docstring (>>> bp_clump_removal?)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/processed/flatted/n0010_flatted.fits'\n",
    "\n",
    "# Let's go\n",
    "image, header = open_fits(path, header=True, verbose=False)\n",
    "display_array_ds9(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove bad pixels procedure on one signe image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "center = np.array([296,296]) # /Define the approximate y and x coordinates of the star where no bad pixel is corrected \n",
    "                             # /in a circle of 1.8 fwhm (given by fwhm_round). \n",
    "fwhm = 30\n",
    "sig = 5.0\n",
    "\n",
    "\n",
    "res = bp_clump_removal(image, center[1], center[0], fwhm=fwhm, sig=sig, verbose=True)\n",
    "#res = bp_annuli_removal(frame, center[1], center[0], fwhm=fwhm, sig=sig, verbose=True)\n",
    "\n",
    "display_array_ds9(cube, res)\n",
    "\n",
    "print ('')\n",
    "print ('DONE !')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "display_array_ds9(cube, cube_cu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Create a cube from fits images <a id='cube'></a>\n",
    "**Summary**: create and save a cube from a set of fits images.\n",
    "\n",
    "From a set of $N$ fits images (registered or not), we create a cube with the shape $N \\times l \\times c$ where $l \\times c$ corresponds to the size of each image in pixels. \n",
    "\n",
    "+ If **header** is *True*, a list with all fits image headers is returned.\n",
    "+ If **save** is *True*, the cube is saved into the current repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20150609/sci/sci_Tint_0p2_coadds_100/flatted/'\n",
    "\n",
    "# Let's go\n",
    "file_list = listing(path_files, selection=False)\n",
    "cube, headers = create_cube_from_frames(path_files, header=True, verbose=False, save=False)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can display the cube with DS9... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "display_array_ds9(cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##Find the VORTEX center from a sky image <a id='center'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "0. [Introduction](#intro)\n",
    "1. [Step-by-step procedure](#proc)\n",
    "   1. [Initialization](#ini)\n",
    "   2. [Define a model (optional)](#usermodel)\n",
    "   3. [Minimization](#min)\n",
    "   4. [Representation](#res)\n",
    "4. [Routine](#routine)\n",
    "5. [Evolution of the VORTEX position](#evo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Introduction <a id='intro'></a> \n",
    "We model the VORTEX signature profile with a 2-D Gaussian profile. It constitutes so far a good approximation and should be sufficient to obtain the position of the center and reduce the effect due to the model choice. For this purpose, we minimize a function of merit which compare the intensity of all pixels in a squared box with the corresponding intensity obtained from the model. The function of merit is a reduced $\\chi^2_r$ defined by:\n",
    "\n",
    "$\\chi^2_r = \\frac{1}{N-n_p}\\sum_{j=1}^{N} \\frac{\\left(I_j - I_j^{model}\\right)^2}{I_j}$,\n",
    "\n",
    "where $n_p$ is the number of model parameters and $N$ the total number of pixels in the box (if $L$ is the size of the box in pixels, we have $N = L^2$). The error associated to $I_j$ is $\\sigma_j = \\sqrt{I_j}$. The 2-D Gaussian profile is defined by:\n",
    "\n",
    "\n",
    "$I^{Gaussian}(x,y) = \\mbox{bkg} + I_0 \\exp{\\left[- \\left(\\frac{(x-x_0)^2}{2 \\sigma_x^2} + \\frac{(y-y_0)^2}{2 \\sigma_y^2}\\right)\\right]}$,\n",
    "\n",
    "where $\\mbox{bkg}$ represents the background, $I_0$ the maximum intensity (taking account of the $\\mbox{bkg}$) and $(x_0,y_0)$ the position where $I = I_0$. Other models are available, such as a cone or a Moffat profile defined by: \n",
    "\n",
    "\n",
    "$I^{Moffat}(x,y) = \\mbox{bkg} + I_0 \\left[1 + \\left(\\frac{(x-x_0)^2 + (y-y_0)^2}{\\alpha^2}\\right) \\right]^{-\\beta}$,\n",
    "where $\\alpha$ is a scale parameter and $\\beta$ the parameter which determines the overall shape of the profile."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Step-by-step procedure <a id='proc'></a> \n",
    "\n",
    "The following steps illustrate the procedure to determine the VORTEX center from a sky image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####Initialization <a id='ini'></a> \n",
    "We first display the image with DS9 in order to roughtly estimage the position of the VORTEX center."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "path_file = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20150609/sky/flatted/n0373_flatted.fits'\n",
    "# Let's go\n",
    "image = open_fits(path_file, header=False, verbose=False)\n",
    "display_array_ds9(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize the center position of a squared box and its size (in pixels). Then, we display a 3-D surface plot of the pixel intensities in the box by using the function **plot_surface** (for more information about this function: ?plot_surface in a new cell). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initializate the center and the size (in pixels) of the squared box\n",
    "center, size = (518,730), 100 # VORTEX center \n",
    "# (521,731), 16 # VORTEX center\n",
    "# (592,606), 31 # Dust\n",
    "\n",
    "# Display a 3-D surface plot of the box\n",
    "plot_surface(image, center, size, output=False, figsize=(16,14), cmap='jet')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We adopt a model for the VORTEX signature and define the associated additional parameters. This module already includes various models:\n",
    "+ Gaussian profile, callable through **gauss2d( )**, 2 additional parameters = $\\sigma_x$, $\\sigma_y$\n",
    "+ Gaussian symetrical profile, callable through **gauss2d_sym( )**, 1 additional parameter = $\\sigma$\n",
    "+ Moffat profile, callable through **moffat( )**, 2 additional parameters = $\\alpha$, $\\beta$\n",
    "+ Cone profile, callable through **cone( )**, 1 additional parameter = radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fun = gauss2d_sym\n",
    "p_additional = [12]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the vector which contains all the initial parameter values is created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# The box\n",
    "box = image[center[0]-size//2:center[0]+size//2,center[1]-size//2:center[1]+size//2]\n",
    "\n",
    "# The center of the box\n",
    "x_ini, y_ini = (size//2,size//2) \n",
    "\n",
    "# The background estimation\n",
    "bkg_ini = np.median(image) \n",
    "\n",
    "# The maximum intensity\n",
    "i0_ini = np.max(box) - bkg_ini\n",
    "\n",
    "# The initial values for all parameters to optimize\n",
    "p_initial = np.array([x_ini,y_ini,i0_ini,bkg_ini]+p_additional)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####Define a model (optional) <a id='usermodel'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have already adopted a model for the VORTEX signature, just skip this step and go to [Minimization](#min). \n",
    "\n",
    "Let us note that you can also define your own model. For instance, we adopt a Sersic profile defined by:\n",
    "\n",
    "$I^{Sersic}(x,y) = \\mbox{bkg} + I_0 \\exp{\\left[- \\left(\\frac{\\sqrt{(x-x_0)^2 + (y-y_0)^2}}{\\alpha}\\right)^{1/n}\\right]}$\n",
    "\n",
    "Then, we only have to define a function **sersic( )**, and the corresponding **p_initial** vector. The only requirement is that the 2 first arguments must the (x,y) grid. Indeed, these arguments are fixed the others will be optimized during the minimization procedure. However, we advocate to organize the arguments as follow: $x, y, x_0, y_0, I_0, \\mbox{bkg}$, [additional parameters] where \n",
    "+ $(x_0,y_0)$ locates the position of the maximum intensity\n",
    "+ $I_0$ is the maximum intensity irrespectively of the background, i.e. $I_0 = I(x_0,y_0) - \\mbox{bkg}$\n",
    "+ $\\mbox{bkg}$ is the background\n",
    "\n",
    "For the case of the Sersic profile, there are 2 additional parameters, respectively $\\alpha$ and $n$. Finally, the code is given by: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Define the Sersic profile as a callable function\n",
    "def sersic(x, y, x0, y0, i0, bkg, alpha, n):\n",
    "    r = ((x-x0)**2+(y-y0)**2)**0.5\n",
    "    return bkg + i0 * np.exp(-(r/(alpha))**(1/float(n)))\n",
    "\n",
    "# *fun* is passed to vortex_center() function. But, we can directly pass sersic or any other.\n",
    "fun = sersic\n",
    "p_additional = [20, 0.7]\n",
    "\n",
    "# Define the vector which contains the initial values for all parameters, i.e. x0, y0, i0, bkg, [additional parameters]\n",
    "box = image[center[0]-size//2:center[0]+size//2,center[1]-size//2:center[1]+size//2]\n",
    "x_ini, y_ini = (size//2,size//2) \n",
    "bkg_ini = np.median(image)\n",
    "i0_ini = np.max(box) - bkg_ini\n",
    "\n",
    "p_initial = np.array([x_ini,y_ini,i0_ini,bkg_ini]+p_additional)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "####Minimization <a id='min'></a> \n",
    "\n",
    "To start the minimization, we call the function **vortex_center( )**. The first 5 arguments are mandatory while the others are optional, including:\n",
    "\n",
    "+ If **Display** is *True*, some figures are displayed during the minimization.\n",
    "+ If **verbose** is *True*, the results are displayed at the end of the minimization.\n",
    "+ ****kwargs** which are options passed to scipy.optimize.minimize( ). For instance, solver_options = {'xtol': ...}.\n",
    "\n",
    "The **vortex_center( )** function returns a tuple of 3 objects:\n",
    "1. the position of the VORTEX center in the original image\n",
    "2. the information returned by the minimization tool\n",
    "3. the box grid as a tuple of 2 numpy.array (useful to represent the model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "solver_options = {'xtol': 1e-04, 'maxiter': 1e+05, 'maxfev': 1e+05}\n",
    "center_vor, minimization_output, grid = vortex_center(image, \n",
    "                                                      center, \n",
    "                                                      size, \n",
    "                                                      p_initial,\n",
    "                                                      fun,\n",
    "                                                      display= True, \n",
    "                                                      verbose=True,\n",
    "                                                      savefig=False,\n",
    "                                                      method = 'Nelder-Mead',\n",
    "                                                      options = solver_options)\n",
    "\n",
    "p_optimized = minimization_output.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####Representation <a id='res'></a> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print ('Representation of the best model (in terms of chi2)')\n",
    "print ('---------------------------------------------------')\n",
    "print ('Model adopted: {}'.format(fun))\n",
    "\n",
    "model = fun(grid[0],grid[1],*p_optimized)\n",
    "plot_surface(model, figsize=(12,10), cmap='jet')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "###Routine <a id='routine'></a> \n",
    "\n",
    "Here is a routine which allows to determine the center of the VORTEX for a set of raw images. Into a loop, each image is optionally preprocessed (flat) and the position is determined (Nelder-Mead minimization). \n",
    "\n",
    "Parameters initialization. The **cards** parameter allows to extract specific header cards from all the images. It can be useful when you want to [determine the evolution of the position of the VORTEX center in function of time](#evo) or pointing position, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Repositories which contain all files to process (or a list of file paths) and, if required, the master flat.\n",
    "path_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/sci/flatted/reg/'\n",
    "path_mflat = None\n",
    "\n",
    "# Box parameters\n",
    "center, size = (512,512), 20\n",
    "# (521,731), 16 # VORTEX center\n",
    "# (592,606), 31 # Dust\n",
    "\n",
    "# Model and additional parameter(s), depending on which model we've adopted\n",
    "fun = gauss2d_sym\n",
    "p_model = [1.5]\n",
    "\n",
    "# Routine parameters\n",
    "preprocess = False\n",
    "cards = ['EXPSTART','DATE-OBS','RA','DEC']\n",
    "verbose = 2\n",
    "\n",
    "# Optimization algorithm options\n",
    "options = {'xtol':1e-02, 'maxiter':1e+04,'maxfev':1e+04}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "center_all, success_all, file_list, header_cards = vortex_center_routine(path_files, \n",
    "                                                                         center, \n",
    "                                                                         size, \n",
    "                                                                         fun,\n",
    "                                                                         preprocess=preprocess, \n",
    "                                                                         path_mflat=path_mflat,\n",
    "                                                                         additional_parameters=p_model, \n",
    "                                                                         verbose=verbose,\n",
    "                                                                         cards=cards,\n",
    "                                                                         options=options)\n",
    "\n",
    "print ('')\n",
    "print ('Convergence reached for all minimization ? {}'.format(success_all))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Informations exctracted from the headers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, filename in enumerate(file_list):\n",
    "    print ('{} = [{:.2f},{:.2f}]'.format(filename,header_cards['RA'][k], header_cards['DEC'][k]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "###Evolution of the VORTEX position <a id='evo'></a> \n",
    "\n",
    "From the routine outputs, we can represent the position of the VORTEX center as a function of time. For this aim, we use the function **timeExtract( )** as well as the matplotlib.pyplot module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# From the observation date and time, we determine the delta time between the first and all other observations\n",
    "t_start = timeExtract(header_cards['DATE-OBS'],header_cards['EXPSTART'])\n",
    "delta_time = [(t-t_start[0]).seconds/3600. for t in t_start]\n",
    "\n",
    "# Then, we illustrate the obtained VORTEX position as a function of time\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(14,7))\n",
    "plt.hold('on')\n",
    "plt.plot(delta_time,center_all[:,0]-center_all[0,0],'.r', markersize=14)\n",
    "plt.xlabel(r'$\\Delta t$ (hour)',fontsize=20)\n",
    "plt.ylabel(r'$\\Delta x$ (pixels)',fontsize=20)\n",
    "plt.title(path_files)\n",
    "#plt.ylim([-0.03,0.06])\n",
    "#plt.savefig(path_files+'DUST_delta_x_'+'Gaussian'+'_'+path_files.split('/')[:-1][-1]+'.pdf')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(14,7))\n",
    "plt.hold('on')\n",
    "plt.plot(delta_time,center_all[:,1]-center_all[0,1],'.b', markersize=14)\n",
    "plt.xlabel(r'$\\Delta t$ (hour)',fontsize=20)\n",
    "plt.ylabel(r'$\\Delta y$ (pixels)',fontsize=20)\n",
    "plt.title(path_files)\n",
    "#plt.ylim([-0.86,-0.74])\n",
    "#plt.savefig(path_files+'DUST_delta_y_'+'Gaussian'+'_'+path_files.split('/')[:-1][-1]+'.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Find the VORTEX center directly from sci images <a id='centerdust'></a>\n",
    "\n",
    "In the case of a dust signature can be seen on all sci images taken during the same run, we assume that its relative position with respect to the VORTEX center is constant. Then, if a sky image is available, we can determine the position of the VORTEX center for all the sci images by:\n",
    "+ determining the center of the dust signature for all sci images\n",
    "+ determining the VORTEX and dust centers in the sky image\n",
    "+ deducing the relative position of the dust with respect to the VORTEX center\n",
    "+ determining the position of the VORTEX center for all sci images\n",
    "\n",
    "This procedure can take advantage of what we have already illustrated so far. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center of the dust signature in all sci images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_sci_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/flatted/'\n",
    "center, size = (592,606), 31\n",
    "fun = gauss2d\n",
    "p_model = [5,5]\n",
    "center_all_dust, success_all, file_list = vortex_center_routine(path_sci_files, \n",
    "                                                             center, \n",
    "                                                             size, \n",
    "                                                             fun,\n",
    "                                                             preprocess=False, \n",
    "                                                             path_mflat=None,\n",
    "                                                             additional_parameters=p_model, \n",
    "                                                             verbose=True,\n",
    "                                                             cards=None)\n",
    "\n",
    "# Illustration of the dust signature in the first image of file_list\n",
    "image = open_fits(file_list[0], header=False, verbose=False)\n",
    "plot_surface(image, center, size, output=False, figsize=(8,6), cmap='jet')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Centers of the VORTEX and dust signature in a sky image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "path_sky_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sky/flatted/n0010_flatted.fits'\n",
    "sky_image = open_fits(path_sky_files, header=False, verbose=False)\n",
    "display_array_ds9(sky_image)\n",
    "\n",
    "# Dust center\n",
    "center_dust, size_dust = (592,606), 31\n",
    "fun_dust = gauss2d\n",
    "p_model_dust = [5,5]\n",
    "\n",
    "center_dust, success_dust, _ = vortex_center_routine(path_sky_files, \n",
    "                                                     center_dust, \n",
    "                                                     size_dust, \n",
    "                                                     fun_dust,\n",
    "                                                     preprocess=False, \n",
    "                                                     path_mflat=None,\n",
    "                                                     additional_parameters=p_model_dust, \n",
    "                                                     verbose=False,\n",
    "                                                     cards=None)\n",
    "\n",
    "# VORTEX center\n",
    "center_vortex, size_vortex = (520,730), 30\n",
    "fun_vortex = gauss2d\n",
    "p_model_vortex = [5,5]\n",
    "\n",
    "center_vortex, success_vortex, _ = vortex_center_routine(path_sky_files, \n",
    "                                                         center_vortex, \n",
    "                                                         size_vortex, \n",
    "                                                         fun_vortex,\n",
    "                                                         preprocess=False, \n",
    "                                                         path_mflat=None,\n",
    "                                                         additional_parameters=p_model_vortex, \n",
    "                                                         verbose=False,\n",
    "                                                         cards=None)\n",
    "\n",
    "print ('Center of the dust signature: {}'.format(center_dust))\n",
    "print ('Center of the VORTEX signature: {}'.format(center_vortex))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Relative position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "relative_position = center_vortex[0]-center_dust[0]\n",
    "print ('Relative position: {}'.format(relative_position))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The position of the VORTEX in all sci frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center_all = center_all_dust + np.tile(relative_position,(center_all_dust.shape[0],1))\n",
    "#center_all = center_all+np.random.uniform(-0.5,0.5,center_all.shape)\n",
    "print (center_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Registration<a id='regi'></a> \n",
    "\n",
    "Thanks to previous sections ([Find the VORTEX center](#center) or [Find the VORTEX center directly from sci images](#centerdust)), we are now able to determine the VORTEX center position. The next stage consists in registering a sequence of images with respect to the corresponding VORTEX center. To this end, we use the function **registration( )** which requires, at least, the 3 following parameters:\n",
    "+ **fileList** : the list of paths to the N sci images to register.\n",
    "+ **initial_position** : a N x 2 array which defines the N positions to register.\n",
    "+ **final_position** : a 1 x 2 array which defines the position where all initial_position are registered in the processed cube.\n",
    "\n",
    "We choose the **initial_position** to be the positions of the VORTEX center and the **final_position** to be the center of the frame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###VORTEX center <a id='regivor'></a>\n",
    "\n",
    "First, we have to determine the center of the VORTEX signature for a set of images. For this aim, select between [Routine](#routine) or [Find the VORTEX center directly from sci images](#centerdust) sections. At the end, we should have the variable **center_all**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Create the **initial_position** array <a id='regiini'></a>\n",
    "\n",
    "If you have determined the position of the center of the VORTEX from a dust signature on all sci images, the matrix **center_all** and the list **file_list** should respectively contain all the positions of the VORTEX center and all the sci image paths. Then, you can skip this step and go directly to the next sub-section [Registration](#regireg).\n",
    "\n",
    "If you have determined the position from few sky images, you have likely less VORTEX center positions than sci images. Then, you have to construct the **center_all** matrix in such a way that center_all.shape = (N,2) where N corresponds to the number of sci images to register. The next cell is dedicated to that part of the job. Some manual adjustments can be done here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Path: images to register\n",
    "path_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/flatted/' \n",
    "file_list = listing(path_files)\n",
    "#for filename in file_list:\n",
    "#    print filename\n",
    "\n",
    "# Initial position from the VORTEX center position\n",
    "part_0 = np.tile(center_all[0,:],(10,1))\n",
    "part_1 = np.tile(center_all[1,:],(len(file_list)-10,1)) # np.array([520.06,730.5])\n",
    "\n",
    "#part_0 = np.tile(center_from_sky[0,:],(9,1))\n",
    "#part_1 = np.tile(center_from_sky[1,:],(10+10,1))#(len(file_list)-9,1))\n",
    "#part_2 = np.tile(np.array([520.06,730.5]),(len(file_list)-9-10-10,1))\n",
    "\n",
    "# Concatenation\n",
    "center_all = np.concatenate((part_0,part_1))\n",
    "\n",
    " \n",
    "print ('Number of files: {},  center_all shape: {}'.format(len(file_list),center_all.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center_all = np.concatenate((np.tile([496,590],(5,1)),\n",
    "                             np.tile([220,590],(4,1)),\n",
    "                             np.tile([494,592],(1,1)),\n",
    "                             np.tile([200,477],(5,1)),\n",
    "                             np.tile([550,478],(5,1))))\n",
    "\n",
    "center_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Registration <a id='regireg'></a>\n",
    "\n",
    "After defining the **final_position**, we run the **registration( )** function. It returns the cube which contains all registered images and (optional) their headers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "path_files = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/GJ504/sci/flatted/' \n",
    "file_list = listing(path_files)\n",
    "\n",
    "# Let's define the pixel coordinates at which all the VORTEX center will be shifted\n",
    "final_position = np.array([512,512])\n",
    "\n",
    "# Registration\n",
    "cube_reg, headers = registration(file_list, \n",
    "                                 initial_position = center_all, \n",
    "                                 final_position = final_position, \n",
    "                                 header=True, \n",
    "                                 verbose=False, \n",
    "                                 display=True,\n",
    "                                 save=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Parallactic angles <a id='paral'></a> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###PAs of a set of images <a id='paral_set'></a>\n",
    "\n",
    "For a set of images, the parallactic angles can be retreived from few header cards:\n",
    "+ *ROTPPOSN*: the rotator physical position\n",
    "+ *PARANTEL*: the parallactic angle telescope\n",
    "+ *EL*: the telescope elevation\n",
    "+ *INSTANGL*: porg to instrument angle\n",
    "\n",
    "combined as follows (Justin Crepp):\n",
    "\n",
    "parallactic angle = *ROTPPOSN* + *PARANTEL* - *EL* - *INSTANGL* ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Define the folder which contains the images\n",
    "root_to_images = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/sci/flatted/reg/'\n",
    "\n",
    "# Let's go\n",
    "file_list = listing(root_to_images, selection=False, ext='fits')\n",
    "parallactic_angles = np.zeros(len(file_list))\n",
    "\n",
    "for k, filename in enumerate(file_list):\n",
    "    _, header = open_fits(filename, header=True)\n",
    "    parallactic_angles[k] = header['ROTPPOSN']+header['PARANTEL']-header['EL']-header['INSTANGL']    \n",
    "    \n",
    "print (parallactic_angles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We save the parallactic angles into a 1-column fits file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pa_filename = '/Users/Olivier/Documents/ULg/VORTEX/Data/Cube_PSF_PA/Keck/HR8799/20091101/pa_HR8799_20091101.fits' \n",
    "write_fits(pa_filename,np.array(parallactic_angles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###PA of a single image <a id='paral_single'></a>\n",
    "\n",
    "The same idea but for only one image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20150609/sci/sci_Tint_0p2_coadds_100/n0377.fits'\n",
    "_, header = open_fits(data, header=True)\n",
    "parallactic_angle = header['ROTPPOSN']+header['PARANTEL']-header['EL']-header['INSTANGL']\n",
    "\n",
    "print ('PA = {}'.format(parallactic_angle))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Calculating the PA directly from $\\delta$, $\\phi$ and $H$. <a id='paral_calc'></a>\n",
    "\n",
    "From the header, one can retreived informations about the declinaison $\\delta$, observation date and hour angle $H$ at the given observation. Given the geographical latitude $\\phi$ of the observer, one can deduce the parallactic angle $q$ from (Meeus, 1998, Eq. 13.1 p94):\n",
    "\n",
    "$\\tan{q} = \\frac{\\sin{H}}{\\tan{\\phi} \\cos{\\delta} - \\sin{\\delta} \\cos{H}}$ .\n",
    "\n",
    "The function **get_parang( )** performs that calculation and requires:\n",
    "+ the header\n",
    "+ the geographical latitude of the observer\n",
    "+ the epoch of the declinaison written in the header (default=None). \n",
    "\n",
    "If **epoch** is None, the header card *EQUINOX* is used if it exists, otherwise 'FK5' is used. The parameter **epoch** can also accept 'FK5', 'todate' or any epoch (float). In the first and third cases, the declinaison is precessed from initial epoch to observation epoch. Furthermore, if **epoch** is not None, the header card *EQUINOX* is not used.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20150609/sci/sci_Tint_0p2_coadds_100/n0377.fits'\n",
    "_, header = open_fits(data, header=True)\n",
    "parallactic_angle = get_parang(header, latitude=19.82636, epoch=None)\n",
    "\n",
    "print ('PA = {}'.format(parallactic_angle))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Crop the cube<a id='crop'></a> \n",
    "\n",
    "This part is dedicated to crop the cube of registered frames. To this end, we use the function **cube_crop_frames_optimized( )** which requires at leat the 3 following parameters:\n",
    "+ cube : the cube to crop\n",
    "+ ceny : the y-coordinate of the center of the cropped cube\n",
    "+ cenx : the x-coordinate of the center of the cropped cube\n",
    "\n",
    "The optimized size of the cube frames are automatically determined in such a way it maximizes the area of non-zero pixel values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cube_path = '/Users/Olivier/Documents/ULg/VORTEX/Data/RAW/Keck/HR8799/20091101/sci/flatted/reg/reg/cube_reg.fits'\n",
    "cube_reg = open_fits(cube_path)\n",
    "\n",
    "cropped_cube_path = '/Users/Olivier/Documents/ULg/VORTEX/Data/Cube_PSF_PA/Keck/HR8799/20091101/cube_HR8799_20091101.fits'\n",
    "\n",
    "final_position = np.array([512,512])\n",
    "cube_reg_crop = cube_crop_frames_optimized(cube_reg, \n",
    "                                           final_position[1], final_position[0], \n",
    "                                           ds9_indexing=True, \n",
    "                                           verbose=True, \n",
    "                                           display=True,\n",
    "                                           save=True,\n",
    "                                           filename=cropped_cube_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## ADI PCA and RDI PCA: brief overview <a id='pca'></a>\n",
    "\n",
    "We propose here a basic call of VIP functions to perform ADI or RDI PCA of the obtained cube and have a quick look to the result. A more in-depth analyze can be found in the **Photo_Astrometry.ipynb** IPython notebook, including PCA, Vortex center determination, Speckle noise and MCMC, ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "path_cube = '/Users/Olivier/Documents/ULg/VORTEX/Data/Cube_PSF_PA/Keck/HR8799/20091101/cube_HR8799_20091101_CleanedUp.fits'\n",
    "path_ref = '/Users/Olivier/Documents/ULg/VORTEX/Data/Cube_PSF_PA/Keck/HR8799/20150609/cube_HD219196_20150609.fits'\n",
    "path_pa = '/Users/Olivier/Documents/ULg/VORTEX/Data/Cube_PSF_PA/Keck/HR8799/20091101/pa_HR8799_20091101.fits'\n",
    "\n",
    "cube = open_fits(path_cube, header=False, verbose=False)\n",
    "cube_ref = open_fits(path_ref, header=False, verbose=False)\n",
    "angs = open_fits(path_pa, header=False, verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the cube(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "display_array_ds9(cube)\n",
    "#display_array_ds9(cube, cube_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vip.pca import pca\n",
    "out_ADI = pca(cube, angs, cube_ref=None, ncomp=3, svd_mode='randsvd', full_output=False)\n",
    "#out_RDI = pca(cube, angs, cube_ref=cube_ref, ncomp=1, svd_mode='randsvd', full_output=False, center='global')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the result(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_array_ds9(out_ADI)\n",
    "#display_array_ds9(out_ADI,out_RDI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write the resulting image into a fits file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = '/Users/Olivier/Documents/ULg/VORTEX/Data/Cube_PSF_PA/Keck/HR8799/20091101/HR8799_ADI_20091101.fits'\n",
    "write_fits(filename,out_ADI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##Work in progress..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}