centroid_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:20:09.310351Z",
     "start_time": "2020-02-12T19:20:08.186641Z"
    },
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#imports\n",
    "from jwst_ta_cent_moments import checkbox, fine_centroid\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import CubicSpline\n",
    "import pandas as pd\n",
    "from astropy.io import fits\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:20:09.828681Z",
     "start_time": "2020-02-12T19:20:09.825543Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def open_im(file):\n",
    "    #import MIRI webbpsf file\n",
    "    \n",
    "    hdu = fits.open(file)\n",
    "    head = hdu[0].header\n",
    "    data = hdu[0].data\n",
    "    \n",
    "    return head, data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:20:11.682300Z",
     "start_time": "2020-02-12T19:20:11.667536Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def centroid_b(im, input_type='image', ext=0, cbox=5, cwin=5, incoord=(0., 0.), roi=None, bgcorr=-1, flat=None, flatext=0, out=None, thresh=0.05, silent=False):\n",
    "    \n",
    "    '''\n",
    "    Implementation of the JWST GENTALOCATE algorithm. Parameters key:\n",
    "    \n",
    "    - infile:       FITS filename\n",
    "    - input_type:   description of input data: 'image' or 'ramp'. If 'ramp'\n",
    "                    then make_ta_image functin is run. If 'image' (default)\n",
    "                    centroiding is performed directly on the data in the\n",
    "                    input file\n",
    "    - ext:          extension number of the FITS file containing the science data (default = 0)\n",
    "    - cbox:         the FULL size of the checkbox, in pixels, for coarse centroiding (default = 5)\n",
    "    - cwin:         the FULL size of the centroid window, in pixels, for fine centroiding (default = 5)\n",
    "    - incoord:      (x,y) input coordinates of the source position\n",
    "    - roi:          size of a region of interest to be used for the centroiding (optional). If not set, full image will be used for coarse                       centroiding. \n",
    "                        * setting an ROI also requires input coordinates\n",
    "                        * the ROI size must be bigger than the cbox parameter\n",
    "    - bgcorr:       background correction parameter. set to:\n",
    "                        * negative value for NO background subtraction (default)\n",
    "                        * 0 < bgcorr < 1 for fractional background subtraction\n",
    "                        * bgcorr > 1 for constant background subtraction number (this number will be subtracted from the entire image)\n",
    "    - out:          enter a filename for output of the fit results to a file (default = None)\n",
    "    - thresh:       the fit threshold, in pixels. default is 0.1 px. consider setting this to a higher number for testing, long-wavelength    \n",
    "                       data or low SNR data to prevent.\n",
    "    - silent:       set to True if you want to suppress verbose output\n",
    "    '''\n",
    "\n",
    "    # Read in data. Create the TA image if requested\n",
    "#     if input_type.lower() == 'image':\n",
    "#         hdu = fits.open(infile)\n",
    "#         im = hdu[ext].data\n",
    "#         h = hdu[ext].header\n",
    "#     elif input_type.lower() == 'ramp':\n",
    "#         im = make_ta_image(infile, ext=ext, useframes=3, save=False)\n",
    "#         # Save TA image for code testing\n",
    "#         h0 = fits.PrimaryHDU(im)\n",
    "#         hl = fits.HDUList([h0])\n",
    "#         indir, inf = os.path.split(infile)\n",
    "#         tafile = os.path.join(indir, 'TA_img_for_'+inf)\n",
    "#         hl.writeto(tafile, overwrite=True)\n",
    "\n",
    "    # Do background correction first\n",
    "    if bgcorr > 0.:\n",
    "        \n",
    "        # if a ROI size was provided, the value to be subtracted as background will be calculated using the pixels in the ROI only. Otherwise, use the full array.\n",
    "        if roi is not None:\n",
    "            im = bgrsub(im, bgcorr, roi, incoord, silent=silent)\n",
    "        else:\n",
    "            im = bgrsub(im, bgcorr, -1, incoord, silent=silent)\n",
    "\n",
    "    # Apply flat field\n",
    "    if flat is not None:\n",
    "        # Read in flat\n",
    "        with fits.open(flat) as ff:\n",
    "            flatfield = ff[flatext].data\n",
    "            \n",
    "        # Flat must be the same size as the data\n",
    "        ffshape = flatfield.shape\n",
    "        dshape = im.shape\n",
    "        if dshape != ffshape:\n",
    "            raise RunTimeError((\"WARNING: flat field shape ({}) does \"\n",
    "                                \"not match data shape ({})!\"\n",
    "                                .format(ffshape,dshape)))\n",
    "        # Apply flat\n",
    "        im = apply_flat_field(im, flatfield)\n",
    "        \n",
    "    ndim = np.ndim(im)\n",
    "    #pdb.set_trace()\n",
    "    \n",
    "    n = [np.size(im, axis=i) for i in range(ndim)]\n",
    "    \n",
    "    # NOTE: in python the x-coord is axis 1, y-coord is axis 0\n",
    "    xin = incoord[0]\n",
    "    yin = incoord[1]\n",
    "    if not silent:\n",
    "        print('Input coordinates = ({0}, {1})'.format(xin, yin))\n",
    "    \n",
    "    # Extract the ROI\n",
    "    if (roi is not None):\n",
    "        \n",
    "        #first check that the ROI is larger than the cbox size\n",
    "        assert roi > cbox, \"ROI size must be larger than the cbox parameter\"\n",
    "        \n",
    "        # now check that the ROI is a sensible number.\n",
    "        # if it's bigger than the size of the array, use the\n",
    "        # full array instead\n",
    "        if (roi >= n[0])|(roi >= n[1]):\n",
    "            print('ROI size is bigger than the image; using full image instead')\n",
    "            roi_im = im\n",
    "            xoffset = 0\n",
    "            yoffset = 0\n",
    "            #xc, yc = checkbox(im, cbox, bgcorr)\n",
    "        else:\n",
    "            roi_im = im[np.round(yin-(roi/2.)).astype(int):np.round(yin+(roi/2.)).astype(int),\n",
    "                        np.round(xin-(roi/2.)).astype(int):np.round(xin+(roi/2.)).astype(int)]\n",
    "            \n",
    "            \n",
    "            #print(\"ROI size is {0}\".format(np.shape(roi_im)))\n",
    "            xoffset = np.round(xin-(roi/2.)).astype(int)\n",
    "            yoffset = np.round(yin-(roi/2.)).astype(int)\n",
    "    else:\n",
    "        #xc, yc = checkbox(im, cbox, bgcorr)\n",
    "        roi_im = im\n",
    "        xoffset = 0\n",
    "        yoffset = 0\n",
    "    \n",
    "    # Perform coarse centroiding. Pay attention to coordinate\n",
    "    # offsets\n",
    "    xc, yc = checkbox(roi_im, cbox)\n",
    "    xc += xoffset\n",
    "    yc += yoffset\n",
    "    if not silent:\n",
    "        print('Coarse centroid found at ({0}, {1})'.format(xc, yc))\n",
    "    \n",
    "    # Iterate fine centroiding\n",
    "    # Take the threshold from the input parameter thresh\n",
    "    iter_thresh = thresh\n",
    "    nconv = 0\n",
    "    while nconv == 0:\n",
    "        xf, yf= fine_centroid(im, cwin, xc, yc)\n",
    "        err = np.sqrt((xf-xin)**2 + (yf-yin)**2)\n",
    "        if not silent:\n",
    "            print((\"Fine centroid found at (x, y) = ({0:.4f}, {1:.4f}). \"\n",
    "               \"Rms error = {2:.4f}\".format(xf, yf, err)))\n",
    "        if (abs(xf-xc) <= iter_thresh) & (abs(yf-yc) <= iter_thresh):\n",
    "            nconv = 1\n",
    "        xc = xf\n",
    "        yc = yf\n",
    "    \n",
    "\n",
    "    return xf, yf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:55:19.101273Z",
     "start_time": "2020-02-12T19:55:19.086100Z"
    }
   },
   "outputs": [],
   "source": [
    "head, data = open_im('MIRI_psf_model_F1000W_0.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:55:19.544686Z",
     "start_time": "2020-02-12T19:55:19.540915Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'F1000W'"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "head['FILTER']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:55:20.773565Z",
     "start_time": "2020-02-12T19:55:20.582994Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0xb1bbb9ef0>"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAD7CAYAAAD6gVj5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAa+0lEQVR4nO3da5Bc5X3n8e9Po1uA5SYZryzJllhNbMtODKwWSHxZrxUc4TgWWwtBxGuzWVWpUms2zq28YreCE8ovQiplnFRU3qgMCVZcgCM7YcpRrPUKnMuWV9FgMCAUhUFm0QA2CMlyMOgy0/99cZ4WrVb3zKPpnunTZ36fqlPd55yn+zw6M/PXcz2PIgIzsyqb0+sMmJlNNwc6M6s8BzozqzwHOjOrPAc6M6s8Bzozq7xpCXSS1knaL2lE0ubpuIaZWS51exydpAHgn4BrgFFgD3BTRDzZ1QuZmWWaOw3feSUwEhEHACTdB6wH2ga6+VoQCzl3GrJiZnX/zJFDEfGGqX7+Z//dufHy4fGstA8/dnxnRKyb6rW6bToC3VLgYMP+KHBVcyJJm4BNAAs5h6u0dhqyYmZ1/zu2/79OPv/y4XH+Yeebs9IOLHlqcSfX6rbpCHRqceyM+nFEbAW2Apyviz0PzazkAqhR63U2pmQ6At0osLxhfxnw/DRcx8xmUBCcjLyqa9lMR6DbAwxKWgk8B2wAfnEarmNmM8wluiQixiTdAuwEBoC7I2Jvt69jZjMrCMb79GlH01GiIyJ2ADum47vNrHdqZza394VpCXRmVj0BjDvQmVnVuURnZpUWwEm30ZlZlQXhqquZVVzAeH/GOQc6M8tTzIzoT34enZllEuOZW9a3TfI4N0kLJN2fzu+WtCId/6ikRxu2mqTLJrqWA52ZZSk6I5S1TSY9zm0LcC2wGrhJ0uqmZBuBIxGxCrgTuAMgIr4UEZdFxGXAx4BnIuLRia7nQGdmWYpxdF0r0Z16nFtEnADqj3NrtB64J73fDqyV1PzlNwH3TnYxt9GZWbZaRmktWSxpuGF/a3piUV3O49xOpUlTS48Ci4BDDWlu5MwAeQYHOjPLUi/RZToUEWsmOJ/zOLcJ00i6Cng1Ip6YLDMOdGaWJRDj3WvtynmcWz3NqKS5wAXA4YbzG8iotoLb6MzsLNRCWVuGU49zkzSfImgNNaUZAm5O768HHoy0yI2kOcANFG17k3KJzsyyBOJEDHTnu9o8zk3S7cBwRAwBdwHbJI1QlOQ2NHzF+4DR+to0k3GgM7MsxYDh7lUCWz3OLSJua3h/jKLU1uqz3wSuzr2WA52ZZTuLzohScaAzsywRYjz6s1nfgc7MstVcojOzKis6I/ozZPRnrs1sxnW7M2ImOdCZWbbx/ClgpeJAZ2ZZujwzYkY50JlZtpp7Xc2syopJ/Q50ZlZhgTjZpSlgM82BzsyyROABw2ZWdfKAYTOrtsAlOjObBdwZYWaVFmQ/VLN0HOjMLEux3GF/howpl0MlLZf0kKR9kvZK+mQ6frGkb0h6Kr1e1L3smlnvdHcB65nUSYV7DPiNiHg7xZM+P5EWoN0M7IqIQWBX2jezPhcUMyNytrKZco4i4oWI+HZ6/8/APop1GBsXnb0HuK7TTJpZOczGEt0pklYAlwO7gTdGxAtQBEPgkjaf2SRpWNLwSY53IxtmNo0i1NUSnaR1kvZLGpF0Rs1P0gJJ96fzu1OcqZ/7SUnfSs1mj0taONG1Om5ZlHQe8BXgVyPih1JeNE+rdm8FOF8XNy9ca2YlU3RGdGcKmKQBYAtwDcX6rXskDUXEkw3JNgJHImKVpA3AHcCNaY3XPwM+FhHfkbQIODnR9Toq0UmaRxHkvhQRX02Hvy9pSTq/BHixk2uYWVkUa0bkbBmuBEYi4kBEnKBYn3V9U5rGZrDtwFoVJakPAo9FxHcAIuLliBif6GKd9LqKYt3FfRHx2YZTjYvO3gw8MNVrmFl5FJ0RXVvAeilwsGF/NB1rmSYixoCjwCLgx4GQtFPStyV9arKLdVJ1fTfwMeBxSY+mY/8d+F3gy5I2As/SZl1GM+s/ZzEzYrGk4Yb9ram5qq5VNGxuwmqXZi7wHuDfAK8CuyQ9HBG72mVmyoEuIv6+TUYA1k71e82snM5yZsShiFgzwflRYHnD/jLg+TZpRlO73AXA4XT8byLiEICkHcAVFMPZWirfgBczK60ac7K2DHuAQUkrJc0HNlA0ezVqbAa7HngwIgLYCfykpHNSAPy3wJNMoD/nc5jZjIuAk7XulI0iYkzSLRRBawC4OyL2SrodGI6IIYo+gG2SRihKchvSZ49I+ixFsAxgR0T81UTXc6AzsyxF1bV7lcCI2AHsaDp2W8P7Y7Rp44+IP6MYYpLFgc7MspVx1kMOBzozy1IfXtKPHOjMLFN3q64zyYHOzLJ5zQgzq7Si19XLHZpZhflR6mY2K7jqamaV5l5XM5sV3OtqZpUWIcYc6Mys6lx1NbNKcxudmc0KDnRmVmkeR2dms4LH0ZlZpUXAWJcevDnTHOjMLJurrmZWaW6jM7NZIRzozKzq+rUzoj9bFs1sxkUUbXQ5Ww5J6yTtlzQiaXOL8wsk3Z/O75a0Ih1fIek1SY+m7X9Odi2X6MwskxjvUq+rpAFgC3ANxYLUeyQNRUTj+qwbgSMRsUrSBuAO4MZ07umIuCz3ei7RmVm2CGVtGa4ERiLiQEScAO4D1jelWQ/ck95vB9ZKmlLd2YHOzLLU57pmVl0XSxpu2DY1fd1S4GDD/mg61jJNRIwBR4FF6dxKSY9I+htJ750s7666mlmeKNrpMh2KiDUTnG9VMmv+9nZpXgDeHBEvS/rXwF9KekdE/LDdxVyiM7NsNZS1ZRgFljfsLwOeb5dG0lzgAuBwRByPiJcBIuJh4Gngxye6mAOdmWWJ1BmRs2XYAwxKWilpPrABGGpKMwTcnN5fDzwYESHpDakzA0mXAoPAgYku5qqrmWU7i6rrJN8TY5JuAXYCA8DdEbFX0u3AcEQMAXcB2ySNAIcpgiHA+4DbJY0B48AvR8Thia7XcaBLkXUYeC4iPixpJUUPysXAt4GPpV4VM+tz3ZwZERE7gB1Nx25reH8MuKHF574CfOVsrtWNqusngX0N+3cAd0bEIHCEYiyMmfW5iK4OL5lRHQU6ScuAnwO+kPYFfIBizAsUY2Cu6+QaZlYe3ZwZMZM6rbp+DvgU8C/S/iLgB2nMC7QeGwNAGlezCWAh53SYDTObCd1qo5tpUy7RSfow8GLq3j11uEXSlrcmIrZGxJqIWDOPBVPNhpnNkEDUanOytrLppET3buAjkj4ELATOpyjhXShpbirVtRobY2Z9qk8LdFMv0UXErRGxLCJWUHT7PhgRHwUeohjzAsUYmAc6zqWZ9d5s7Yxo478Bv57GviyiGAtjZlUQmVvJdGXAcER8E/hmen+A4skENls0P1CiX1usbVJlLK3l8MwIM8sSQK3mQGezVb0EVy/Z1V9dsquWAFyiM7Oq69f/uxzorHua/womexhsv/7VzGZ9+iNzoDOzTOUcOpLDgW42aVfCmq6SVXPb3XRfz6Zfn/7oHOjMLE9AuNfV+s7ZlqwmK5m1O+8SXIU40JlZ1fXp/1kOdLPJVEtWqaSmgYHTv258/PTvbdcmZ9XhQGdmleYBw1ZJ9ZLc3HnF60DxDIhIJbf6r3zbkp1VTr/+aMv3hDwzK6+a8rYMktZJ2i9pRNLmFucXSLo/nd8taUXT+TdLekXSb052LZfo7HVzija45ra4utqJk8WbqJ1+QnPS507/BT+jpNeOx9n1DXXpR5NWD9wCXEOx5MIeSUMR8WRDso3AkYhYJWkDxcJbNzacvxP465zruURnZnlyn0WXFwyvBEYi4kBaDvU+YH1TmvUUC2xBseDW2rQAF5Kuo1i0em/OxVyiMzS3+DXQ/PnFgVSii2PHi9exk6edr7/GiWK53jhepGNO0/ekEl09nUtq/U5n0xmxWNJww/7WiNjasL8UONiwPwpc1fQdp9KkBa+PAoskvUbxgN9rgEmrreBAZ2ZnI///qkMRsWaC8zkLabVL8zsUa0e/osyhTA50VTbZc+Hq51MJTvOL3lVqRfpaKskNXHA+ACcu/1cAvPKmosR27gvF+QWPHABg/MiR9D2p5Je+P1IbHjHe0T/HSqA2eZJMo8Dyhv1WC2nV04xKmgtcABymKPldL+n3gAuBmqRjEfFH7S7mQGdmebo7jm4PMChpJfAcxQJbv9iUZohiga1vUSy49WAUY5veW08g6beBVyYKcuBAN7vVS3q100t8MTZ22v7Y6hUAPH1j8evy3ncVHWN/9523ATB4/C0A6Fs/LD5/svh8fdzdmb20TSVNz6joG93qdU1tbrcAO4EB4O6I2CvpdmA4IoYoFtbalhbaOkwRDKfEgc7M8nWxPykidgA7mo7d1vD+GHDDJN/x2znXcqCzUyWueknsVIkszYh49V8uBOCqn3gKgC++5W8B2FArfn1GlwwCcN689OtUHz/n0UtWEg50ZpatW1XXmeZAV2WT9bam81E7fe4qc1JvaRr/9mPfOwbA7sdXAfAfVTttf9X30ni7U21zaWbFnFSia+519Xi6/hRkT+8qGwc6M8vXp/9HOdDNRm1KVPWSWP11PJXoBh4vxskNfvFSAJ5501uB10ty8574bpE+tfVpftGmVy/RqV5C7N4YLOsRV13NrPoc6Kx0JpsZUU82p6ndJc2QmHPeeQDUXnmlSPetxwA4L/XG1p9O8npJbv7p33OymDkRtcy/DrfdlV+f/ogc6Mwsi8JVVyujzBJSvcRV72Wtl/DqTzWpl+xOjY+rpRLcvKZfn3Q80nPrTj2PruY5rpXhXlczq7pZWaKTdCHwBeCdFLX3/wzsB+4HVgDPAL8QEUc6yqVNr1Tiqh1PJbX6XNd6L2z9eXULFqTj9XFxqSSYnlt36gnELsFVV58Guk7n6PwB8PWIeBvwLmAfsBnYFRGDwK60b2b9Ll5vp5tsK5spl+gknQ+8D/hPAOlxyCckrQfen5LdA3yT4mmgVnb1Elq9RDfWNPd1rD4Hdu5p6WupRHdGSS6z19f6SJ/+KDsp0V0KvAT8iaRHJH1B0rnAGyPiBYD0ekmrD0vaJGlY0vBJjneQDTObKarlbWXTSaCbC1wBfD4iLgd+xFlUUyNia0SsiYg181jQQTas66TTnw1XG4faOHH8OHH8OLVXXy22116j9tprp863/bxZj3US6EaB0YjYnfa3UwS+70taApBeX+wsi2ZWGt1bBWxGTTnQRcT3gIOS3poOrQWe5PXHH5NeH+gohzbzIiZuV6ufb5dusvPWn2ZjZ0TyX4EvSZpPscbiL1EEzy9L2gg8yyRPCDWzPlLCIJajo0AXEY8CrZY0W9vJ95pZSc3GQGdms4coZ49qDj/U37qv3uvq3tdq6XIbnaR1kvZLGpF0xogNSQsk3Z/O75a0Ih2/UtKjafuOpH8/2bUc6MwsX5d6XSUNAFuAa4HVwE2SVjcl2wgciYhVwJ3AHen4E8CaiLgMWAf8cVrgui0HOmvPJTJr1r3hJVcCIxFxIM2qug9Y35RmPcXsKiiGr62VpIh4NSLqiw8vzLmiA52ZZTuLquvi+syntG1q+qqlwMGG/dF0rGWaFNiOAosAJF0laS/wOPDLDYGvJXdGWHtTHQPnsXPVlf+jPRQRrUZk1LWqKjR/e9s0aaLCOyS9HbhH0l+nBa9bconOzPJEV+e6jgLLG/aXAc+3S5Pa4C4ADp+WpYh9FNNP3znRxRzozCxf99ro9gCDklamCQcbKGZVNWqcZXU98GBERPrMXABJbwHeSvHsy7ZcdTWzbN2a3hURY5JuAXYCA8DdEbFX0u3AcEQMAXcB2ySNUJTkNqSPvwfYLOkkUAP+S0Qcmuh6DnRmlq+Lza8RsQPY0XTstob3x2gxhTQitgHbzuZaDnRmlqekTybJ4UBnZllEOZ9MksOBzsyyOdCZWfU50JlZ5TnQmVmllfTpwTkc6MwsnwOdmVVdvz5404HOzLK56mpm1eYBw2Y2KzjQmVmVeWaEmc0KqvVnpHOgM7M8bqMzs9nAVVczqz4HOjOrOpfozKz6HOjMrNLCU8DMrOL6eRydlzs0s3wReVsGSesk7Zc0Imlzi/MLJN2fzu+WtCIdv0bSw5IeT68fmOxaHQU6Sb8maa+kJyTdK2lhWnNxt6SnUibnd3INMysPRd426fdIA8AW4FpgNXCTpNVNyTYCRyJiFXAncEc6fgj4+Yj4CYp1XyddEWzKgU7SUuBXgDUR8U6KtRk3pMzcGRGDwJGUWTPrd7mLV+cV6K4ERiLiQEScAO4D1jelWQ/ck95vB9ZKUkQ8EhHPp+N7gYWSFkx0sU6rrnOBH0urZp8DvAB8IGWKlMnrOryGmZWEankbsFjScMO2qemrlgIHG/ZH07GWaSJiDDgKLGpK8x+ARyLi+ET5nnJnREQ8J+n3gWeB14D/BTwM/CBlql3mAUj/8E0ACzlnqtkwsxl0Fr2uhyJizURf1eJYc1lwwjSS3kFRg/zgZJnppOp6EUXRciXwJuBcivp224yddjBia0SsiYg185iw1GlmZRB0szNiFFjesL8MeL5dmlRrvAA4nPaXAX8BfDwinp7sYp1UXX8G+G5EvBQRJ4GvAj8NXJgy1S7zZtanutUZAewBBlPn5XyK9v2hpjRDFJ0NANcDD0ZESLoQ+Cvg1oj4PzkX6yTQPQtcLekcSQLWAk8CD6VMkTL5QAfXMLMy6VJnRGreugXYCewDvhwReyXdLukjKdldwCJJI8CvA/UhKLcAq4DfkvRo2i6Z6HqdtNHtlrQd+DYwBjwCbKWItPdJ+kw6dtdUr2Fm5dHtAcMRsQPY0XTstob3x4AbWnzuM8BnzuZaHc2MiIhPA59uOnyAouvYzKokwg/eNLNZoD/jnAOdmeXr17muDnRmlicAV13NrPL6M8450JlZPlddzazy3OtqZtXm5Q7NrOqKAcP9Gekc6Mwsn9eMMLOqc4nOzKrNbXRmVn2e62pms4GrrmZWaV7A2sxmBZfozKzy+jPOOdCZWT7V+rPu6kBnZnkCDxg2s2oT0bcDhjtZBczMZpvureuKpHWS9ksakbS5xfkFku5P53dLWpGOL5L0kKRXJP1RzrUc6MwsX5cCnaQBYAvFovergZskrW5KthE4EhGrgDuBO9LxY8BvAb+Zm20HOjPLU2+jy9kmdyUwEhEHIuIEcB+wvinNeuCe9H47sFaSIuJHEfH3FAEvi9vozCzbWfS6LpY03LC/NSK2NuwvBQ427I8CVzV9x6k0ETEm6SiwCDh0VpnGgc7MsuW3vwGHImLNBOfV+gJnnSaLq65mlifoZmfEKLC8YX8Z8Hy7NJLmAhcAh6eSdQc6M8vXvTa6PcCgpJWS5gMbgKGmNEPAzen99cCDEVMb3+Kqq5ll69Y4utTmdguwExgA7o6IvZJuB4YjYgi4C9gmaYSiJLfhVD6kZ4DzgfmSrgM+GBFPtrueA52Z5evigOGI2AHsaDp2W8P7Y8ANbT674myu5UBnZnkiYLw/54A50JlZvj6dAuZAZ2b5+jTQTdrrKuluSS9KeqLh2MWSviHpqfR6UTouSX+Y5qY9JumK6cy8mc2gAGqRt5VMzvCSPwXWNR3bDOyKiEFgV9qHYt7aYNo2AZ/vTjbNrPcCopa3lcykgS4i/pYzB+k1zkG7B7iu4fgXo/B/gQslLelWZs2sh4KiMyJnK5mpDhh+Y0S8AJBeL0nHW81fW9rqCyRtkjQsafgkx6eYDTObUV18TNNM6nZnRPbctDTBdyvA+bq4fHfGzM5UwiCWY6oluu/Xq6Tp9cV0PGf+mpn1pczSXAmD4VQDXeMctJuBBxqOfzz1vl4NHK1Xcc2szwVQq+VtJTNp1VXSvcD7KZ4vNQp8Gvhd4MuSNgLP8vo0jR3Ah4AR4FXgl6Yhz2bWKyUsreWYNNBFxE1tTq1tkTaAT3SaKTMrI08BM7OqC4gSjpHL4UBnZvlKOOshhwOdmeWrahudmRlQBLkS9qjmcKAzs3wu0ZlZtQUxPt7rTEyJA52Z5ak/pqkPOdCZWT4PLzGzKgsgXKIzs0qLcInOzKqvXzsjNMWFr7ubCekl4EfAoV7nZQKLcf464fx1phv5e0tEvGGqH5b09ZSPHIcionkJhp4pRaADkDQcEWt6nY92nL/OOH+dKXv+ym6qz6MzM+sbDnRmVnllCnRbe52BSTh/nXH+OlP2/JVaadrozMymS5lKdGZm08KBzswqrxSBTtI6SfsljUja3OO8LJf0kKR9kvZK+mQ6frGkb0h6Kr1e1ON8Dkh6RNLX0v5KSbtT/u6XNL+HebtQ0nZJ/5ju40+V6f5J+rX0s31C0r2SFvby/km6W9KLkp5oONbyfqUV9v4w/a08JumKmcpnP+t5oJM0AGwBrgVWAzdJWt3DLI0BvxERbweuBj6R8rMZ2BURg8CutN9LnwT2NezfAdyZ8ncE2NiTXBX+APh6RLwNeBdFPktx/yQtBX4FWBMR7wQGgA309v79KdA8uLbd/boWGEzbJuDzM5TH/hYRPd2AnwJ2NuzfCtza63w15OcB4BpgP7AkHVsC7O9hnpZR/PJ/APgaIIpR83Nb3dMZztv5wHdJHV0Nx0tx/4ClwEHgYoopkF8DfrbX9w9YATwx2f0C/hi4qVU6b+23npfoeP0Xr240Hes5SSuAy4HdwBsjLcadXi/pXc74HPApoD7DehHwg4gYS/u9vIeXAi8Bf5Kq1l+QdC4luX8R8Rzw+xTrEb8AHAUepjz3r67d/Srt30uZlSHQqcWxno95kXQe8BXgVyPih73OT52kDwMvRsTDjYdbJO3VPZwLXAF8PiIup5jD3Otq/imprWs9sBJ4E3AuRXWwWc9/B9so08+6b5Qh0I0Cyxv2lwHP9ygvAEiaRxHkvhQRX02Hvy9pSTq/BHixR9l7N/ARSc8A91FUXz8HXCip/jSaXt7DUWA0Inan/e0Uga8s9+9ngO9GxEsRcRL4KvDTlOf+1bW7X6X7e+kHZQh0e4DB1Os1n6JheKhXmZEk4C5gX0R8tuHUEHBzen8zRdvdjIuIWyNiWUSsoLhXD0bER4GHgOtLkL/vAQclvTUdWgs8SUnuH0WV9WpJ56SfdT1/pbh/DdrdryHg46n39WrgaL2KaxPodSNhFA2qHwL+CXga+B89zst7KKoCjwGPpu1DFO1gu4Cn0uvFJbhv7we+lt5fCvwDMAL8ObCgh/m6DBhO9/AvgYvKdP+A3wH+EXgC2AYs6OX9A+6laC88SVFi29juflFUXbekv5XHKXqPe/p72A+bp4CZWeWVoepqZjatHOjMrPIc6Mys8hzozKzyHOjMrPIc6Mys8hzozKzy/j+GQFgUzFInrAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T19:55:22.066257Z",
     "start_time": "2020-02-12T19:55:22.005468Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input coordinates = (0.0, 0.0)\n",
      "Coarse centroid found at (54, 54)\n",
      "Fine centroid found at (x, y) = (54.0040, 53.9883). Rms error = 76.3621\n"
     ]
    }
   ],
   "source": [
    "cent = centroid_b(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T20:04:18.641266Z",
     "start_time": "2020-02-12T20:04:18.636975Z"
    }
   },
   "outputs": [],
   "source": [
    "def compute_centroid_diff(data, head):\n",
    "    \n",
    "    offset = head['R_OFFSET'] #y offset of psf in pixels\n",
    "    \n",
    "    origin = (54,54) #center pixel that psfs offset from\n",
    "    \n",
    "    true_cent = np.add(origin, [0, offset])\n",
    "    \n",
    "    cent = centroid_b(data,silent=True)\n",
    "    \n",
    "    diff = np.subtract(cent,true_cent)*0.11 #diff in arcsec\n",
    "    \n",
    "    rms = np.sqrt(np.mean(diff**2))\n",
    "    \n",
    "    return rms*1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T20:04:21.111616Z",
     "start_time": "2020-02-12T20:04:19.732872Z"
    }
   },
   "outputs": [],
   "source": [
    "f560 = []\n",
    "f1000 = []\n",
    "f1500 = []\n",
    "f2550 = []\n",
    "\n",
    "for f in os.listdir():\n",
    "    if f.endswith('.fits'):\n",
    "        head, data = open_im(f)\n",
    "        \n",
    "        filt = head['FILTER']\n",
    "        \n",
    "        if filt == 'F560W':\n",
    "            mas_acc = compute_centroid_diff(data, head)\n",
    "            f560.append(mas_acc)\n",
    "            \n",
    "        elif filt == 'F1000W':\n",
    "            mas_acc = compute_centroid_diff(data, head)\n",
    "            f1000.append(mas_acc)\n",
    "            \n",
    "        elif filt == 'F1500W':\n",
    "            mas_acc = compute_centroid_diff(data, head)\n",
    "            f1500.append(mas_acc)\n",
    "            \n",
    "        elif filt == 'F2550W':\n",
    "            mas_acc = compute_centroid_diff(data, head)\n",
    "            f2550.append(mas_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T20:15:23.067077Z",
     "start_time": "2020-02-12T20:15:23.063393Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.081403171713944"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(f560)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T20:15:23.272722Z",
     "start_time": "2020-02-12T20:15:23.268484Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.66258781182637"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(f1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T20:15:23.496343Z",
     "start_time": "2020-02-12T20:15:23.492291Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.662985304370933"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(f1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-12T20:15:23.728559Z",
     "start_time": "2020-02-12T20:15:23.724682Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.751846507303288"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(f2550)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}