Fresh GFA Target Functions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def find_points_in_tel_range(telra, teldec, ra, dec, radius = None):\n",
    "    \"\"\"Return a list of indices of points that are within a radius of an arbitrary telra, teldec.\n",
    "\n",
    "    This function is optimized to query a lot of points with a single telra and teldec.\n",
    "\n",
    "    radius is in units of degrees. The return value is a list\n",
    "    that contains the index of points that are in each tile.\n",
    "    The indices are not sorted in any particular order.\n",
    "\n",
    "    if tiles is a scalar, a single list is returned.\n",
    "\n",
    "    default radius is from desimodel.focalplane.get_tile_radius_deg()\n",
    "    \n",
    "    Note: This is simply a modified version of find_points_in_tiles, but this function does not know about tiles. \n",
    "    \"\"\"\n",
    "    from scipy.spatial import cKDTree as KDTree\n",
    "    import desimodel.focalplane\n",
    "    import numpy as np\n",
    "\n",
    "    if radius is None:\n",
    "        radius = desimodel.focalplane.get_tile_radius_deg()\n",
    "\n",
    "    # check for malformed input shapes. Sorry we currently only\n",
    "    # deal with vector inputs. (for a sensible definition of indices)\n",
    "\n",
    "    assert ra.ndim == 1\n",
    "    assert dec.ndim == 1\n",
    "\n",
    "    points = _embed_sphere(ra, dec)\n",
    "    tree = KDTree(points)\n",
    "\n",
    "    # radius to 3d distance\n",
    "    threshold = 2.0 * np.sin(np.radians(radius) * 0.5)\n",
    "    xyz = _embed_sphere(telra, teldec)\n",
    "    indices = tree.query_ball_point(xyz, threshold)\n",
    "    return indices\n",
    "\n",
    "def _embed_sphere(ra, dec):\n",
    "    \"\"\" embed RA DEC to a uniform sphere in three-d \"\"\"\n",
    "    import numpy as np\n",
    "    \n",
    "    phi = np.radians(np.asarray(ra))\n",
    "    theta = np.radians(90.0 - np.asarray(dec))\n",
    "    r = np.sin(theta)\n",
    "    x = r * np.cos(phi)\n",
    "    y = r * np.sin(phi)\n",
    "    z = np.cos(theta)\n",
    "    return np.array((x, y, z)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def on_gfa(telra, teldec, ra, dec, buffer_arcsec = 100):\n",
    "    \"\"\"\n",
    "    Checks if a target is on any of the 10 GFAs given telra, teldec and an array of RA and Dec pointings,\n",
    "    as well as a parameter for degrees of tolerance one would like to allow. When using \n",
    "    desimodel.footprint.find_points_in_tiles(tiles, ra, dec, radius) with this function to \n",
    "    check what points are on the GFAs, the default radius parameter should be set to 1.651 (degrees),\n",
    "    so that boundary GFA area actually encompasses points normally outside of the tile. \n",
    "    Parameters:\n",
    "    telra: The telescope's arbitrary RA pointing\n",
    "    teldec: The telescope's arbitrary Dec pointing\n",
    "    ra: An array of RA values for locations in the sky\n",
    "    dec: An array of declination values for locations in the sky\n",
    "    buffer_arcsec: A value in arcseconds on the sky of how much tolerance \n",
    "    one would allow for seeing if a target is on the gfa.  \n",
    "    Returns: \n",
    "    targetindices: a list  of targets with their respective indices in the \n",
    "    RA and Dec list passed in that fall on certain GFAs denoted by the index \n",
    "    in the gfaindices list.\n",
    "    gfaindices: a list equal in length with the targetindices list with the gfa location 0-9 as each element\n",
    "    \"\"\"\n",
    "    import desimodel.footprint\n",
    "    # If any calculated area is under the threshold area, it is mathematically impossible\n",
    "    THRESHOLD_AREA = 469.7\n",
    "    MIN_TOLERANCE = 0.001\n",
    "    \n",
    "    #inrangeindices = desimodel.footprint.find_points_in_tel_range(telra, teldec, ra, dec)\n",
    "    inrangeindices = find_points_in_tel_range(telra, teldec, ra, dec, 1.651)\n",
    "    if not inrangeindices:\n",
    "        return [], []\n",
    "    inrangeindices = np.asarray(inrangeindices)\n",
    "    \n",
    "    targetx, targety = desimodel.focalplane.radec2xy(telra, teldec, ra[inrangeindices], dec[inrangeindices])\n",
    "    \n",
    "    x_tolerance, y_tolerance = degrees2xytolerance(buffer_arcsec)\n",
    "    \n",
    "    targetindices = []\n",
    "    gfaindices = []\n",
    "    \n",
    "    # x and y hold the 40 new GFA coordinates\n",
    "    x, y = shift_gfa_points(x_tolerance, y_tolerance)\n",
    "    # The area boundary's value is the area of the gfa plus some tolerance. \n",
    "    AREA_BOUNDARY = retrieve_minimum_boundary(x_tolerance, y_tolerance) + MIN_TOLERANCE\n",
    "    \n",
    "    targetx = np.asarray(targetx)\n",
    "    targety = np.asarray(targety)\n",
    "    # Method to check if point is inside the rectangle\n",
    "    for gfaid in range(0, 40, 4):\n",
    "        # a1 through a4 are edge lengths of the rectangle formed by corners of the GFAs\n",
    "        a1 = np.sqrt((x[gfaid] - x[gfaid + 1])**2\n",
    "                    + (y[gfaid] - y[gfaid + 1])**2)\n",
    "        a2 = np.sqrt((x[gfaid + 1] - x[gfaid + 2])**2\n",
    "                    + (y[gfaid + 1] - y[gfaid + 2])**2)\n",
    "        a3 = np.sqrt((x[gfaid + 2] - x[gfaid + 3])**2\n",
    "                    + (y[gfaid + 2] - y[gfaid + 3])**2)\n",
    "        a4 = np.sqrt((x[gfaid + 3] - x[gfaid])**2\n",
    "                    + (y[gfaid + 3] - y[gfaid])**2)\n",
    "        # b1 through b4 are the line segments from each corner to the target location\n",
    "        b1 = np.sqrt((x[gfaid] - targetx)**2\n",
    "                    + (y[gfaid] - targety)**2)\n",
    "        b2 = np.sqrt((x[gfaid + 1] - targetx)**2\n",
    "                    + (y[gfaid + 1] - targety)**2)\n",
    "        b3 = np.sqrt((x[gfaid + 2] - targetx)**2\n",
    "                    + (y[gfaid + 2] - targety)**2)\n",
    "        b4 = np.sqrt((x[gfaid + 3] - targetx)**2\n",
    "                    + (y[gfaid + 3] - targety)**2)\n",
    "        # Calculating areas of triangles using Heron's Formula\n",
    "        u1 = (a1 + b1 + b2) / 2.0\n",
    "        u2 = (a2 + b2 + b3) / 2.0\n",
    "        u3 = (a3 + b3 + b4) / 2.0\n",
    "        u4 = (a4 + b4 + b1) / 2.0\n",
    "        area1 = np.sqrt((u1 * (u1 - a1) * (u1 - b1) * (u1 - b2)).clip(0))\n",
    "        area2 = np.sqrt((u2 * (u2 - a2) * (u2 - b2) * (u2 - b3)).clip(0))\n",
    "        area3 = np.sqrt((u3 * (u3 - a3) * (u3 - b3) * (u3 - b4)).clip(0))\n",
    "        area4 = np.sqrt((u4 * (u4 - a4) * (u4 - b4) * (u4 - b1)).clip(0))\n",
    "        targetarea = area1 + area2 + area3 + area4\n",
    "        \n",
    "        \n",
    "        assert np.all(targetarea > THRESHOLD_AREA)\n",
    "        \n",
    "        if(any(targetarea < AREA_BOUNDARY) and all(targetarea > THRESHOLD_AREA)):\n",
    "            newtargetindices = np.where(targetarea < AREA_BOUNDARY)\n",
    "            targetindices.extend(newtargetindices[0])\n",
    "            gfaindices.extend([int(gfaid / 4)] * len(newtargetindices[0]))\n",
    "    return inrangeindices[targetindices], gfaindices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def retrieve_minimum_boundary(x_tolerance, y_tolerance):\n",
    "    \"\"\"\n",
    "    Used as a helper function to the on_gfa function to find the minimum boundary\n",
    "    area for a point to lie inside a certain GFA given an tolerance in x and y in mm\n",
    "    Parameters:\n",
    "    x_tolerance: tolerance in x in mm\n",
    "    y_tolerance: tolerance in y in mm\n",
    "    Returns: \n",
    "    targetarea: the minimum boundary area for the procedure to check if a point is inside the GFA\n",
    "    \"\"\"\n",
    "    import desimodel.footprint\n",
    "    import desimodel.focalplane\n",
    "    \n",
    "    targetx = 116.279135121\n",
    "    targety = -372.885546514\n",
    "    #6.644525362152656, -9.055425745149217 GUARANTEED TO BE IN GFA (RA, DEC)\n",
    "    #x, y = desimodel.focalplane.radec2xy(7.11, -10.53, targetx, targety)\n",
    "    # If any calculated area is under the threshold area, it is mathematically impossible\n",
    "    THRESHOLD_AREA = 469.7\n",
    "    MIN_TOLERANCE = 0.001\n",
    "    # The area boundary's value is the area of the gfa plus some tolerance. \n",
    "    \n",
    "    # x and y hold the 40 new GFA coordinates\n",
    "    x, y = shift_gfa_points(x_tolerance, y_tolerance)\n",
    "    \n",
    "    targetx = np.asarray(targetx)\n",
    "    targety = np.asarray(targety)\n",
    "    # Method to check if point is inside the rectangle\n",
    "    for gfaid in range(0, 4, 4):\n",
    "        # a1 through a4 are edge lengths of the rectangle formed by corners of the GFAs\n",
    "        a1 = np.sqrt((x[gfaid] - x[gfaid + 1])**2\n",
    "                    + (y[gfaid] - y[gfaid + 1])**2)\n",
    "        a2 = np.sqrt((x[gfaid + 1] - x[gfaid + 2])**2\n",
    "                    + (y[gfaid + 1] - y[gfaid + 2])**2)\n",
    "        a3 = np.sqrt((x[gfaid + 2] - x[gfaid + 3])**2\n",
    "                    + (y[gfaid + 2] - y[gfaid + 3])**2)\n",
    "        a4 = np.sqrt((x[gfaid + 3] - x[gfaid])**2\n",
    "                    + (y[gfaid + 3] - y[gfaid])**2)\n",
    "        # b1 through b4 are the line segments from each corner to the target location\n",
    "        b1 = np.sqrt((x[gfaid] - targetx)**2\n",
    "                    + (y[gfaid] - targety)**2)\n",
    "        b2 = np.sqrt((x[gfaid + 1] - targetx)**2\n",
    "                    + (y[gfaid + 1] - targety)**2)\n",
    "        b3 = np.sqrt((x[gfaid + 2] - targetx)**2\n",
    "                    + (y[gfaid + 2] - targety)**2)\n",
    "        b4 = np.sqrt((x[gfaid + 3] - targetx)**2\n",
    "                    + (y[gfaid + 3] - targety)**2)\n",
    "        # Calculating areas of triangles using Heron's Formula\n",
    "        u1 = (a1 + b1 + b2) / 2.0\n",
    "        u2 = (a2 + b2 + b3) / 2.0\n",
    "        u3 = (a3 + b3 + b4) / 2.0\n",
    "        u4 = (a4 + b4 + b1) / 2.0\n",
    "        area1 = np.sqrt(u1 * (u1 - a1) * (u1 - b1) * (u1 - b2))\n",
    "        area2 = np.sqrt(u2 * (u2 - a2) * (u2 - b2) * (u2 - b3))\n",
    "        area3 = np.sqrt(u3 * (u3 - a3) * (u3 - b3) * (u3 - b4))\n",
    "        area4 = np.sqrt(u4 * (u4 - a4) * (u4 - b4) * (u4 - b1))\n",
    "        targetarea = area1 + area2 + area3 + area4\n",
    "        \n",
    "        assert np.all(targetarea > THRESHOLD_AREA)\n",
    "        return targetarea"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def degrees2xytolerance(buffer_arcsec):\n",
    "    \"\"\"\n",
    "    Used as a helper function to the on_gfa function to find the tolerance in x and y\n",
    "    given a tolerance in arcseconds\n",
    "    Parameters:\n",
    "    buffer_arcsec: a tolerance in arcseconds for checking if a point is on the GFA\n",
    "    Returns: \n",
    "    x_tolerance: tolerance in x in mm\n",
    "    y_tolerance: tolerance in y in mm\n",
    "    \"\"\"\n",
    "    # Uses the center of a given GFA from DESI-0530-v13 Excel Spreadsheet to find the tolerance \n",
    "    import desimodel.io\n",
    "    import scipy.interpolate\n",
    "    platescale = desimodel.io.load_platescale()\n",
    "    fn = scipy.interpolate.interp1d(platescale['radius'], platescale['radial_platescale'], kind = 'quadratic')\n",
    "    fn1 = scipy.interpolate.interp1d(platescale['radius'], platescale['az_platescale'], kind = 'quadratic')\n",
    "    # Center of a given GFA from DESI-0530-v13 Excel Spreadsheet\n",
    "    x = 333.738\n",
    "    y = 217.766\n",
    "    radius = np.sqrt(x**2 + y**2)\n",
    "    # Platescales are in units of microns per arcsecond\n",
    "    r_ps = fn(radius)\n",
    "    az_ps = fn(radius)\n",
    "    x_tolerance = buffer_arcsec / (10**3) * r_ps\n",
    "    y_tolerance = buffer_arcsec / (10**3) * az_ps\n",
    "    return x_tolerance, y_tolerance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def shift_gfa_points(deltax, deltay):\n",
    "    \"\"\"\n",
    "    Used as a helper function to the on_gfa function to find the new \n",
    "    GFA locations after incorporating a tolerance in x and y\n",
    "    Parameters:\n",
    "    deltax: tolerance in x in mm\n",
    "    deltay: tolerance in y in mm\n",
    "    Returns: \n",
    "    Returns the 40 new GFA locations in x and y\n",
    "    \"\"\"\n",
    "    import numpy as np\n",
    "    x = [-125.10482863, -129.83038525, -159.04283509, -154.31646944]\n",
    "    y = [-370.01790486, -384.56223777, -375.05643893, -360.51151824]\n",
    "    point1 = [x[2], y[2]]\n",
    "    point2 = [x[1], y[1]]\n",
    "    vector1 = [(point2[0] - point1[0]), (point2[1] - point1[1])]\n",
    "    vector2 = [1, 0]\n",
    "    # Angle between vector1 and vector 2 using dot product\n",
    "    angle = np.arccos((np.dot(vector1, vector2))/(np.sqrt((vector1[0]**2) + (vector1[1]**2))))\n",
    "    \n",
    "    shiftmat = np.zeros(shape=(2,2))\n",
    "    shiftmat[0] = [np.cos(angle), -np.sin(angle)]\n",
    "    shiftmat[1] = [np.sin(angle), np.cos(angle)]\n",
    "    reverseshift= np.zeros(shape=(2,2))\n",
    "    reverseshift[0] = [np.cos(angle), np.sin(angle)]\n",
    "    reverseshift[1] = [-np.sin(angle), np.cos(angle)]\n",
    "    \n",
    "    # Shifts the initial coordinates to be parallel to the vector [1, 0]\n",
    "    coord = np.zeros(shape=(2,1))\n",
    "    oldxcoord = x\n",
    "    oldycoord = y\n",
    "    for i in range(4):\n",
    "        coord[0] = oldxcoord[i]\n",
    "        coord[1] = oldycoord[i]\n",
    "        newcoord = np.matmul(shiftmat, coord)\n",
    "        oldxcoord[i] = newcoord[0]\n",
    "        oldycoord[i] = newcoord[1]\n",
    "        if(i == 0 or i == 1):\n",
    "            x[i] = newcoord[0] + deltax\n",
    "        else:\n",
    "            x[i] = newcoord[0] - deltax\n",
    "        if(i == 1 or i == 2):\n",
    "            y[i] = newcoord[1] - deltay\n",
    "        else:\n",
    "            y[i] = newcoord[1] + deltay\n",
    "    oldxcoord = x\n",
    "    oldycoord = y\n",
    "    for i in range(4):\n",
    "        coord[0] = oldxcoord[i]\n",
    "        coord[1] = oldycoord[i]\n",
    "        newcoord = np.matmul(reverseshift, coord)\n",
    "        oldxcoord[i] = newcoord[0]\n",
    "        oldycoord[i] = newcoord[1]\n",
    "        x[i] = newcoord[0]\n",
    "        y[i] = newcoord[1]\n",
    "    \n",
    "    rotatemat = np.zeros(shape=(2,2))\n",
    "    rotatemat[0] = [np.cos(np.radians(36)), -np.sin(np.radians(36))]\n",
    "    rotatemat[1] = [np.sin(np.radians(36)), np.cos(np.radians(36))]\n",
    "    return find_new_gfa_coordinates(x, y, rotatemat)\n",
    "\n",
    "def find_new_gfa_coordinates(x, y, rotatemat):\n",
    "    \"\"\"\n",
    "    Used as a helper function to the on_gfa function to find the new\n",
    "    GFA coordinates given a list of x coordinates, y coordinates, and a rotation matrix\n",
    "    Parameters:\n",
    "    x: a list of x coordinates for the GFAs\n",
    "    y: a list of y coordinates for the GFAs\n",
    "    rotatemat: a matrix for rotating the respective coordinates\n",
    "    Returns: \n",
    "    x_all: a complete list of the 40 GFA x coordinates\n",
    "    y_all: a complete list of the 40 GFA y coordinates\n",
    "    \"\"\"\n",
    "    import numpy as np\n",
    "    x_all = np.zeros(shape=(40,1))\n",
    "    y_all = np.zeros(shape=(40,1))\n",
    "    coord = np.zeros(shape=(2,1))\n",
    "    gfacoord = np.zeros(shape=(4, 2))\n",
    "    oldxcoord = x\n",
    "    oldycoord = y\n",
    "    counter = 0\n",
    "    for j in range(10):\n",
    "        for i in range(4):\n",
    "            coord[0] = oldxcoord[i]\n",
    "            coord[1] = oldycoord[i]\n",
    "            newcoord = np.matmul(rotatemat, coord)\n",
    "            oldxcoord[i] = newcoord[0]\n",
    "            oldycoord[i] = newcoord[1]\n",
    "            gfacoord[i] = [newcoord[0], newcoord[1]]\n",
    "            x_all[counter] = newcoord[0]\n",
    "            y_all[counter] = newcoord[1]\n",
    "            counter += 1\n",
    "    return x_all, y_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHzhJREFUeJzt3X9wHGed5/H3V5LlkNyyASXEP4Rx2IWtheWMfMJHc/wY\nl52QJRxO4LIV6m4Vgx0n5kdVtoC9cGwoVXEUuWxgDdwRW7aJrapluRxsIJBsEsd4krDXi3GseHES\nOAJoFzl2ElSw1Ca2ZGu+98fTI02PZ6yxR6NpaT6vKlVPdz8987ht90fPj+4xd0dERKSordkVEBGR\nbFEwiIhIioJBRERSFAwiIpKiYBARkRQFg4iIpCgYREQkRcEgIiIpCgYREUnpaHYFzsVFF13ky5cv\nb3Y1RETmlMcee+xX7n7xdOXmZDAsX76cAwcONLsaIiJzipn9Uy3l1JUkIiIpCgYREUlRMIiISIqC\nQUREUhQMIiKSomAQEZGUOTldVSQljiGfh64uGB2FXA6iqNm1EpmzFAwyt5SGwNAQHDsG994LJ0+G\n/WbQ2Qn79ikcRM5Rw4LBzFYAW4F/AwwD/9ndf1uh3BXAF4F2YIe739qoOskcF8ewZg2MjUGhULmM\ne9g/OKhgEDlHjWwx7AA+7u4Pm9kHgU8At5QWMLN24H8BlwEjwA/N7B53f7KB9ZK5oLx7qKsLvvlN\nOHEC3Im7YXBFKNpzFEbPh9wwRCPJ8QcPhvdQOIictUYGw2uBR5LXe4AHKAsGYBXwtLv/HMDMvg6s\nAxQMrSqOw2/7O3fCyZPE3ZBfDl0vhot/Vw8MLYbtK2GidOqEQ2cB8ruScNi/H1avVpdSrSoFscZr\nWlYjg+EJwkX+W8A1wCsrlFkK/LJkfQT49w2sk2RZsauopFWwej2MtQFWUs6TdStZb4NxC62IyVaD\nupRqMzAAH/oQTExMBnFuGKIjBh0dsGED9PQoMFpIXcFgZg8Biyrs+hTwQeBLZnYLcA8wXudnbQI2\nASxbtqyet5KsyudhfDyMExAu8mPthADwkmVbsqzg2AXlG441qLLzRBzD5s3ESwoMroCdK+GUQUcB\nNgw5PUdPMnp4K137y7rr2tvhYx+DCy9UWMxDdQWDu6+dpsjlAGb2WuDKCvuPkG5JdCfbKn3WADAA\n0NvbW+WyIHNaLhdmFB0/DlS4yNdg0QszW6V5b3CQeEmBNdfBiXbwpCV20mBrbyhiHra3OSycgL27\nIRqZgNtuC7PA3KGtDRYuhL17FQ7zQMNucDOzVyTLNuAvCDOUyv0QeI2ZXWpmncC1hNaFtKIoCheW\nVauAGi7yzmTLwQrQOQF9h8rKLKrUoJVS+eWhZeblV4Okq664vdAWwuOmK2DzlRB3M9m6o1AIrb18\nfnYqLQ3VyDuf329m/w/4MfAMcCeAmS0xs/sA3P0U8BHCwPRTwF3u/kQD6yRZF0WhT5twke+cCBd9\nAJKlFaBjAq56Cm48ANu+A5/9XsnAM0zdz9DXN9t/grkjjuHYMbpehEKxm44zL91g/9LQmli9PgkH\nCC2Gzs7Q6pM5r2GDz+7+RcL9CeXbnwHeVbJ+H3Bfo+ohc9DoKLS1EY0UyO8qm5X0YoWpqQCXXw63\nvE8DpLUqGegf/Q+hm6hQHLspH9OptATG28PfTTQC9PbCli063/OE7nyW7MnlQn/18eNEI2UBAFP9\n2kULF0J/vy5KZ6NkoD83HMYOxkrDYToexh5yw8n6kiU6//OIgkGypzjWULyf4dSpMAtm48b0tMmh\noVC+r08XpbNVHOgfGyMaKbB391TLbGgx7OgJ94m4EbrwSpeJ9/y4JLQ1ljOvmPvcm+DT29vr+s7n\nFlG88UrdQjPvDM+dirshfyl0vTDVhTe0GHb2wKk2WFB6M2FnZ3gf/f1knpk95u6905ZTMIjIpEp3\nQA8NnXYnem4YoqPtcP31arHNIQoGEZk55a0LUCDMQbUGg8YYRGR6UaQQaCH6BjcREUlRMIiISIqC\nQUREUhQMIiKSomAQEZEUBYOIiKQoGEREJEXBICIiKQoGERFJUTCIiEiKgkFERFIUDCIikqJgEBGR\nFAWDiIikKBiyIo7hc58LSxGRJtL3MWRBHMOaNeHL2Ts7w/cd69n3ItIkDWsxmNkKM4vN7Edm9h0z\ne2mVcsNJmcfNrPW+li2Oob8fTpyAiYmw7O9Xy0FEmqaRLYYdwMfd/WEz+yDwCeCWKmVXu/uvGliX\nbIpjWL0axsZKvkvXiR58EL73PXj3u2HRIn2FoojMqkYGw2uBR5LXe4AHqB4MrSeO4aabJkNh9XoY\nb4OOAmwYgr5Dp4i+9a1Qdvt2fem6iMyaRg4+PwGsS15fA7yySjkHHjKzx8xsUwPrkw1xDJs3wzve\nAfv3E3fDTVfAWDt4G5xsh629ISji7uSYiQnYti2MQ6iLSUQarK4Wg5k9BCyqsOtTwAeBL5nZLcA9\nwHiVt3mrux8xs1cAe8zsx+7+SHmhJDQ2ASxbtqyeajdPcZD5xAlwJ+6GNdfB8faychaCYnAFRCPJ\nNvcwOJ3Pq9UgIg1VVzC4+9ppilwOYGavBa6s8h5HkuVzZnY3sIqpLqjScgPAAEBvb6/XUe3mGRxM\nhUJ/Dk60E9pt0/2JzMKMpVyu4dUUkdbWyFlJr0iWbcBfAFsrlLnAzH6n+JoQJIcbVaemimPYvp14\nqbP5SnjHB+DBV4MbU6FQtuw5mqy3tcENN2gaq4jMikYOPr/fzD6cvP5b4E4AM1sC7HD3dwGXAHeb\nWbEuX3P3+xtYp+ZIBprjxROsXg9jbYAlP1552TYBo+cnx7/nPXDHHc2ouYiUiuPQndvVBaOjoQU/\nD39Za1gwuPsXgS9W2P4M8K7k9c+BFY2qQyYUxxWOH2fwyjB2QGkroQJzWFiA3HCyYVGlYRwRmTVx\nHLqCd+6EkyfDNjPo6IANG+bdjEHd+dxo+XwYNK5FIUxX3TgEfYeSgefOzvCPTkSao2zSyCT3EBLb\ntsHu3fOqq1fB0Gi5XLi4Hz9O3yH4ak+4X6G01WCF0vsXkkDo6IAbN86730RE5oxiK+HgwdSkkfxy\n6HoxdPXmhiEa8aknFvT3z4v/r+Y+9yb49Pb2+oEDc+jpGcV/YNu3Ey+eSP3DSv0DO2LQ3g4bFQgi\nTRXH4V6jpNso7g7Tx3euhJPJ+KABHRNlv9AtXAj79mX2/66ZPebuvdOWUzDMotLfQH74w3SzFGDV\nKtiyJbP/qERaxubNsDVMpCzeb3SiPZlFWD5ZBOicgPyuJBwy/P+41mDQY7dnUxSF2UVbtsB5552+\nf+XKTP5jEmll+eUwnjyZoCIL+2+6Inlawf798La3wcDALNZyZikYmiGKwkDVjTfCggVTN69pkFkk\nG/r6wv1DhG7ezgmgcOZD9i+Bt6+HzVdCvHgitDo2b56Tj7FRV1KzFedFz9P50CJz1sAAfOhDMDEx\nNcbQE55nBlS9D8mA807B3t3JuOF552VmxpK6kuaKKIJPfjIT/2hEpMSmTfDoo3DjjUTHOrjjXnh4\nF6w6cubD3MLMw/xy0s84m0MUDCIi1RTHBR95BK66imgEttxf1rVU9iibtgJ0Fm9QbWubk884UzCI\niEwnisJso7Y2opEwA+nyn4enFGBhefnPYNt34L9/L+lGGgHWrs1MN9LZ0A1uIiK1yOXCfQrHjxON\nQH8eHn0VjHtoIfTnSx6TD+GepDl6w5taDCIitSibTRiNhJbBZ/aVtBCK2tvhK1+Zk6EAmpUkInL2\nijer7tgRvmGx+MSCnp5MP3W11llJ6koSETlbURR++vrm5XRzBYOIyLkqBsQ8ozEGERFJUTCIiEiK\ngkFERFIUDCIikqJgEBGRFAWDiIikKBhERCSlrmAws2vM7AkzK5hZb9m+T5rZ02b2EzN7Z5XjX25m\ne8zsp8nyZfXUR0RE6ldvi+Ew8F7gkdKNZvY64Frg9cAVwFfMrP30w7kZ2OvurwH2JusiItJEdQWD\nuz/l7j+psGsd8HV3H3P3XwBPA6uqlNudvN4NXFVPfUREpH6NGmNYCvyyZH0k2VbuEnc/mrw+BlzS\noPqIiEiNpn1Wkpk9BCyqsOtT7v7tmaqIu7uZVX3Uq5ltAjYBLFu2bKY+VkREykwbDO6+9hze9wjw\nypL17mRbuWfNbLG7HzWzxcBzZ6jHADAA4bHb51AnERGpQaO6ku4BrjWzhWZ2KfAaYH+Vctclr68D\nZqwFIiIi56be6apXm9kIEAH3mtkDAO7+BHAX8CRwP/Bhd59IjtlRMrX1VuAyM/spsDZZFxGRJtI3\nuImItIhav8FNdz6LiEiKgkFERFIUDCIikqJgEBGRFAWDiIikKBhERCRFwSAiIikKBhERSVEwiIhI\nioJBRERSFAwiIpKiYBARkRQFg4iIpCgYREQkRcEgIiIpCgYREUmZ9juf55U4hnweurpgdLTycmgo\nlO3rgyhqanVFRJqhdYIhjiGXg/Hx2srv3AkPP6xwEJGW0zrBMDh4WijE3ZBfDl0vwuj5U8vcMEQj\nJ8MxCgYRaTGtEwzHjqVW425YvR7G2gCb2m5A5wTs2wVR2TEiIq2gZQefB1fAWDtTZyAJB7ewfXBF\ns2omItJcLRsMxy6ob7+IyHzVOsGwaFF69YVpir9w+jEiIq2grmAws2vM7AkzK5hZb9m+T5rZ02b2\nEzN7Z5Xj+83siJk9nvy8q576nFFPT3r1aPLC00srhDGGvkOnHyMi0grqHXw+DLwX2Fa60cxeB1wL\nvB5YAjxkZq9194kK7/FX7n57nfWYXvH+hOLq4uSFEUIhWb7pGdhyP0Qjpx8jItIK6goGd38KwMzK\nd60Dvu7uY8AvzOxpYBUQ1/N5dSmbYVRtDGHJb5NQqHCMiEhD1HLzbekyl2voVPpGTVddCvxDyfpI\nsq2Sj5pZH3AA+Ji7/7pSITPbBGwCWLZsWd0VrDbGMN3Yg4jIjIpjWLMGxsagUAAzcK++bGuDhQth\n796GhcO0Ywxm9pCZHa7ws24GPv8O4NXAG4GjwOerFXT3AXfvdffeiy+++Ow/qWwgue8QLJwAqzS2\nUOUYEZEZl89PhQKEi/+ZloVCKJ/PN6xK07YY3H3tObzvEeCVJevdybby9362+NrMtgPfPYfPqk1f\nH9x5ZzihhO6ifbsgfyl0vVB6x3NSfsGCcIyISCN1dU2FQqL8qQypaxOE8l1dDatSo7qS7gG+ZmZf\nIAw+vwbYX17IzBa7e3F+0NWEwezGiCLYty/Vjxd1dRHpIXoi0kxlk1zKn8qQehrDSPXjZlJdwWBm\nVwNfBi4G7jWzx939ne7+hJndBTwJnAI+XJyRZGY7gK3ufgC4zczeSJgXNAzcUE99phVFutiLSLaU\nTXKZfCpDMlOy9GkMqWBo4OSYemcl3Q3cXWXfZ4HPVti+seT1n9bz+SIi8021GZOz+TSG1rnzWUQk\ni2p8KsNp2xs4OUbBICLSTDU+lWFye5XjZpKCQUSkmc70VIaS5eT2KsfNJAWDiEgz1fhUhtO2N3Dw\nWcEgIpIhWXgqg4JBRKSZpnkqA4WwnnoqQ4XjZlLrfLWniEgW9fXBV786+Z30k09lWF7lzmezhj+Z\nQcEgItJMURSeyDDdUxnmwdNVRUSkVhl7KoPGGEREJEXBICIiKQoGERFJUTCIiEiKgkFERFIUDCIi\nkqJgEBGRFAWDiIikKBhERCRFwSAiIikKBhERSVEwiIhIioJBRERS6goGM7vGzJ4ws4KZ9ZZs7zKz\nfWb2r2b2P89w/MvNbI+Z/TRZvqye+oiISP3qbTEcBt4LPFK2/QRwC/DxaY6/Gdjr7q8B9ibrIiLS\nRHUFg7s/5e4/qbD9BXf/PiEgzmQdsDt5vRu4qp76iIhI/Zo9xnCJux9NXh8DLqlW0Mw2mdkBMzvw\n/PPPz07tRERa0LTBYGYPmdnhCj/rZrIi7u6An2H/gLv3unvvxRdfPJMfLSIiJab9ak93X9vAz3/W\nzBa7+1EzWww818DPEhGRGjS7K+ke4Lrk9XXAt5tYFxERof7pqleb2QgQAfea2QMl+4aBLwDrzWzE\nzF6XbN9RMrX1VuAyM/spsDZZFxGRJpq2K+lM3P1u4O4q+5ZX2b6x5PUosKaeOoiIyMxqdleSiIhk\njIJBRERSFAwiIpKiYBARkRQFg4iIpCgYREQkRcEgIiIpCgYREUlRMIiISIqCQUREUhQMIiKSUtez\nkqSKOIZ8Hrq6YHQUcjmIombXSkSkJgqGmRTHMDgIO3fCyZNhmxl0dMCGDdDXp4AQkcxTMNSjtGUw\nNDQZCHE35JdD14swer6TGz5JtHVr2K+AEJGMUzCcqziGNWtgbAwKhanN3bDmOhhrh4JBm8PCCdi7\nG6KRk7BtG+zeDXv3KhxEJJM0+Hyu8vmKodCfgxPtUEjObKEtrPfnwn7cw3H5/OzXWUSkBgqGc9XV\nVbGlsOfV4AZ4ssPD+oO/B6vXJ+FQKMBvftOESouITE/BcK5GR6Ft6vTll8N4O3hxk52+HGuHm65I\nwuEv/xI2bw5dUiIiGaJgOFe5HCxcOLU6DJ0TQKHaAcH+pfC2D8BAj8PWreF9FA4i2RTH8LnPtdz/\nUQXDuYqiMIB8442wYAHRSBhgXvXMNMcZTLTBh65MWg7j42GKq4hkS3GCyS23hGULhYOCoR5RBHfc\nAQ8/DKtWEY3AlvvDLKQzthwMJix0P4lIRg0OwokTMDERlv39LRMOCoaZEEWwZQt0dhKNwL5dcONj\nsKA0IEoGo3FYWAjdTwD09MxyhUWkqjgO438DA8RLnc+9FeKlDg8+CO94x+S++dzFVNd9DGZ2DdAP\n/CGwyt0PJNu7gG8AbwJ2uftHqhzfD1wPPJ9s+m/ufl89dWqaKApTUPv7ifbsIRpx+g6Faap7Xh0G\npW0C3nQUVh6FvkMQjRAGsEdHm1x5ESl/csHAytDlWzBoL8DGIeg7lNysCuGpBp2dsG/fvLsnqd4b\n3A4D7wW2lW0/AdwC/FHycyZ/5e6311mPbIii0Nx89FE4fpxoBPrz8OirYNyhsxC6mqKRpLxZGMDO\n5ZpXZxEJobB6NYyNEXfD4ArY9u+SqecGpwy29sKdPaFHIBph6p6kwUEFQyl3fwrAzMq3vwB838x+\nv573n5OKg9LJbx7RyEn27g7jCbnhklBob4frr9fjMUSyYHBwMhTWXBduSi2GwqRkyvngipL/xwDH\njs1yZRsvC4/E+KiZ9QEHgI+5+68rFTKzTcAmgGXLls1i9c5BFIWfvj4YHCTasYPoyEQIg6veDYsW\nKRBEsiKO4eBBIPwCN1a8H8nPeNSURYsaVbOmmTYYzOwhoNKf/FPu/u06P/8O4DOEv4LPAJ8HPlip\noLsPAAMAvb29tf6VNVdpQOTzevy2SJaUjSnE3bB/SRhTSE0WKSqE7uC+QyXbOjvD/+95ZtpgcPe1\njfpwd3+2+NrMtgPfbdRnNVUxIEQkG0rGFCDcU7R6fWgtAKELyZnqSnK4/OdhzHBy0sjatWFMcR7+\n327qdFUzW1yyejVhMFtEpLGSMYXJ1RVJKFiFsh6ekjwZCsVJI/M0FKDOYDCzq81sBIiAe83sgZJ9\nw8AXgPVmNmJmr0u27zCz3qTYbWb2IzP7R2A18Gf11EdEpBHe8+MkFBYsgBtumPePzTf3udFdX6q3\nt9cPHDjQ7GqIyFwVx2HMb3w8rHZDbj2Mt5FqNZjDggLkd0G0ZFW4kXUOB4KZPebuvdOVy8KsJBGR\n2VW8IbVkWnl+V+k3L04tc8MQHW2Hu+Z2KJwNBYOItKayWYNRVxfRN78Jf78n3LxWZAY3XN8yoQAK\nBhFpdaWzBt/whvDkguK3M7a1hYHmeTgl9UwUDCIiRcUnF+Tz4VsaR0db8v4jBYOISCndd6THbouI\nSJqCQUREUhQMIiKSomAQEZEUBYOIiKQoGEREJEXBICIiKQoGERFJUTCIiEiKgkFERFIUDCIikqJg\nEBGRFAWDiIikKBhERCRFwSAiIin6PgZpnDhOf+FJ+bIFvwBFZC5QMEhjxDGsXh2+IrG4qTt82Xpu\nGKIjBh0dsGED9PTA0FAo1NensBBpsrqCwcyuAfqBPwRWufuBZPtlwK1AJzAOfMLdv1fh+JcD/xtY\nDgwDf+Luv66nTpIRg4OToRB3w+AK2NkDp9qgowAbhpyeoycZPbyVrv0wen4SGNu3w/XXh7BQ66Kx\nylt0OseSqLfFcBh4L7CtbPuvgP/o7s+Y2R8BDwBLKxx/M7DX3W81s5uT9f9aZ50kC558EgihkFsP\n4+3JdoOTBlt7y8o7tDt87P9OcOHhreS+C9FIsq/4hex79+rCdTbO1JU3NAQ7d8LJk6GsGXR2wr59\nOsdSXzC4+1MAZla+fahk9QngJWa20N3HSFsH5JLXu4E8Coa5L47h+98HQkthvB0wwEvKFNeLyzaY\ncLjtrdDmsHAC9u5OwqFQgPHxcJHTRas2cQxr1oRWW6EQLvzuU8ty7qHs4KDOsczKGMP7gIMVQgHg\nEnc/mrw+BlxS7U3MbBOwCWDZsmUzXkmZQfl8uBidrSQkCm0w7mE8Ihph6rfZXG5GqzmvlLYOhobg\n4EE4cWIqBNyTMR6n68XQdVdc5oZLWmci1BAMZvYQsKjCrk+5+7enOfb1wP8ALp/uc9zdzazCrzKT\n+weAAYDe3t6q5SQDuromX/Ydgjt7QqvBAQqkWw/lS8AK0FkIFywA1q2DP/9z/SZbTdlAf3GQv6tn\nKgCGFsPOlaEbj5IGvgGdE7BvVxIOL33p7NdfMmfaYHD3tefyxmbWDdwN9Ln7z6oUe9bMFrv7UTNb\nDDx3Lp8lGTM6OtllEY2Ei05+ebhAffN1sOfV4G2ku5KKFyuHNz0DW+4v+S120SKFwpmUDfSvXg9j\nbaQCYFJZF54bjLWHLr9oBLj9drjqKp3vFteQriQzuxC4F7jZ3f/+DEXvAa4jzGC6DjhjC0TmiFwO\nFiwI4wKEC07xIv+G5+DRV8EYUDCwibKQAFYeVdfGWTl2bPLl4IpwoT9tDKd8jKf8LS5IXhQKGmeQ\nuqerXg18GbgYuNfMHnf3dwIfAX4f+LSZfTopfrm7P2dmO4CtydTWW4G7zGwD8E/An9RTH8mIKAr9\n3YODYb14n8KOHUQjp9i7e6oFMXo+/GYhfP4tISgWFEL306TOznBvg1S3aKqnd/ICf7Zv8cIM1UXm\nBfNKMxQyrre31w8cONDsasjZqjR9cmgItm8nXjyRvvmtvR02btQNb7WI49BKGx9n85XJVODpWgzJ\nNvMQxvldSSutrS3MKNM5n5fM7DF3L58sfno5BYM0XRynWxe62ersxTHcdhvxgW+lxhgsGUcwh/YC\nbByCnqNVZiW1t8NXvgKbNjX3zyINU2sw6JEY0nxRpBCoVxTBokWnDfZXDAAzuOwy+NP36e5yqUjB\nIDLPlA72Tyre2Fa8i7y/XyEgVSkYROaLvj7YsQNOnQrr7e167pScEwWDyHwRRfDII1PjNRq4l3Ok\nYBCZTzReIzNA3+AmIiIpCgYREUlRMIiISIqCQUREUhQMIiKSomAQEZGUOfmsJDN7nvA01plwEeE7\nquca1Xt2qd6za67WG7Jd91e5+8XTFZqTwTCTzOxALQ+VyhrVe3ap3rNrrtYb5nbdi9SVJCIiKQoG\nERFJUTDAQLMrcI5U79mles+uuVpvmNt1BzTGICIiZdRiEBGRlJYJBjMbNrMfmdnjZnba94KaWc7M\n/iXZ/7iZfboZ9SxnZhea2TfM7Mdm9pSZRWX7zcy+ZGZPm9k/mtnKZtW1VA31ztz5NrM/KKnP42b2\nWzO7qaxM5s53jfXO3PkGMLM/M7MnzOywmf2NmZ1Xtj9z5xtqqncmz3fN3L0lfoBh4KIz7M8B3212\nPSvUazewMXndCVxYtv9dwN8Rvu79zcAPml3nGuudyfNdUr924Bhh3nfmz3cN9c7c+QaWAr8AXpKs\n3wWsz/r5rrHemTvfZ/PTMi2GucjMfhd4O7ATwN3H3f03ZcXWAYMe/ANwoZktnuWqptRY76xbA/zM\n3ctvpMzc+S5Trd5Z1QG8xMw6gPOBZ8r2Z/V8T1fvOa2VgsGBh8zsMTPbVKXMW5Lm6t+Z2etns3JV\nXAo8D9xpZkNmtsPMLigrsxT4Zcn6SLKtmWqpN2TvfJe6FvibCtuzeL5LVas3ZOx8u/sR4Hbgn4Gj\nwL+4+4NlxTJ3vmusN2TsfJ+NVgqGt7r7G4E/Bj5sZm8v238QWObu/xb4MvCt2a5gBR3ASuAOd+8B\nXgBubm6ValJLvbN4vgEws07gPcD/aXZdzsY09c7c+TazlxFaBJcCS4ALzOy/NLdW06ux3pk732ej\nZYIhSXnc/TngbmBV2f7fuvu/Jq/vAxaY2UWzXtG0EWDE3X+QrH+DcMEtdQR4Zcl6d7Ktmaatd0bP\nd9EfAwfd/dkK+7J4vouq1juj53st8At3f97dTwJ/C7ylrEwWz/e09c7o+a5ZSwSDmV1gZr9TfA1c\nDhwuK7PIzCx5vYpwbkZnu66l3P0Y8Esz+4Nk0xrgybJi9wB9yeyNNxOatUdns57laql3Fs93ifdT\nvTsmc+e7RNV6Z/R8/zPwZjM7P6nbGuCpsjJZPN/T1juj57tmHc2uwCy5BLg7+XvqAL7m7veb2Y0A\n7r4V+E/AZjM7BRwHrvVkekGTfRT466Sb4OfAB8rqfR9h5sbTwIvAB5pV0TLT1TuT5zv5xeEy4IaS\nbZk/3zXUO3Pn291/YGbfIHS7nAKGgIGsn+8a65258302dOeziIiktERXkoiI1E7BICIiKQoGERFJ\nUTCIiEiKgkFERFIUDCIikqJgEBGRFAWDiIik/H8/odn77WYsTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1113bc710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "#telra and teldec = 7.11 -10.53\n",
    "#tileid is 23658\n",
    "mindec = -12.5\n",
    "maxdec = -7\n",
    "minra = 5.4\n",
    "maxra = 8.9\n",
    "ra = []\n",
    "dec = []\n",
    "import desimodel.focalplane\n",
    "import desimodel.footprint\n",
    "import desimodel.io\n",
    "tiles = desimodel.io.load_tiles()\n",
    "\n",
    "while mindec < maxdec:\n",
    "    startra = minra\n",
    "    while startra < maxra:\n",
    "        ra.append(startra)\n",
    "        dec.append(mindec)\n",
    "        startra += .02\n",
    "    mindec += .02\n",
    "\n",
    "ra = np.asarray(ra)\n",
    "dec = np.asarray(dec)\n",
    "\n",
    "points = desimodel.footprint.find_points_in_tiles(tiles[8254], ra, dec, 1.651)\n",
    "\n",
    "#plt.plot(ra[points], dec[points], 'b.')\n",
    "targetindices, gfaindices = on_gfa(7.11, -10.53, ra, dec, 100)\n",
    "plt.plot(ra[targetindices], dec[targetindices], 'r.')\n",
    "targetindices, gfaindices = on_gfa(7.11, -10.53, ra, dec, 0)\n",
    "plt.plot(ra[targetindices], dec[targetindices], 'g.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "([], [])\n"
     ]
    }
   ],
   "source": [
    "def on_tile_gfa(tileid, targets, buffer_arcsec = 100):\n",
    "    \"\"\"\n",
    "    This function takes a tileid, a table of targets, and an optional\n",
    "    buffer_arcsec parameter to return the indices of targets lying on the GFA\n",
    "    as well as the GFA locations from 0-9\n",
    "    Parameters:\n",
    "    tileid: (int) DESI tile ID, used to lookup telescope (RA, dec)\n",
    "    targets: table with columns RA, DEC\n",
    "    Options:\n",
    "    buffer_arcsec: (float) additional buffer region around GFA to include\n",
    "    Returns: \n",
    "    targetindices: list of indices for targets that are covered by GFA number \n",
    "    in corresponding gfaindices \n",
    "    gfaindices: list of indices corresponding to 0-9 GFA location\n",
    "    \"\"\"\n",
    "    import desimodel.footprint\n",
    "    telra, teldec = desimodel.footprint.get_tile_radec(tileid)\n",
    "    return on_gfa(telra, teldec, targets['RA'], targets['DEC'], buffer_arcsec)\n",
    "\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "sweepname = 'sweep-280p050-290p055.fits'\n",
    "sweep = Table.read(sweepname)\n",
    "import desimodel.io\n",
    "tiles = desimodel.io.load_tiles()\n",
    "print(on_tile_gfa(23658, sweep, 120))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RELEASE BRICKID BRICKNAME OBJID ... SHAPEEXP_E2 SHAPEEXP_E2_IVAR TILEID GFA_LOC\n",
      "------- ------- --------- ----- ... ----------- ---------------- ------ -------\n",
      "   4000  585540  2814p502   398 ...         0.0              0.0   2199       8\n",
      "   4000  584615  2816p500  3947 ...         0.0              0.0   2199       8\n",
      "   4000  584615  2816p500  4206 ...         0.0              0.0   2199       8\n",
      "   4000  584615  2816p500  4447 ...         0.0              0.0   2199       8\n",
      "   4000  584615  2816p500  3687 ...         0.0              0.0   2199       8\n",
      "   4000  584615  2816p500  4046 ...    0.127714          9212.78   2199       8\n",
      "   4000  584615  2816p500  4047 ...  -0.0190994      1.01114e+06   2199       8\n",
      "   4000  585540  2814p502   563 ...         0.0              0.0   2199       8\n",
      "   4000  587373  2802p507  1302 ...         0.0              0.0   2199       9\n",
      "   4000  587373  2802p507  1368 ...         0.0              0.0   2199       9\n",
      "    ...     ...       ...   ... ...         ...              ...    ...     ...\n",
      "   4000  597979  2807p537  1112 ...         0.0              0.0  42563       9\n",
      "   4000  597979  2807p537  2242 ...         0.0              0.0  42563       9\n",
      "   4000  597979  2807p537  1391 ...         0.0              0.0  42563       9\n",
      "   4000  597979  2807p537  1490 ...         0.0              0.0  42563       9\n",
      "   4000  600515  2800p545  3009 ...  -0.0156742      8.62941e+06  42564       7\n",
      "   4000  600515  2800p545  3475 ...  -0.0459809      1.61356e+07  42564       7\n",
      "   4000  600515  2800p545  1624 ...         0.0              0.0  42564       7\n",
      "   4000  600515  2800p545  2433 ...         0.0              0.0  42564       7\n",
      "   4000  600515  2800p545  1568 ...  -0.0061286      2.08762e+06  42564       7\n",
      "   4000  600515  2800p545  3632 ...   -0.019829      1.62225e+06  42564       7\n",
      "   4000  600515  2800p545  1152 ...         0.0              0.0  42564       7\n",
      "Length = 1005 rows\n"
     ]
    }
   ],
   "source": [
    "def get_gfa_targets(targets, rfluxlim = 1000, tiles = None, buffer_arcsec = 100):\n",
    "    \"\"\"\n",
    "    This function takes a table of targets, as well as optional parameters \n",
    "    including a minimum flux in the r-band, a list of tiles, and a buffer in arcseconds\n",
    "    and returns a table of targets on the GFA satisfying a minimum flux_r\n",
    "    Parameters\n",
    "    ----------\n",
    "    targets: table with columns RA, DEC, FLUX_R\n",
    "    Options:\n",
    "    rfluxlim: (float) r-band flux limit; default 1000 = rmag 15\n",
    "    tiles: table of tiles, default to desimodel.io.load_tiles()\n",
    "    buffer_arcsec: (float) additional buffer region around GFA to include\n",
    "    Returns subset of input `targets` with additional columns:\n",
    "    TILEID: (integer) DESI tile ID\n",
    "    GFA_LOC: (integer) GFA location [0-9]\n",
    "    Note that the same target could be repeated with different TILEID, GFA_LOC\n",
    "    \"\"\"\n",
    "    if(tiles is None):\n",
    "        import desimodel.io\n",
    "        tiles = desimodel.io.load_tiles()\n",
    "    import desimodel.footprint\n",
    "    points = desimodel.footprint.find_points_in_tiles(tiles, targets['RA'], targets['DEC'])\n",
    "    alltargetindices = []\n",
    "    tileidlist = []\n",
    "    gfaidlist = []\n",
    "    # Checks if the flux_r meets a minimum threshold\n",
    "    brightindices = np.where(targets['FLUX_R'] > rfluxlim)\n",
    "    #print(brightindices)\n",
    "    #print(len(brightindices))\n",
    "    if(len(brightindices[0]) == 0):\n",
    "        return np.array([])\n",
    "    counter = 0\n",
    "    for lists in points:\n",
    "        if lists:\n",
    "            tileid = tiles[counter]['TILEID']\n",
    "            targetindices, gfaindices = on_tile_gfa(tileid, targets[brightindices[0]], buffer_arcsec)\n",
    "            tileidlist.extend([tileid] * len(targetindices))\n",
    "            alltargetindices.extend(targetindices)\n",
    "            gfaidlist.extend(gfaindices)\n",
    "        counter += 1\n",
    "    validtargets = targets[brightindices[0]][alltargetindices]\n",
    "    tileidlist = np.asarray(tileidlist)\n",
    "    gfaidlist = np.asarray(gfaidlist)\n",
    "    validtargets['TILEID'] = tileidlist\n",
    "    validtargets['GFA_LOC'] = gfaidlist\n",
    "    return validtargets\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "sweepname = 'sweep-280p050-290p055.fits'\n",
    "sweep = Table.read(sweepname)\n",
    "\n",
    "import desimodel.io\n",
    "tiles = desimodel.io.load_tiles()\n",
    "valid = get_gfa_targets(sweep)  \n",
    "print(valid)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "tiles = np.zeros((4,), dtype=[('TILEID', 'i2'), ('RA', 'f8'), ('DEC', 'f8'), ('FLUX_R', 'f8')])\n",
    "tiles['TILEID'] = np.array([23658] * 4)\n",
    "tiles['RA'] = [0.0, 1.0, 2.0, 3.0]\n",
    "tiles['DEC'] = [-2.0, -1.0, 1.0, 2.0]\n",
    "tiles['FLUX_R'] = [50, 50, 50, 100]\n",
    "targets = get_gfa_targets(tiles)\n",
    "print(targets)\n",
    "print(targets.dtype.names)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvX2cFNWZ9/093fMCGg3rBIM6ImLME9mbBcQPphSxCC6G\nJN6ik80TNY5RNwPxJSF7G5AYdnmWhAlo9iYrMTIGXUbRfZJnhNtsHGXttYQ4lVUQDCqJ8WXE4cWX\nSYgkwrx0n+ePUzVV3dPD9Lx2V8/1/Xz6093V1TPn6qr61TnXOdd1Ka01giAIQvSJ5bsBgiAIwuAg\ngi4IglAkiKALgiAUCSLogiAIRYIIuiAIQpEggi4IglAkiKALgiAUCSLogiAIRYIIuiAIQpFQMpz/\n7GMf+5ieMGHCcP5LQRCEyLNjx473tdZje9tvWAV9woQJbN++fTj/pSAIQuRRSr2Vy37ichEEQSgS\nRNAFQRCKBBF0QRCEIkEEXRAEoUjISdCVUs1Kqd1KqV1Kqe3etuVKqX3etl1Kqc8NbVMFQRCEY9GX\nVS6ztdbvZ2z731rruwazQYIgCEL/GPEuF9eF2lrzLAiC0FeWrHU5+4ZalqztLiJ1jS6Xfq+Wusbh\nEZhce+gaeEoplQTWaa3rvO23KqWqge3A/9Ja/3EoGjlUuC7YNnR0QGkpOA5YVr5bJQhCVFiy1mX1\nwTlwejurD5bB2gTzp1s4Dhw6wfss3s6WZ8tobEyw+CprSDUm1x76TK31VGAecLNSahbwE2AiMBU4\nAPww2xeVUjVKqe1Kqe3vvffeYLR5QIR75PX10N4OWpvnRYukpy4IQu7c99/1ED8KsSTE2rnncYeL\nLoI77oA7f+5AvL3rs807HWbPHlqNyUnQtdb7vOd3gU3ADK31O1rrpNY6BdwHzOjhu3Va6/O01ueN\nHdtr5OqQ4rowe7b5sWfPhoMH0z9//nmYM0dEXRCE3lmy1uWPEx4ApY0PQ8f5826bZBL0aS76hL2Q\nKoFkHFJl0GzT1mY6kkNFry4XpdTxQExrfdh7PRf4Z6XUKVrrA95uVwAvDV0zB05do8uKBx3axtoA\ntE1wePlPNhCMf/yeurheokldo0vDDoeq6TY184rnALquOSdtuzjPy6gcN/84VFRAays8vNuBsztB\nASkFO2+AFgsqXbjOuFpIxeGFr8GL1eazISYXH/rHgU1KKX//h7XWTyilHlRKTcXcm5qBBUPWygFS\n1+iy4Nk5cHY7nBUHFMQ6+X2yDH6f6PqhlYKyMnPhCNHBdWH1Iy6bTwz8lZAoaHHIFdc1o8b2dnNu\nJhLFJepd12a8jS1ujK2//TEPfasm383qhn8c2toglYJYDKi0YWIZ6HbTA3+x2uw8wQlcLRr403ji\nByxSyszVVVcPXTt7FXSt9RvAlCzbrx2SFg0BDTuc4AcmZTbGtDkQExxosSgthRtvND92MV0wxY5/\noR05z4HZ/kXUTsMOJ/KC7rqwfHkgIsUweswcbZhrsw1iKdApNv7xFk5bO5lVtxSWkY5jfv+UJx+p\nFLDXgg0JoyHNdlfHsDJp05IMhL50v83ae0yvfqhHWcOabTFfVE23Ta9Nt4P2eui6s8uvpZQR85/8\nJN8tFfpKfT0cOQK8acOs4CKqmm5H2lWRrUcY9dGjP4fljzaeftq7Nt0Y6JRxXagkq3/mcFaZRU0B\nddRt27TZPx5dtFiofRZaB5u+MMVi2oUJ1iccTm23WfzQ0K5sCTMiBN301BJdfjqA9QmHnZtsUgcs\nykYN7TBIGBpcF+67z3vTYnpLJ0x2+PrnbCaPsSLrqsjsmcdicMklZptvQ1T8zmHq641NYJ5Xr4ZN\nmywatv2YLSU3mV56qhTetFn4fZf1rzrcOKcw7LMscw75PvSdO+H++yGZNJ/7zz4186y8tHtECDp0\n/4Fr5lm4VdHtwQnm2KVdSC0Wh1ssfvRf8MH1RsyTSfN8+1qXN+52mKhsfnDL8PWY+kpmz5zTXTjT\ngekVLE+0UnXIBvD8ztGeL/jFL4y9x/15Mnw0bnrpaDh5N3reIp6Lt/NcAdlnWek6UV0dCPyttwbx\nLPnsHI4YQc9G5gESooVtm95r2hAYI+Bgeubt7UCly9YzjQC2JMu46OoE2x4uTFFP89We7kL1HFLx\nNraoFHTE2NJUyvGHp8BHfb9zdOYLqqvNiMq/CWtt7H21zfHWamvQSZjUEJpULFz7wvoxeXJhdA5H\nfOi/EF0sC847r/v2eBxe+cDl41+s5ZNzXCqmO2kBHslKB8cZ7tbmhu+rjcchdqbjtdufyE9BvI2/\nfPQ5UClIxbrmC6KAZcE995hebCwG5eXw8svwSqMNybKu9dqnHa5Kez91jJ3vpveKZcHSpfnvII7o\nHroQfW68EZ57Ln1b57igR874MtQTa+CcYMI03mIX7ORi2Fd76ATbhJNrvzeuAG26YckYJ31wCbWX\nLi/I3mtP1NQEvdlDh4wfHYLVIvOn2cw43+KOeyejxzuovTZjFkbHvnwjgi5EGn8lxJ13wuuvm2F8\n+jrgdvRxrcx6M8Eb2vOhF6i7xScYyluc1Wgm88ceX8HOAzt5pfx+UElIlUVOzH18+y69NLSxxSK2\n32LxavN21AqL9n1W5Ff2DDci6ELk8Xt9XZOJzbYZsns98pIWmx+ssrAKWcV7IHMyv66xOnKrW3qi\nqgq2bAne33Zb4LLwRyn59klHDaXDCyiHmPPOO09v37592P6fMLIIh2Y3vuTyapvDJ8vtIc9wJ/Sf\nujpoaDDiXkjrzgsNpdQOrXWWGaOM/aIs6FEOHBEEQciVXAU9si6Xujq45RazBKq8PFqBI4IgCENB\nZAQ93BsHuPlm6Ow0r9vaop/jQsidKEZJCsJwEAlBr6szAp5Kmd74ddelB5PE4zITPlIIsvNlj5IU\nN5wwkil4QXdd41oJ98bBCHtbmwlQWLtWLt6RQlrmTN3O+oRD6y6rK7/G+vXmXCkpkeyZwsij4AU9\nM1+HUjBtWvBeLtiRRVrmzFQZOzfZPPdG9/06OmDdOtiwQeZXBkKmq1NGP4VNwQu6bZtQYb9nDkGP\nPd+JcIYacR90p2aexdbfrqHxrQbGH67ixTe9CjEZOakhqEB100qXN1IOE2M293xHljBm0tN5Fi6u\nEY+bzlRnZ/SyV44kCl7QLQuuv970trQ2vXW/x+5n0RuViN4EWfgi2n3ITPJNHWMz5rDV1RuKavpX\nCNYXT50KY8aki8VAblR1jS4bWxfBR9v5w0e2wbnAZxcZN0yyzISQh0X9NJddU4zPfVeyjIuuSrDt\nERF1n2NVRHIcOHo0uO6UCm6S9fVB+P6uXYW1jtwvAA8jbwRf8IIO5qA88EB6Lx0Isuh1RCuNaPgi\nUuNdOq/2JvkOlkF9gtErLK67Lj39a5RW8dTVwQKvIOGWV0zvueQ+m60bjQG2HaQa7atdmT70sqkN\ntIfe+xWoAGbMgN+Odfgg9HnydAfHEUH3CYt25mqxQ4cICjdUuugzHdRbNiXvWNx/vzmG/udbtpjU\nC6tW5cGIEK5rzi8/4+b69fDMM9G5dgZKJLIt+r10U9Y0RDhnR8yk2YwCforUZBI6T3PSbOAMp+vG\n5Wfdi1I+C9c1eVWAoFjuZ5bRefUcbl/rUl9vbA/39PpC1XQ7LRPfJJ2emY9mGzCT5WvWwGWT0/eP\nvz3wxFyuC7W15jnqhEU7lTLvfXbt8l74x3H2MvS1c/irv3Hp7ITMmMS77sr/b+I45kbj09EBixYF\n7SqmY5eNSPTQwfTS169PP1jxt22SoZwdUUkjattGqJNJjACFbKDZJpUyE79+Av2o+ND9YK+uY5SR\nJGtbi8NJJw3MkMzqU5PHWFx01WSSpztdPvRYzJQTNEmgLFia4Bc7PR/6ANwtrmuyA/7iF0bMiiGg\nrUu0PX76UzjrLOM+6cq1knEcD5Y7kMpudL5Hkn6O/PBCiuefN9s//Wl49llz7HSlS9nZDhdPsHny\npxE+gJlorYftMX36dN1f1q3T2jsUGrSOxcy2dY836bkrVup1jzf1+2/ng/nzQ/ZUNmlmrjTP3rb5\n87VuatJ65UrzXOg0NWldUpJ+jKhs0twxWrMsbp5Pb9KVlWY/pbQuLx8c25qatF640PxmCxcOze+1\nbp3W8Xi6fUqZ4xNlMq8r/zFrlvktTzwxy3H0ztO5c7W+5hpzLcZiWo8eXRjnarZjlfW8/Edjz9wb\nC6DRvQBs1zlobGRyuZx/fve81/Nvddlf5hVijVACprCP+VjE4+Y5CpOitbVwxx3dh+E9rUD5xCeM\nu6WQbfJxXZg1K4iFCLN4sZn0ragYnqruQ0FdHXz72/DBB8fYKeM4lpebIs+WVZirsfyJ0QceCFx8\nXcyshc8sMyOOZJzR/72CD59cmre25kLR5XIZNSpjw+kum0/0ihiMKuOXX0nwzDBW1+4vrgtf/3pu\n+4ZX8+R7KNsb2Ya6gBHxlu4Nf+012Ly5sG3ycZzuZe7AzOn8y7+Yz/xizlF0w0yeDEeO9LKTdxwn\nTIDPLkxfPVKIpRz9Np14YmhOxyfDzXnR6XYeWjg0RGJS1HXh178O3isFk+Y5aZOJHacVblmxMD2J\nQ08oFY1JUcuCyy7r+fNYljPt0UeHrj2DiW0boc6clFcqEHMwz/7NN0pknpPdFh94xOPw8MPB/ESh\n47rmhttt1NhiKiSV/GoFcw8misqHHglBD0eLKmXcFd+cb6evbnjTpqIin63MDds2YelhpkwJ3CtK\nBY+yMmNrVHp8ixebNitlBLykJHj+yU9g7tz0/a+8Mj/t7Ct+WbgFC4yw+zbddlvwHsxzFG6+mYTr\nmI4eDffeC/PnB8FE8bh5v21bNM5Dn8wbVTxujlc8DqNbLbZ+f2lRiTlEJB96T8EPV3zDZfNOB5pt\nYvstvvc9U6i10PFXS+zfb/KN1NQUT4h1b3YsWWJ65ldemf81y/0h018cLqoRVR86ZPeDF6JvvC/4\nuuHnfPrxj4N6plGzqegKXPR0wkU5mlIQhKEl6jcln6IT9J4olgMmCILQE0W3yqUnCnGGXRAEIR9E\nYlJUEARB6B0RdEEQhCIhJ0FXSjUrpXYrpXYppbZnfPa/lFJaKfWxoWmiIAiCkAt98aHP1lq/H96g\nlDodmAvsHdRWCYIgCH1moC6X/w0sBoZvqYwgCIKQlVwFXQNPKaV2KKVqAJRSlwP7tNYvDlnrBEEQ\nhJzJ1eUyU2u9Tyl1MvCfSqnfAt/BuFuOiXcDqAEYP358vxsqCIIgHJucBF1rvc97flcptQm4GDgT\neFGZTD6VwAtKqRla64MZ360D6sAEFg1i2wVhUOkprD8ctObXSi2kGpp9JVxzc9q0aKcs6ImRGnDY\nq6ArpY4HYlrrw97rucA/a61PDu3TDJyXOWkqCFEhnEaipATmzYNf/jKocv/007B8vcuW3znQbLNl\ni1GJqIh6OOfMN76RXp9XKZOeulhSZ2SmBFmzpjhvWtnIpYf+cWCT1xMvAR7WWj8xpK0ShGEmXOc1\nmTS52n3a2uCmlS67psyB09pNls8NCRoarEgIuuvC7NnGvmw56/VpLkcmONy43Gb98sKvKRAmW088\nfCzb2oL6A1HMVd9XehV0rfUbwJRe9pkwWA0ShHzgp5DtqdDD7zuctLqaTHCoqoqGMtTXBz3ybgVI\n/ALQ8Xb2JMuYdU2CrRujIeqZo6rrrzeFN/xj2daWnj63ra3wC8UMFIkUFQTMRb5mTfbiDvE4nF1q\np+XfP7vUjkTvHODgwWN8GC4AHWunszIahWKge0983Toj8GCOZWbewXg8ernq+0rkk3Plg5E64VLs\ntLYaQdenBfUz4wcs7rkHJk+2mPnlBKnxDrG9Nhv+PToHfty4Y3yYUY6tpMWOjOj5PfGjR4MS0H4v\nfO/edEFXCtauLf7rVQS9j4ST5sfj5iSJSk9NODa2DfEJLqmrvFq1yTK+nEzQ2mqxezeUvmPR/rZF\naVm+W9o3pk3r+TO1zyL+cIJP/q3DJ8ttFkfE3QJBJanVq4M5j1TKTPzuzYhdv/zykXGdiqD3AdeF\n5cuDHkEqBTffbKqgQPFVfBlpWBacfYnDKyFf+cPbHGLPWl31Q7U2K1+i5IttbTWToZm1bONx+NrX\noLrawoqKMRlYFsyYAY89FhTqbm01vvT774eODigtNeURRwIi6DlSV2dmyzMvilTKTDpt2JC+TGrn\nTnjgAXPxq/Eu517pcOMcm5p50bxwRgqfLLd5JeSC0G/YJJNGKPwam1GrG+oXufY7ImDs+NrXTK3X\nqOPb519/fgfKcUZehyryFYuGA9eFmTO7izmY2fWPf9zUB9U6uPA7O72LJ7SKgGQZ6y5MiKgXMK4L\ns65x6ax0iO+1KTloda1Fj/J6Zj+YyO9kFFvJxmIfDY+YikXDQWb18DCdnbBvn3kdiwXrfLvuk+FV\nBLqdhh2OCHoBY1mwdaOF41iRL9gdxq/sVV1dHPZkIpXLDEUv6IMRqm3bpife2Rlsi3kLPsNCP3Ei\nfPvbcOutZvgHdFtFUDXd7l8jhGEjUxyKSShE+Iqbohb0ujpYsMC83rIFXn8dVq3q+/DMsmDrVjOb\nvn+/+d6YMXDokNnm8+1vm5vGzp1w773exhaL2IMJzvs78aELgjC0FLWgNzSkv7/rLjjrLFi0KJhA\nydWPaFmwaVP37Wed1X0EUF1tJknb2kxP/sd3WNTUiJALgjC0FLWgV1WZnrmP1kZ8/eiy9vaBLz+r\nqenuyvHXxxajr1IQhMKlqAW9psa4We6804h5aakR+W3b0pc4QeBrnzrVuFMGKsTiqxT6SjGu1Ahn\neYzqCqEoUdSCDjB/PvzoR0bAlTJBQLeudnn0BYcrz7XZvdvixhthz2ET7r3lYRtaLMrKhid4ZLgu\n4mIUCzA34vXr4dRTTfBIVG3zI5CPfMxFNTpcfaHNQ7VW2udRO351dXDTTUFCsFhsZGQ8zCta62F7\nTJ8+XQ+UpiatV640z7mwcqXW8bjJ9BCPaz3/1ibNHaM1/xg3z5VN5pG5Da1nzMj9//TXltGjTbtG\njx66/5X5f9at69tvWGj458DixX4GD/MoLY2mTU1NWs+dq7udh9fc3tT1eenEJs1FK3XpxKZI2Jh5\nbPxHLGZsDdvQ12t6JAJs1zlobKR66JmJ63O50/sJfPzvvNrmdEuDCnTf1mLx/PPm/w2kRxEecm58\nxmXPUYdzRtn84BaL1auDdK1tbSatwPLlA++9ZPbmMrPS3XyzubyiElwS/g0bG+EXv+ieSQ9MmPft\na11GJRyqpkdjRZF/Th89ClzopJ2HG7c6zKqzaHzJpePq2RBrpyNVxupHnmZTAR+0JUvSV3+FSaXg\nqaeM2zORMNv6ek0LPRMpQQ8LU64TmpkTlKsfsXlFx0AngZhZJw5pa8X9bVqbC62+PveTbMnawJ0z\nf7rVlcgrdWoQMfpesowLv5xA7w3+aCplJnCfecZUxxnIDSTzAunKSjfWRZ/h0Nlsw9vWoEwKDyV1\njS7rEw47N9l0jNkNkxpgbxUkewgoqHTZeuYc6Ghny7NlQOFH5TqOOT+0plvMAm/a3HwzjP1qPZzU\nBgpQbexK1QOFZVd4DurOO4+9byoVXL8wuIsURjqREvTM3nau+TTCE5S/u2c3nNRh3sQ6+MjE3fx5\naw1sSHSlTKUlOKO0NuHSudReXLLWZfU7s+H0dla/U8aWlU/T3m6Z4KOMiFE93oG9lkkNMMUr8Phi\nNW0tVp9uIJlku+ktXWrmDVYfDFIQqAcTlL1n5fwbDjd1jS4Lnp0Dx7fD1TGIe8fsrNCypUkNnEMV\nh50aWlro9hvf+XOHyWMKO3ugbYcSZ7VY3c7DVAw6M4pSjD5u+Nt5LDLjPXojMx9Of65pITuRKHDh\nulBba14nErBiRf+HZntOXm5eeIUM/nzeckpKMBfThxVgL4dz69K+09EBC7/v8p3Ha7nwyy5XXGHa\nlMkDL9RDvA1iGuJt7NL1lJR4UaV+78srkECzbcT8qzacd695fHW22TYA/JtePJ5+gfzHbietkMFx\nkxzWrDGf1dZmtyefNOxwQu31xNwvPvHpNXDZAjhrC3vOWsD+T3+F0lJQe+203/i1p2wuvtgITiHz\nhS8Eib8yKS+H68+tNnalFCTLWDS7evgbeQwy4z2yUunCzFqodLn88uD69UfQA7mmhYCC76FncyEs\nXTqAP1j+p27vzz4bDn+yjpapXjfD7wW+UEMsBrrSRV9rerc6WcbmDQkaZ1vdXCPHHwfvZfy7KVPg\n3HPhxBMt7nk0wZ8rnGAUMLPWiJV/IcfaURMdqqv7f1aHc0Tv3w+7d5ttutmGscFw/i8v2dz6X0ZE\nCjFZU9V027hNdDtor4fu+82194MpQEPqrzeSemMWsz5SA28meEM7tPzKhrctOghSHBeKbT7hc5tK\nl9jUepKTH4BYJyTLGNuY4P+sNaltz2p0aNhRmHMDmfEe3chIUMcHibR0vbLEd/AoeEHvj9/8WJzR\ndgVvlW0MxGHPFezZA6Mu9LoZnkiUTmngW5fU4Djwm486HM2YNG3fZ3Vry9LPV7Og6X4j0qlSeLGa\n7fth1y4jnB0dFqRCX3jLNvspL/FLqoxv/5094JN79+4g4f9zz5nnRVUWC1akD+c7PF3UuvD8lzXz\nLNasSbDnQ8e092TPh/5KldnhsgXmGHrHi0kN/OrhGpSyiMct4knwPRWpVGHZ5tN1bp/iwlfmQPwo\nKO3Z1M7Y85wu4auZZxWckPv48R533dVDErsMV9j+ModCmwMoFgpe0PvrN++J5h8+xKiroa2yEX4/\nDzY/BMDRnVVw2pYuoe/8TRU/3OWtoa20YUb6pGm2tpgLzmF9wuHwbpvf7Tf+8w7PY+Cn173kEtOr\naW21ePkDhy376jn+OHNDGIyLNnMI3NBgVs+UHLTo9OYHYjGTcCzcQy80/+WiKosFC7zfo8XilAM1\nHDjgfTh+K0wJ3ZgPTCV1Qa0R/wMWl10G//EfRmDKywvPNggVpj7T8QRPG3tSClJlfHO+nd8G9oFV\nq0zMR7h6UBcZk703zrHz0MKRQSTyoQ92UEVtLXznO+nblAI9rS7oBb6QsZKi0uXEv3GYGLP5dKVF\ndfWx2xIeTvv+0eFybYQnqcAUz21thWXL6CrWcMklRuShsANWwtkyId2uM66s4+0TGih9fyptU+/u\nGtKXPpLgmYeMMYVsG5jzZPUjLptPnAOxdtBxJrXdwDdnD87NPR/4wV47dgRBRfEJLtOrJEFdf8k1\nH3okBH2wCa/9VQquugoefdRbXuiVsVIqOBnBiPK2bX0ThvCNCIZXXDLTBvdnDX8hki0d8qXfq2VL\nxzIzpE/GmfHhCv77roFMtAw/dY1uwfrI+4tfVAPotQMkHBsR9F7I7PVn5pyoqIBvfCMoBn3PPdEv\nMhvF8PFc6FriGDNDeqkKJRQbIuiDQLEKYDFSjD1cQfARQRcEQSgSchX0SAQWDYS6RpdLv1dLXWPh\nRM74gVKFFswjCEK0Kfhli30l7CbZfcjzrcZNbo/XX0+w6pbhG45nc9lkrn654QaZMBIEYXCIhKDn\n6ssOi6VSUD7HgfO9gAaOcueT9cyfPjy5PXpaVRIOlEomzZLC+++PprDLHIMgFBYFL+hhYSwpgeuv\n71n4Vj/icuQ8B/5SAce10vleBaRKQCVBafTU+1n9SDXHPwaNexzmnZNeRGAw8bPopVLm2Y9U7Mp8\neDTIEt3eboR9w4boLCcslmWQglBMFLwPPTOX97p1Rkgy/c91jV5wxuzvmrDwz3wX5i2CV+eZ3B8K\nUEm2/amejfE5/GHKMjbG5/CVpUPjyK6oCMKgUynzHozorVkDn/pU+v5aB8I/lCxZ63L2DbUsWTsw\nu7OlZMhE5goEYXgp+B66bZueuR/ko7XJ6119n8O3DwVL1ILsfCkTPh1LmVDjv4yD5KiusOOjR/Ey\nIqaANn6x06G2Nkgj21MgRF/dC62tQVrUWMy89//OokVeQYMM4vHBD1EPt3vzDi+F7untrD5YBmv7\nP6fQU0qGS79Tx7bWBibHq9j9bzXSgx9h5DOYTshR0JVSzcBhTL6jTq31eUqpFcDlQAp4F/iq1nr/\nUDRSay80X4MabzIfvhZvZ0GoiEGQnc8T61QMUmWo3dXoA9NgUgMzTqjicBnsUZ7oqxQfvFPBd79r\nbhpaB3lX7rsvCCbqb6Wk8vLuguf3bDNXi8bjsHbt4J74me2umO/A2UGSpEdfcFiF1S9feGbhEMuC\nSdfWseesBXAKPMcWOAd4oSZr4q++zIsUu0AUy1xEPtNdCIa+9NBna63fD72/U2u9DEAp9Q3gH4GF\ng9k4MCd6MhkkthozzeEPocxtDTucUCa6BA07HMYeX8F7f2ml6tM2W5OwMb4I4u08l9zGXx2+zqRj\njaUgGYPRrV0JtMIim0zCLbeYtKv19YHPO+xe8FPU3nhj9yjSbIIH6T3beBw+9zkYNy63CdG+FkQO\n+/GPHoX42zZMDJIknf9xmyuuCEq6lZSYydlcinn4Nvqf19XBHtIzVjKpAbWzhpIS87fC0biLFvV+\ng4yqQPRU6T5b8JPrwuzZwW8xkGpV+SbshvPdjYWYybOY6bfLRWv9Qejt8QR57waVzKH9319iG3eB\nJ0pTx9hd+/rC7l9Qk8fAne/UwumhSkEAyfJu5ebAiEY4f0syacT8gQcCsY/HwXnd5TuPO11paP0U\ntdlEPfMk7kno6+pMsqyqquAmcvBgIPa7d6cnpvrlL025umNdJGE/vtaw91mL+NsJTrnA4eIzbH72\nQ4uOyXVwlUlI1v5CDffea/bva4X2hgbg/SqTS74rNbHJqNXZaW56r71m2qOUefZvNKtXw4wZcOiQ\nSTU8dSqMGQN792YXiL6WBRxseho1QHC+dHQE7rZ4HCqmuhyc6y2hbYpT+8sbmBqrhrct2trMd9va\n4PbbjX253rTzRbZRRWZnpZAzeRYruQq6Bp5SSiWBdVrrOgCl1PeBauBPwOxsX1RK1QA1AOPHj+9z\nA30B9AXu1YRF7IUEqfEONNv86D2Ls8pg506z/7RpQe8vdobLxz6zF1JxLy1pGTXnV7NvXzUbtzrw\npp1Wbk4m9VaFAAAgAElEQVSpwLUDUFpqnjs7g/ac9DcuW8bNgdO8ZP0bEtBi0dCQe66XTKHPLOEV\nbgMYgZgyJf1vdHT03usJ+/HB+5stFjdNtti7FyPml3Uv6gHpdR9zEZWqKtiywPsBJjUw5kAVf3qh\nBq2NIO/ZE+zrJz/z27R5c3rKVb9YQmlp4Arz2+R/Z/364V3mmW10ERateDzdZefj37gOljtBTnCV\npLliHc3JDagXEoRzg2/dGnw3l5t2PshceTZvXtDxWLMmGEXOmxdcl8IwobXu9QGc5j2fDLwIzMr4\nfCnw//T2d6ZPn677Q1OT1mVl/iK/7o94PP11LKY1lU2aO0Zr/jGu+W65Hnv9Qr347qa0v7lwodYz\nZnj7o7VS5uG/XrjQ7Dd6dLAPM1eav7kczbK4eY/W69b1yzS9bp3WJ53Us21+W+bP92yauVJT2aRL\nS03bevvdwm2Pxcx733a+MlfzTxhb/gnznu779sWWuXO1Xrz42MerpETrWbOC3/pYj1mztC4t7b6v\nUlqvXNm/37yvhH/HWCz9HOnW5tAx6rb9jtGaf1TBb/6PMfObn96Udu6FH8NlY19YuTL9mgufM13X\nCea4lZebfft6LgnpANt1DlqdUw9da73Pe35XKbUJmAGE+hJsBB4H/mmQ7jNpOE73nk+YTDdJaSmk\nuooGJEHDtInj01Z0+L3kY/lp/R5gImHcIU89BanMyuzNNosX9y8TY2be8p6IjXd5qrwevhqUJ/vW\nuPQyXtkIu3cy/bkA997a3UWilPkd/v7v+94Drqkxj9ra9GOSZos3Yvj1r81v7LsmeuKNN4I5lDDD\nOYwPz0X4xGLmER69ZZZam7EnwXOPBgU62JAwBcGnPeCVHkzBxKfgjG1c9GaCT1darF4d/LnS0sJ0\nVVRUBMcxfFwyj2NHRzDaFD/68NCroCuljgdiWuvD3uu5wD8rpc7WWv/e2+1y4LdD1UjbNid3u1ep\nLexGgHTfd1kZ3H03NL5kszkkvFXT7ax/O9OnDd19g5ZlBH3bNmg/YMFDCSqmO3xqlM0Pftb/yNP1\n63v+LBaDmTPhgxNddk2Zw58zypM5zQ65lPHqqV6jZUHT3TVc9yPYe0ID4w9X8cauGpJe/2r8+P5f\nfJnzHrfeanzjxx1nJmD9KNmvfc38n0OHSBOyMC0twe9RUtK3SeTBwrbNOeafc0qZAiETJ5qbsj8v\nwJkOOjRh/1rSIXyMRr1v0fa4hX6x2hQjn/hU1/La/eUOq1ZZnHVW7hPf+Vgd4y+79QulZIp6mHjc\nHDPxow8jvXXhgYkYN8uLwMvAHd72BuAl4DfAL/DcMsd69NflonXgJli40Azpw0O9xYuDz8LDunWP\nN+m5K1bqdY8PzlivqckMNwdr6Dh/fnY3w6RJwf/4xPUhF88/YYbsd4zW828d3PGr71YYrOFxtt+q\np/+xcuWx3S9KGVdOf9q0eLHWlZXGdTMQm9atM66isCsq055rbvfcKsvimjtGm/cZrsHS0gwXjLdv\n2B14LPzfdd26wT1euRJ2t8Tj5hzO5n6B4HoczGtmpEKOLpdedxjMx0AEPczKlel+4eHwM4ZvDoN1\nkjY1hS7wjAvBZ/HdoQv/u2Wazy/U8QlNg3aBhG0ZjouvJ6E/ls89l/mC8N/yb+7XXKO7CepAbOup\n7eFtmZ2I8E07HjdzNl03r9Ob9EmXr8xJzNc93qQnXLuyy98ejwd/Jxbr/w2vr2S7Kfu/eVjYy8tF\nxAeTohb0we5N9sa6x0MTrHeM1qUTm3QsZnps/Z0M9fEvBn/iL9uFsPjuJl151Uo99QtN3UYhA/2/\n/qRVSYnW55xjBCgfF2JTU9DbC0869uU37u3GkI9JxsxztT8968zzj8qmrJOv8fjAz8dcbcp24w/f\nTEXMB5dcBT2yBS6G0394/m21PHd8ULOSp1fAr0zNynjc+IIH6tMdTnvq6kzQVGdndv9nSUmwfC4f\nPtqeJnF7o7YW7rijZ59uLAa3/cjFaXY4td1m8VXDl3kzW7nDXG3LrJnKC1+DqRu6Jl/9pbPQv9q3\nQuEjFYsGCdeFWde4dF4d1KyMPZggtTd9/fqoUYUdvejjujBrVsbqjCzMnw9PPhmtbIqua0TSnzzP\nZNbVLlvPDFahlD6S4JmHhkfUB0JmzVR2XQfT78vawVDKrJwaPz76qQSEAKlYNEg4DqTe8pacPb2C\n2IMJbvuyRcz/5Spd9IW1HB3rDnmmxMHAcXoXc4CmJjhyJMhyGQXbLMuscIrHu38Wj8Mb2gmWssba\n6TjNiYRdNfMs1l2YYNI7Kzjn+QQzyqtNzzwZzxrt/MADsGxZ9qykQnFT8NkW841te8uzWiwzrI3B\nq696S9hC6451soxDJ6RH/RUiXfYcY+03wLvvBq9Tp7o81OxQ0Vj4BZj9rJZhYjHjhmnZZqflsind\nZ0dmKV2Qr8iI9IVfTqC9aGlaTAfjf/5Ps6TzvvvMjTjfKRKE4Ud66L1gWfDjH5t18H5+k/1+TskJ\nTlqPb9chJ48tzQ3LgttuS98WyzwLKl2YWWuevZvWK+OWseDZOQVVmzUb/hr4eNw8z59vtqdSBME9\nzgpmvRkNd0s2HAd42zJuFs93XlMDmzaZuRx/hKK16a1LL33kID30HKipMQmzumqV7sYk5MqIGu0p\neKnQWLWKbgEsu3ebBFuvHHZp+UzgZ2bXdaGI2yC7ZaGSGShWX58xGmmxiO23+Oz3ottz9W9aflIv\nP6oZjE033GAKwWht3GsSoTlyEEHPkXDEpf+8fr1F+4sJyv4vhxvnFL47Iowfpu9jWeamdeESJ03A\njzsOPozYTSt8rPyCJWHKy6MdtWhZJs1uT8VYqqtNOcPMXPxC8SOrXIQuamvhO/d48wLeior5HySY\nN49uebyjgr/ypaOj/zlqokixFM0QDLJsUegzXcv+TnZhgkO8xWbbw9H0M4cRcROijgi60C9ct+eh\nvCAI+SFXQRcfupBGT9kZBUEofGTZoiD0Adc1cw2yFFAoRKSHLuSFKPm1fTfUwYPQ2Ni9UHW24s+C\nkA9E0IUBExbn3Ydc1id6Tn7li6NfSDkWM4Fb/an4NBz0lB/Gr8Cz+5CXZyXezpZny9i6NcGsWdFd\nFSREGxF0YUCES/ip8V4Ss+PbYVQZv7g6kbZKJps4plIm8+PkyYXZU89W/lCpYH33ogbH2Out29+4\nq56N8Q1dAg8JEXVh2BAfujAgHMcIdDIJnZVOWiqEZKWTFthTX589E2IyWbjJv/zyh346hNh4l8sv\nh+uug82bYfvP7fREWZD2GzTscPLXeGHEIT10YUCE64fqt2xSGQW0+R/BvgcPpn9XKfMo5MhNy4JF\nP3RZfdC4VVLJMn75SILkY5aXUsDLDzPBCbIeTt0QqchaoXgQQRcGRDh3yqFDFnc+GGQBLH3H6sox\nAiYTYJjLL4cZMwp/YnTXISctHULHaQ68EWqwn4kTUxzkH8Yl2HVIfOjC8COCLgwYX4znzAHVZhHb\nZ3HZZd2r1ldXw/33G590aWnvVe0LharptvGHh9LuJr0UxEqZlAL/8A8wZox/c7Io9DTKQnEigi4M\nCr4vPZUyAjdjRnextiyzX1SWK/qYXnaia+XK5M9b/S6TJwhDiYT+C4NCeLVLVErWCUJUkNB/YVjJ\nzEMuYi4Iw48IehFQKFGXA80DI4nBBGFgiKBHHNeF2bNN9Zp4HO65p3CjLo9FZtDR+vXwzDMi6oLQ\nFySwKOLU1welyJKnuCx4sJavLI1W5qi6OvjSl9KDjjo6gt66JMQShNyQHnrE6NG94hVzJt7OxmQZ\nv69K8N8Nhd+9rauDBQt6/lwmWwUhd6SHHiF8cfvOd+CCC+DSS2HaNO/DCU5ayPlz7zosWZLP1uZG\nQ0PPn02blp5awE+IJQhCdkTQI4TjwJEjwfstW+DOO703zXZ6TpFmm4cfzkMj+8jUqdm3KxWs8S4r\nM/MDUvBYEI6NuFwiRDYxe+0170VLRk6RFouJs4avbf1lzJjs27WGQ4dkOaQg9AUR9AhhWTB3rumZ\nZyWUUyQehx/8YPja1l/8bIaZKWoBdu0yz1IWTxByIyeXi1KqWSm1Wym1Sym13dt2p1Lqt0qp3yil\nNimleuhrCYPJk0/CNdfASScZcR892oh3SYkJt1+8GFauhG3boiGClgVr15r2Z1JVNfztEYQok1Po\nv1KqGThPa/1+aNtc4L+01p1KqVUAWutjTsNJ6P/gUyhBRQPFt+PQIdMzr6qK5np6QRgKhjz0X2sd\nHvj/Gvhif/+W0H+KxR1RLHYIQj7JdZWLBp5SSu1QSmXrN90ANGb7olKqRim1XSm1/b333utvOwVB\nEIReyFXQZ2qtpwLzgJuVUl3rJ5RSdwCdwMZsX9Ra12mtz9Nanzd27NgBN1gQBEHITk4uF631Pu/5\nXaXUJmAGsFUp9VXgC8AcPZx5eIuEYvF/+9Q1ukHO8DFW0dhWbMfJp1jtGsn0KuhKqeOBmNb6sPd6\nLvDPSqnPAouBi7XWHw5xO4uKukaX9QmHnZtsUm9ZXSHtEN1sg0vWBnU3tzxbhnowAW8b255+GnYf\nCsQ+SmXZMlMPrFkDO3eaz/zyelEUxWx2tbZK0Y6ok0sP/ePAJqWUv//DWusnlFKvAeXAf3qf/Vpr\nvXDIWlok1DW6LHh2DhzfDleVwYYE7Qcs6uvhpz+Fzk6z3wMPGCGMwkXlunDXzx2w2yCWAtpMXdG9\nFm1tcNNKl11TArGHRGREPZx6oK0NbrrJvAZzvJQyx8y/cUXheIGxq63NVJg6csTYpXVQVq+kxCwn\nlZVG0aJXH7rW+g2t9RTv8dda6+972z+htT5daz3Ve4iY50DDDict5woTHLSGV14JxByilbfEcSD1\n5wpQKTN9rlJQfqjr811/ctJsbtjh5Kmlfce2jbgpZR6pVPBZZ6cJiNLaiKM/uooCFRXptiSTwXut\njV033ywZLqOG5HIZZqqm291yrqRSJhAoTCwWnbwltg3q+FbQCpS38YJ/MRkgAd6002yumm7np6H9\nYPfu4EarFMTGuzCzFipdYhG+elpbe/igMrAvmYxOp0IwSOj/MOMXHP72jx0+eNHuCtXX2oi4X2T5\nnnuiM3y3LPjbT9ps0SWgOzxRT5m8Ml46gnFPJhh/scONc6LjQ3dd00v1XSydp7jwFeM6IlnGJQcT\nJP7NIpUy6Qt8n3oU8JOehXPQh1MwkyxDb0hQURGNYyUYIiPo/ox81CdtXBdad1lM3GexqyX9M6VM\n+P6NN0bLd+m68HS9BZPXwudvBlKQKjdJwjwObrf4w28sJkconN9x0t0SaSmKdTv/+aqDTlrE43D3\n3dE6Hy3L2Ld6NTz2mOlQMNFBh+xjgkNDg8XkydGybSQTCUH3Z+T9SZxYDMrLo1fsINOOTJJJeP55\nM8yP0kXkOF5yrRdq4N3JaRkfw3R0mH2jYpdtm/Osrc2cc+NOqqBFK0jGIFWGftMGzHHzV75ECcsy\nuX8aG01PPbbXJpksM2LuuQOf2m/cgVG71kYqkRB0f6WBL4KpVDBpGKWTLNOObGgdPdtsO3AXhTM+\nzphhhM7PpFhaGp15AUhP3XvoBJfVBxd5E79xeGJNtxtWFHEcM0egNfC2xfwPEuwvczi82+Z3+63I\nXmsjlUgIuu/vC/fQo1jsINMOn9JS+Pzn4fHHTW8varZZFvzkJ7BwoScMBGubIbpr6yHIMXPp9xzP\n3ZKCpGLS9FZefylYxx0l/3mYsC+9rAwWX2VhWZYZTW4LtkfpfBzJ5JRtcbAYSLbFYvKh+3aEA1Qs\nK/qRe64bbfE+Fl3xAzHjjlh3YaJoomF7Ou+ifj4WE7lmW4yMoAtCvgmnNojKSh2hOBBBFwRBKBJy\nFfQIh0YIgiAIYUTQBUEQigQRdEEQhCJBBF0QBKFIEEEXBEEoEkTQBUEQigQRdEEQhCJBBF0QBKFI\nEEEXBEEoEkTQhay4LtTWSgkyQYgSkci2WKgUa/KizIrwkgtbEKKBCHo/KWbRC1e6P3LEVLXZtKl4\nb2CCUCyIoPeTsOgVWwGAioogrznA5s1w8cXQ1GTsLS0tLnsFoVgQQe8DS9a6PPqCw0mjKth/qBVO\nt4nttdCnuTzU7FDR2D2tambK1ULv5bouLFrUvarS1q3B6/Z2k/e8ENsvCCMZEfQcWbLWZfXBOXB6\nmylDNjYGE8tNKbLPLuKVeDsLni0DEl2i3lUUId7OlmfL2Lo1wc9+aJFM5r8mak/FKBzHVFTqotLt\nsUaoIAiFhaxyyZGHm5ygBBmY51g7nNMQVIKPtdOww+n6TsMOJ+2zh5916Ogwvd+2NiOe+cB1zQjh\n3nvNY/bsYDVLRUWod17pwnVz4DPfhesvhnPrAHMzimrJNUEoZkTQc2SisiFZZiq+A6RM5XdeqfK2\nxyFVRtV0u+s7VdPttM90c/BZPJ6/Oo2OExRuhmAOADKq109wIN7m3bw64PM3M/9Wl6efFneLIBQi\n4nLJkR/cYnHRVQmSpztwpIITx7XywW9seNuC9yZTOdNh2bXpPnTzOkHDDoepY2zuft+iLWaKXK9d\nmz9RtG0zsdnebt77RYBdF37609COzbapcK9ToABSfHiygyVqLggFiZSg6wPhCU3o+7LFQpoQzeZD\nr62FO+5IX+Fy9v9dx+8/eTOQglQ56y5MSD1NQRhmpKboMFBIAj0Y+L51v+deXg5PPw27D0lxZEHI\nJyLoQr/oafWLIAj5I1dBFx+6kIZliYgLQlTJaZWLUqpZKbVbKbVLKbXd2/Z3SqmXlVIppVSvdw5B\nEARhaOnLssXZWuupoW7/S8CVwNZjfEcoAIo1c2JdHVx6qXkuJor1eAlDT79dLlrrPQBKqcFrjTDo\nhJOIqfEu517pcOOc6E9ufmWpy8ZfOdBss2WLxeuvw6pV+W5V/6mrg4YGGDsWfvYzCiKaWIgeuQq6\nBp5SSiWBdVrrnPtESqkaoAZg/PjxfW/hMOFfUFVVUFOT79b0n8yVN11JxE5x4eo5PBdv57mMFAVR\no67RZWN8Dnym3QRubUhw110W8+dHU/zq6mDBgtAGL93CkWab1astZswonpVUwtCSq6DP1FrvU0qd\nDPynUuq3WuucXC2e+NeBWeXSz3YOCf6KjldeCZJPbdlCZHt72VL62rZ5feRMJ0hDoNu58+cOk8dY\nkRIJ/2b1ULMD4wJbmOCg91mRzQDZ0BB646dbiJub1eYNCR57zJLeupATOfnQtdb7vOd3gU3AjKFs\n1HDguiaHyb33pmcSBLjrrmj6L/3EWslkkCvGsowQzJ9qp6UheD1hM2dOdOz0b1bf/S688ridZgsf\nVsDMWg6dEBFjMqiqCr2Z4KTl/2GCQyqVnp5BEHqi1x66Uup4IKa1Puy9ngv885C3bAjwe3gvf+Cy\n6Y162i4BXqzulkVQ62jm+w4n1kqlzHvf5sVXWcw7lODOnzu8nrDRey3a44VtZ9h95LuOUinM8dqQ\nMOL3YQXMW4SOt7P6YBlnNRa2KylbMFpNjRkV3nUXpJptc7PS7eZm1WwTiwXpGQThWOTicvk4sMmb\n/CwBHtZaP6GUugK4GxgL/FIptUtrfenQNXVg+D28Ix9z4TobzvHCIac9AP/2dJqojxoVzYuntdXk\niUmlzPPOnXBTrUuy0iFeZ7PtYYv6r1nM+Xdojxe2SPgjKN999MVvuaQucOBN2xwr/zGzNs2V1LDD\nKVhBP1aVK9/Ft3p1cLM66bDN319tMWaM+NCF3OhV0LXWbwBTsmzfhHG/RAK/h8cZjskc6C/O8Ya1\nvqCXlMCaNYV98fSUcsC2zcoIXzB+3eKSvMb4Y5PJMm5fm+CZjRaJROGnLKivD/Kyt431JkFnt8Ms\nMwnadQPO6NGGs10WGr1Vudq1y3vh3az+ANz9qvjOhdwZMelz/clB3rIhVWrW7Wi6hrU+WpueLhTe\nOmfXha9/3fRcly2jmw/c95evWGGe3z/BSfPHvqGdrv2WLo2QSGTxK8fj5uYVP2BR8nCCGR+uKPjE\nYf45GO9hdDR1avfviO9c6AsjIvTf79GuWQOtrRY/a3LYpb2EJSEfulLmQquogCuuMLU0wax8gfwu\nZ/SH60ePBtkQwxOf4V770qXm86t32Kw+GPRer77Azk/j+0F1Ndx/v8nbHnvbJhnqhcf22txzD0ye\n7NtsRSKlr3/D7Wl0NGaMOQf94+ufj4XqFhMKj6JPzpXNbwmm6LFf5MGPjSopgW99C+6+21S7D5df\nm3SixcsvD37bcnF9LFnr8tOnHP7wgm3yr4dYt84IW9jf/K//akYZtg2bd5g6qFeea7PqlsIXvTDh\n32f3IZf1CYdT220WXxWt5Za5Ej5X43G44QZJkCYYiirbYl/T1IaDhFpbjXsimTQXyYoVpgeeFsjh\nEY/D9Onw/POgT0tfD8yGBPPPs1i8eHAusGNNkIXpqmUaagctVtfN5pzjbMYetbotvVTKTO4Op/9V\nMjUOnGJLySwMDkWTbTFX4fMJR91t2QInXVJHcsEaKDlC8p2pvPzBYt5zev4DO3Z4Q96w39YLXtm8\n2eLJJwdHJMMTZEePGiHM9jcffcGB09PbEYtB6loj8nuSZezZkAC6L70Mu2SGGn9Vij+RuX493Hij\nCHtfkWyXwkAo+EnRzGCZ1Y+4XPq9WuoasweRpEXdnVvHHy5cAGP3wF81w6c2s7H0Yg59pPt3/eV+\nyaS3wV894QeveBOngzVJZdtmRABGfO+7L/vk65XnZrTjLZtTL3S6TRJmQ2szGhkOulYReXR0wL3/\n4TLz9lqWrI1mwI8gRI2CF/S0YJmpdWweM4stHd9lQZPNyTd8vZtYpEXdTfLUXYUesQ6ef89h8WKY\nNAnOOQfmz+8+GTW61WLxuAST3lnBGc8kKH3H6nF1Qq6Es+hZlvGR+v77ZBJuusmsYgmvXFl1i8U1\nyQTKWQH1CUa/b5nJzSw3m0y0hkWLhicatKIisAXoCmFP2ctYfXBOjzdgQRAGj4L3oXfVuTzNhesv\nDtaQa0ArSI5i8bhE2oSf70P/TWkdB8/LcJanSuHfnmHl162u1SBht05JCVx/fXdXwUB9m7lMzkLP\nvu/M/1/XGEwSzvsfFjt3mv1OPBEefdREHmodzBv4tg4FXUFbR0IbZ9bCZ5aZUUQqxqgDlzDjL8uZ\ndKIlbhhB6CNF40OvqAj5tFUyEHOAmAbdzqMvOKwK+ZBravwlhjVMuhb2nLgGSo/AwanQtJjR71tp\nvezelpP5+wxEhLIFlSxdCmvXws03Q2en2c/3fS9fbh7+/8z8/zXzrB7XXM+fn37zGOplb5nuFiBw\nWdEGKsXRU55ia3IbWzckeOABi6efFlEXhMGm4AW9K5y92YZkOdBmeubEMJXoy4yfuQdeebAre6/p\n5U7JLtpDPRnlB5WERdZ1jX0//rEJ07//fiPsqRQ89RRs29a/CdhcblCDSdi2WMwbcfj5VuzlMPEp\niKW6JnXbI5wZURAKmYJ3uYRdFanTXPQZDjTbxGNwygUOV18QnfXVYbcJZHfBLF9uxDyVGh53yWAR\ntm3zZli92vvATwcb85JNbUhQ/p700AWhL+TqckFrPWyP6dOn6/7Q1KT1woVal5VpHYtpXVKi9bp1\n/fpTBcPKlVrH41qDeV650mxvatJ69GizbfRo8z6KrFun9YwZno2VTZqLVuqpX2jSCxdG1yZByBfA\ndp2Dxha8ywWCyjvJpOm5KhXkW8kXA50kzeaCgeF3lwwVNTXmGO3YAbRYxA9YfCkiow1BiCqREHTo\nWQDzgesGq1NKS+GZZ8x2xzGTuH7Y/bHE+FjCXSzBJf09ZhItKQj9o+B96GH6cqH7+x46ZNKSTp3K\ngPNK+3/ziSdCVY4qXUad49Dxe5vUWxZaD1/YfRSEr7c2+p9XVEDjSy6vtjn8/imbZLOFUnDZZQxa\nugVBiCpFlculN+oaXRp2OFRNN9Xss2Um9CkpMatK+po5MW1yNuX93Yz6j2l5uoEZM0z4ey499p7s\nWp9wGHXA5qQPLcaNg2nTzN+rqDBBQ7mmRChE/N+0rQ1Sp/b8W5aXI5OowoimaNahZyMzC9+CZ40Q\nbPGq2e98zMoq5mCWBd5yi8lQ2BeBCK8j7yJLvpewoD/3nHnEYvS5yG9XUq7j2+HMdIFTyqyASaVI\nqzcZJcFzXbOip63NiwQ+xm8ZRfsEIR8UfOh/JuFiwRdfDD/a7KTlNal9xGHdOi+ydGat6UVn0NER\nWlaXI74/OC28vYd8L5n0tciv68JdPw/b1WbWc3u2aB0EIqnxLqkLa3niJTdyBZ/95ZmAqQ2qFSRj\n3X7LWExyggtCLkRO0P1kXalTXTrOr+WVHRVpotr8tB2kvv3MMvNc6TJunBEGn82b+1aJyJ/EXLDA\nTIQCQfDM0yu6uVvC9LXIr+OA7rpZxExQzsSnumzxSZ3qoq+dg7aXsfXMOVz8lcIX9bo6k1bh6NGQ\nmFe6MG8RqBQQhyfWED9gfst4HO65R3rngpALkXO52DZGAK4N/K0nPbeGT/xNK4d32+zJUjiYCQ7V\nF5joxOeeC/5WQ0PffOn+6pPq6iDv94knWvzsZxbNLd33Vwouv9z40vviQ7dtKF1u0d5DpGXXjSPD\nTdFxmoPjFG7xhyVLehgZddmRgqQi9pFW7rmn/3MPgjBSiZygWxZ86rMOr4SE7A9HWnlx7VK++U3Y\ns4VuhYNpthnzOTNBGRb0tMyMfWxDWGRWrTJi9eijcP755tmfrOzvCg2lMMLtLIcztqXZ0kWGnaX7\n7IJ1Tbgu3Hln+raTTzarkDr22uiQHbd90c5ruT9BiCqRE3SAb863WfBsumC3t5tliYsXw113WaQ2\nJEzP7y2b0a1WWk/Pr2Y0mKKxapV5wMCXEzpO4COnxYL6BOpMh9heGw5YlJTDvHkwbpzFieMSOM1e\nabaHCrd37jjdJ6kvuMAcL8exOHRCgl2HgpVKgiD0ncguW1yy1uXOnzvoN2xosdKWtoXXNkdx2N5T\nbR6ZNdkAAAc7SURBVEko/HXnPeG6pt1+VsZ43CQfi5odgpAPRsQ69GKuYRmFoKG+UszHSxCGkhEh\n6ELxU4w3NkHoK0UdWCSMDOrqTBBYMhkEZu3ePTRzIIJQDIigjwCi2Mt1Xfj6911S5zvQbHN0n8Xt\ntwc5dLZsMc8i6oIQIIJeBBxLsLPVMo2CqK9+xCUVijXQGxJs25be8B/9SARdEMKIoEec3gS7vj5I\nUpatVmmhsr/M6RYcpjMicffsga9/PUhYFqURiCAMBSLoESdb8WnLMpkaf7TZYc8TNloblRtordLh\n5MY5Ns+FYg3UXpvSMpOHx5/H1xruvde8Vsrc0CQrozCSEUGPONmKSHRlahzXDtemZ2qMSnZGE1yU\noGGHw9QxNmMWWlRUwMaNoVz0IfwRSH19YdslCENJToKulGoGDgNJoFNrfZ5S6iTg/wUmAM3Al7TW\nfxyaZgo9kVn5CLxMjXb2VLR+T7ZQUwSEqZlndUWN+q6lI0fy3ChBKGD6km1xttZ6amgt5O1AQmt9\nNpDw3gt5wLJMrU6/9mqQqTFI6+sL+YIFhe9uyYbvWsqGUsGNyo+oFYSRyEBcLpcDtvd6A+AASwbY\nHmGA2DaMWmFx9EGTy+aiSptrllmRnzT0XUtdBTE8/ApUUbdPEAaDnCJFlVJvAn/CuFzWaa3rlFKH\ntNZjvM8V8Ef/fcZ3a4AagPHjx09/6623BrP9QhaiuO48F8I5enbuNNskhYAwEhjU0H+l1Gla631K\nqZOB/wRuBR4LC7hS6o9a67861t+R0H9BEIS+k6ug5+RD11rv857fBTYBM4B3lFKneP/sFODd/jdX\nEARBGCi9CrpS6nil1An+a2Au8BLwGHCdt9t1wP8ZqkYKgiAIvZPLpOjHgU3GTU4J8LDW+gml1PPA\nz5RSNwJvAV8aumYKgiAIvdGroGut3wCmZNneCswZikYJgiAIfacv69AFQRCEAkYEXRAEoUgY1opF\nSqn3MP72/vAx4P1BbE6+KSZ7iskWEHsKmWKyBXK35wyt9djedhpWQR8ISqntuazDjArFZE8x2QJi\nTyFTTLbA4NsjLhdBEIQiQQRdEAShSIiSoNfluwGDTDHZU0y2gNhTyBSTLTDI9kTGhy4IgiAcmyj1\n0AVBEIRjUBCCrpQ6XSn1tFLqFaXUy0qpb3rbpyqlfq2U2qWU2q6UmhH6zlKl1GtKqd8ppS7NX+u7\n01d7lFIV3v5/VkqtzW/ru9MPe/5WKbVDKbXbe/5Mfi0I6IctM7xtu5RSLyqlrsivBen059rxPh/v\nnW+35afl3enHsZmglDoSOj735teCdPqpa3+jlHK9/XcrpUb16Z9qrfP+AE4BzvVenwC8CkwCtgDz\nvO2fAxzv9STgRaAcOBN4HYjn244B2HM8MBNYCKzNd/sHwZ5pwKne6/8B7Mu3DQOw5TigJPTdd/33\nhfDoqz2h7/1/wM+B2/JtwwCOzQTgpXy3exDtKQF+A0zx3lf0VdcKoki01voAcMB7fVgptQc4DdDA\nid5uHwX2e68vB/5da90GvKmUeg2T0tcd1ob3QF/t0Vr/BfiVUuoTeWhur/TDnp2hr78MjFZKlXvH\nK6/0w5YPQ18f5e1XMPTj2kEpNR94E/jL8Lb22PTHlkKmH/bMBX6jtX7R+05rf/5pQT0wd929nsHn\neK/fBvZhoqUA1gJfCX1nPfDFfLe9v/aE9v0qBdhD76893v5fBJ7Kd7sHYgtwPubG9Gfginy3eyD2\nAB/BdHw+AiyngHro/bBlAuamtAt4Brgo3+0eoD2LgAeBJ4EXgMV9/j/5NjTD6I8AO4Arvff/ClR5\nr7/kC0NUBD1Xe0L7F7Sg98Oev8a4w87Kd9sHaou3/RzgOWBUvtvfX3uAu4Avea8LUtD7YEs5UOG9\nnu4J5In5bv8A7LkNM3L6GMbV5wJz+vS/8m1syOhS7870D6FtfyJYWqmAD7zXS4Glof2eBKx829Bf\ne0KfF6yg99UeoBLjM7ww320fjGMT2u+/gPPybUN/7QG2Ac3e4xDwB+CWfNswSMfGifix+TKwIbTf\nMuDbffl/hbLKRWF62Xu01v8S+mg/cLH3+jPA773XjwFfVkqVK6XOBM7G9JwKgn7YU9D01R6l1Bjg\nl8DtWutnh7OtvdEPW85USpV4r88APoURw4Kgr/ZorS/SWk/QWk8A1gArtdYFsbKqH8dmrFIq7r2e\niNGBN4avxcemHzrwJDBZKXWcd85dDLzSp3+a7zuYdyeaiZko+A3GH7YLM/s7EzNUeRH4b2B66Dt3\nYIbzv8ObMS6URz/tacb0lv4MtACT8m1Hf+0Bvkvg2/QfJ+fbjn7aci3Gf74L49ecn28bBnquhb67\nnAJyufTj2FRlHJvL8m3DQI8N8BXPppeA1X39nxIpKgiCUCQUhMtFEARBGDgi6IIgCEWCCLogCEKR\nIIIuCIJQJIigC4IgFAki6IIgCEWCCLogCEKRIIIuCIJQJPz/XLDVqr4ycdUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12eb89d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_gfa_targets(targets, rfluxlim = 1000, rfluxmax = 10000, tiles = None, buffer_arcsec = 100, tileid = None):\n",
    "    \"\"\"\n",
    "    This function takes a table of targets, as well as optional parameters \n",
    "    including a minimum flux in the r-band, a list of tiles, a buffer in arcseconds, \n",
    "    and a tileid number and uses blue to plot the targets on the GFA satisfying a minimum flux_r\n",
    "    on either a specific tileid in the targets table, or all the targets on any tileid returned \n",
    "    by get_gfa_targets. It also uses red to plot all targets on the GFA above a maximum flux_r.\n",
    "    Parameters:\n",
    "    targets: table with columns RA, DEC, FLUX_R\n",
    "    Options:\n",
    "    rfluxlim: (float) r-band flux limit; default 1000 = rmag 15\n",
    "    rfluxmax: (float) r-band flux maximum; default 5000\n",
    "    tiles: table of tiles, default to desimodel.io.load_tiles()\n",
    "    buffer_arcsec: (float) additional buffer region around GFA to include\n",
    "    tileid: (int) a unique indentifier for a tile pointing\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    if(tiles is None):\n",
    "        import desimodel.io\n",
    "        tiles = desimodel.io.load_tiles()\n",
    "    valid = get_gfa_targets(targets, rfluxlim, tiles, buffer_arcsec)\n",
    "    if(tileid is None):\n",
    "        plt.plot(valid['RA'], valid['DEC'], 'b.')\n",
    "        brightindices = np.where(valid['FLUX_R'] > rfluxmax)\n",
    "        plt.plot(valid[brightindices[0]]['RA'], valid[brightindices[0]]['DEC'], 'r.')\n",
    "    else:    \n",
    "        indices = np.where(valid['TILEID'] == tileid)\n",
    "        plt.plot(valid[indices]['RA'], valid[indices]['DEC'], 'b.')\n",
    "        brightindices = np.where(valid[indices]['FLUX_R'] > rfluxmax)\n",
    "        plt.plot(valid[indices][brightindices[0]]['RA'], valid[indices][brightindices[0]]['DEC'], 'r.')\n",
    "    plt.show()\n",
    "    \n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "sweepname = 'sweep-280p050-290p055.fits'\n",
    "sweep = Table.read(sweepname)\n",
    "\n",
    "plot_gfa_targets(sweep, tileid = None)\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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}