Tolerance On GFA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 392,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#x = [-125.10482863, -129.83038525, -159.04283509, -154.31646944]\n",
    "#y = [-370.01790486, -384.56223777, -375.05643893, -360.51151824]\n",
    "#plt.plot(x, y, 'b.')\n",
    "#plt.show()\n",
    "def shift_gfa_points(deltax, deltay):\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",
    "    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",
    "    #print(x_all, y_all)\n",
    "    return x_all, y_all\n",
    "    \n",
    "#x, y = shift_gfa_points(50, 50)\n",
    "#plt.plot(x, y, 'g.')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def degrees2xytolerance(degrees_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 = degrees_tolerance / (10**3) * 3600 * r_ps\n",
    "    y_tolerance = degrees_tolerance / (10**3) * 3600 * az_ps\n",
    "    return x_tolerance, y_tolerance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def is_on_gfa_test(targetx, targety, degrees_tolerance):\n",
    "    import desimodel.footprint\n",
    "    #import desimodel.io\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",
    "    x_tolerance, y_tolerance = degrees2xytolerance(degrees_tolerance)\n",
    "    \n",
    "    #AREA_BOUNDARY = 469.8037221 + MIN_TOLERANCE\n",
    "    targetindices = []\n",
    "    #gfatable = desimodel.io.load_gfa()\n",
    "    \n",
    "    # x and y hold the 40 new GFA coordinates\n",
    "    x, y = shift_gfa_points(x_tolerance, y_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",
    "        #gfaarea = a1 * a2\n",
    "        #gfaarea2 = a2 * a3\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",
    "        #use np.where!!!!!!!!!\n",
    "        # This if statement may be redundant and thus unnecessary\n",
    "        if(any(targetarea < AREA_BOUNDARY) and all(targetarea > THRESHOLD_AREA)):\n",
    "            print(\"The target is on GFA\", int(gfaid / 4))  \n",
    "            newtargetindices = np.where(targetarea < AREA_BOUNDARY)\n",
    "            targetindices.append(newtargetindices[0])\n",
    "    \"\"\"targetindices is an array 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 targetindices array.\"\"\"\n",
    "    #print(gfatable['X'], gfatable['Y'])\n",
    "    return targetindices\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 395,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The target is on GFA 0\n",
      "The target is on GFA 0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHdNJREFUeJztnX2sZVdZxn8vM7QWjEI7V+jH6ECpGvqFzA22KtoIDqV+\nDFpNqDXxAxkiGTVSrbaN9TZKo8FIKB9lJiAQSjCxiB+YCo6xapxqva3tdKAFpzqUjkBvQY2g2FJf\n/zgbezq959y9z917r7Xe9TzJypyz9ln72c9aZ5939vOcs6+5O4IgCELdeErqAxAEQRDSQ8VAEARB\nUDEQBEEQVAwEQRAEVAwEQRAEVAwEQRAEVAwEQRAEVAwEQRAEVAwEQRAEYGvqA2iLbdu2+Y4dO1If\nhiAIQlG44447Hnb3pY1eV0wx2LFjB6urq6kPQxAEoSiY2SfbvE42kSAIgqBiIAiCIKgYCIIgCKgY\nCIIgCKgYCIIgCKgYCIIgCMCWlZWV1MfQCvv371/Zs2dP53G33QY33QRbt8L27Yv397kvcY/Pcc/+\n27j36pv4z//eyrN2PpF81rah+6NzR9eX+pxpi+uuu+7TKysr+zd8obsX0Xbu3OldcfCg+0knuW/Z\nMvn34MHF+vvcl7jH5zi076B/kZP8Ubb4FznJD+17nHzWtqH7o3NH15f6nOkCYNVbfMaGtoluvRUe\neQQee2zy7623Ltbf577EPT7H5z5wKyfwCFt5jKfyCJ/7wOPks7YN3R+dO7q+1OfMEAhdDC66CE44\nAbZsmfx70UWL9fe5L3GPz3HKpRfxCCfwKFt4lBM45dLHyWdtG7o/Ond0fanPmUHQ5vIhh7aITeTu\nvm+f+65dk38309/nvsTdnuODVx70dz7vev/glQc31X/r5fv8H07e5bde/mTyWduG7o/OHV1f6nOm\nLWhpEyX/kG/blBnE5J63LYK3XCt3dH2pz5kuaFsMQttEJfnatXLP2xbBW66VO7q+1OfMEAhdDEry\ntWvlnrctgrdcK3d0fanPmUHQ5vIhh6bMoBzurn7+vG0RvOVauaPrU2ZQUDEoyTuPwp2jvytuZQbK\nDJQZFOOdR+HO0d8VtzIDZQYbI3QxKMk7j8Kdo78rbmUGygxaoM3lQw5NmcGwHIv4+X35/IuMKclb\nrpU7ur5554wyg8yKQa2+fV8+/7xtJfm74ta6jsmtzCBD1Orb9+Xzz9tWkr8rbq3rmNxVZgZm9utm\ndsjM7jKzj5jZaVPbzjOz28zso2Z2j5l9VdO/s3l+xMxuMDPbrIhZqNW3L+nePWNwiFvrOiZ3lZkB\n8DVTj38OeHvzeCtwCDi/eX4KsKV5fDtwAWDALcDL23ApM4h7754xOMStdR2Tu+rMALgKuLF5fAlw\n0zqvORW4b+r5ZcC+NvtXZhDz3j1jcIhb66rMYITMwMxeb2afAi4Hrm26vxFwM/uwmd1pZlc2/acD\nD04Nf7Dpm7XvPWa2amara2trnY+tJN9+Ee4I3usYHOLWuioz2BgbFgMzO2Bmh9dpuwHc/Rp33w68\nD9jbDNsKfAeTAvEdwA+a2Uu6Hpy773f3ZXdfXlpa6jq8KN++1nv3jMEhbq2rMoMWaHP50KYBXw8c\nbh6/EnjP1LZfBX6JkW0i9zwzA927Jy9/V9zlcuTKXV1mAJw19fhngZubx88E7gSexuQq4QDwvb5+\ngHxJG64omUFp/mcEDnFrXZUZDJ8Z/GZjGR0CdgE/31xt/BvwO8A/AHcBd7r7nzZjXgu8AzgC3N8U\nhEGQY2ZQmv8ZgUPcWldlBhtjU8XA3S9193Pc/Tx3/353Pza17SZ3P7vZfuVU/2rTd6a7720q1yDI\nMTMozf+MwCFurasygxZoc/mQQ8shMxj6O/25+p8ROMStdR2Te9a5HzYzGLOlzgxq9T8jcIhb65oD\nd/TMIGv06dHV6n9G4BC31jUH7tCZQe7o06Or1f+MwCFurWsO3MoMCssM5n3Xv1b/MwKHuLWuOXAr\nMygkMyjNg4zAHV1frdzR9SkzKAx9/QZg3rZa/c+SOMStdc2BW5lBQpR0X/9auaPrq5U7uj5lBoXZ\nRO6TS6vrr3/yJdas/kP7Dvpf7rr+CZd9G20buj86d3R9tXJH17cId9fPo1n9XYBsIkEQBKE12lSM\nHJoC5Jjc0fXVyh1d3yLcCpATQgFy/tzR9dXKHV2fAuTCoAA5f+7o+mrljq5PAXJhNpG7fnRWAnd0\nfbVyR9e3CHdfN7DrAvSjM92orgTu6Ppq5Y6ur09uZQYjQDeqy587ur5auaPr65NbmcEI0I3q8ueO\nrq9W7uj6+uRWZpBpZjDPo8vRg4zAHV1frdzR9fXJrcxAmUE4/zNHDnFrXXPmVmYwApQZ5M8dXV+t\n3NH1KTMoDMoM8ueOrq9W7uj6lBkUZhO5KzMogTu6vly53/qqK/21u57nb33Vla36S9NXErcyA2UG\n4fzPHDnE/WTu9755n590Db7lWvyka/D3vnnf3P7S9JXErcxgBCgzyJ87ur5cuf/qzg/wyBZ47Cnw\nyFMmz+f1l6avJG5lBiNAmUH+3NH15cr9XS+8lBMegy2PwQn/O3k+r780fSVxKzMoMDOYdV+RHD3I\nkrij6+uTu6vPP6vf3f2Kn7ncv/UnT/YrfubyVv3R5zYld9d7FnUBygz6zQy6/g2EWv3PHDmicHf1\n+Wf1LzIm+tzmyK3MoEf0mRl0/RsItfqfOXJE4e7q88/qX2RM9LnNkVuZQY/oMzPo+jcQavU/c+SI\nwt3V55/Vv8iY6HObI7cyg4wzg1nbcvQgS+KOrm+eb9/V6+/q88/qX2RMjnMbnVuZgTKDcP5nrfr6\n9O2VGdTFXVRmYGa/bmaHzOwuM/uImZ3W9D/VzN5jZveY2b1mdtXUmJ1N/xEzu8HMbHPXNrOhzCB/\n7uj6+vTtlRnUxV1aZvAGdz/P3V8AfAi4tun/EeBEdz8X2Am8xsx2NNtuBF4NnNW0izd5DDOhzCB/\n7uj6+vTtlRnUxV1sZgBcBdzYPL4M+BNgK3AK8AngZOBU4L6pMZcB+9rsX5lB+f5nafr6+k5/n759\nysyg6z2L5m2r9T0VOjMAXg98CjgMLDV9TwV+D1gDvgjsafqXgQNTY18MfKgNjzKDsv3P0vTl6Nvn\nyF3aupbEnV1mYGYHzOzwOm13c2VxjbtvB94H7G2GvQh4DDgNeA5whZk9t+tVi5ntMbNVM1tdW1vr\nOlyZQQHcuerL0bfPkbu0dS2JO7vMwN1f6u7nrNP+6LiXvg/4iin5o8Cfufuj7v4Q8LdMrgqOAWdM\njTmj6ZvFvd/dl919eWlpqYsuQJlBCdy56svRt8+Ru7R1LYm7qMwAOGvq8c8CNzePfxl4V/P46cDH\ngPOa57cDFwAG3AJc0oYr58xg6HsW9bmvHLn73Nci9+gZ+jv9UTKDrr8/mLetpPdUSu5Zny1dwBiZ\nAfABJlnBISaB8elN/1cDvw98tCkEvzQ1ZrkZcz/wFsDacOWaGQydJfS5rxy5+9xXrb69MoO6uLui\nbTHY1FdL3f1Sn1hG57n797v7sab/C+7+I+5+trs/393fMDVmtRlzprvvbQ52EIyRGQydJfS5rxy5\n+9xXrb69MoO6uIeC7k20ycxg6Cyhz33lyN3nvmr17ZUZ1MU9GNpcPuTQcs4Mhv79QZ/7ypF7kTF9\n+fyLjMnRt1dmUBd3F6B7EykzKIF7kTE5eue1ctf6vk19znRB22IQ2iZSZpA/9yJjcvTOa+Wu9X2b\n+pwZAqGLgTKD/LkXGZOjd14rd63v29TnzCBoc/mQQ1NmENv/VGZQLve833BEft+mPmfaAmUGygxK\n4O5zX7X69jlyR3/fpj5nuqBtMQhtEykzyJ+7z33V6tvnyB39fZv6nBkCoYuBMoP8ufvcV62+fY7c\n0d+3qc+ZQdDm8iGHVmJm0Nc9ixYZU5r/2de+avXtlRnUd860BcoM0mYGfWUJi4wpzf9UZlAud63v\n29TnTBe0LQahbaKUmUFfWcIiY0rzP5UZlMtd6/s29TkzBEIXg5SZQV9ZwiJjSvM/lRmUy13r+zb1\nOTMI2lw+5NBKzAz6+v3BImNK8z+VGYzP3dffeK71fZv6nGkLlBkoMyjF/1RmMD53SesqbmUGm4Yy\ngzL8T2UG43OXtK7iVmawaSgzKMP/VGYwPndJ6ypuZQbKDDYxpjT/s699db1n0bxtOWYGffn87mWt\nq7iVGSgzWGBMaf7n0By5+vZ9fac/xzmP/r5NPbdd0LYYhLaJlBmU4X8OzZGrb9/Xd/pznPPo79vU\nczsEQhcDZQZl+J9Dc+Tq2/f1nf4c5zz6+zb13A6CNpcPOTRlBrH9z6E5UmYG83z7vr7Tn+Ocj8FR\nK3cXoMxAmUEp/mfkzEDrGlNf6rntgrbFILRNpMygDP8zcmagdY2pL/XcDoHQxUCZQRn+Z+TMQOsa\nU1/quR0EbS4fcmiRMoOuf+dg3rYo/meOmYHu3ZM/R63cXYAygzwzg65ZwrxtUfzPHDODCLqjr2ut\n3F3RthiEtolyzAy6ZgnztkXxP3PMDCLojr6utXIPhdDFIMfMoGuWMG9bFP8zx8wggu7o61or92Bo\nc/mQQ4uUGXT9/cG8bVH8z6E5FvmufwTd0de1Vu4uQJmBMoNS/M/I+mrljq4v9dx2Qdti0ItNZGZX\nmJmb2bapvqvM7IiZfdzMXjbVv9PM7mm23WBm1scxrAdlBmX4n5H11codXV/quR0Cmy4GZrYd2AU8\nMNX3fOCVwNnAxcDbzGxLs/lG4NXAWU27eLPHMAvKDMrwPyPrq5U7ur7UczsI2lw+zGvAzcD5wFFg\nW9N3FXDV1Gs+DFwInArcN9V/GbCvDY8yg9j+Z2R9tXJH15d6btuCMTIDYDfwpubxdDF4C/BjU697\nJ/DDwDJwYKr/xcCH2nApM4jrf0bWVyt3dH2p57YL2haDDW0iMztgZofXabuBq4FrN3VpMp97j5mt\nmtnq2tpa5/HKDMrwPyPrq5U7ur7UczsENiwG7v5Sdz/n+Ab8M/Ac4G4zOwqcAdxpZs8GjgHbp3Zz\nRtN3rHl8fP8s7v3uvuzuy0tLS121KTMoxP+MrK9W7uj6Us/tIGhz+dCm8USb6GzgbuBEJgXjn4Et\nzbbbgQsAA24BLmmz/xoyg1n3LHKP739G1lcrd3R9qee2LRj7dwbTxaB5fg1wP/Bx4OVT/cvA4Wbb\nWwBrs//omcE87uj+Z2R9tXJH15d6brugbTHo7XYU7r7D3R+eev56dz/T3b/J3W+Z6l/1idV0prvv\nbQ52EJSUGczjju5/RtZXK3d0fanndgjo3kSZZAbzuKP7n5H11codXV/quR0EbS4fcmg1ZAbzuKP7\nn5H11codXV/quW0LdG8iZQal+J+R9dXKHV1f6rntgrbFILRNpMygDP8zsr5auaPrSz23QyB0MVBm\nUIb/GVlfrdzR9aWe20HQ5vIhh6bMILb/GVlfrdzR9aWe27ZAmYEyg1L8z8j6auWOri/13HZB22IQ\n2iZSZlCG/xlZX63c0fWlntshELoYKDMow/+MrK9W7uj6Us/tIGhz+ZBDqz0zmHXfoij+Z2R/t1bu\n6PrG4J53v7K2QJlBnMyg699AKM3/jOzv1sodXd8Y3PM+E7qgbTEIbRNFyQxmbYvif0b2d2vljq5v\nDO55nwlDIHQxiJIZdP0bCKX5n5H93Vq5o+sbg3veZ8IgaHP5kEOrPTPo+neTc/Q/o/u74q5L3zzu\nWV5/1/55nwltgTIDZQY5+Z/R/V1x16VvHvcYnwld0LYYhLaJlBnk439G93fFXZe+edxjfCYMgdDF\nQJlBPv5ndH9X3HXpm8c9xmfCIGhz+ZBDU2agzEDcZXGXpq8vn999nM+EtkCZgTKDErzXMTjErXXd\naExJnxVd0bYYhLaJlBnk772OwSFuretGY0r6rBgKoYuBMoP8vdcxOMStdd1oTEmfFYOhzeVDDk2Z\nwfrbhr5nUZ/7iuItizsfjnm+fV/f6c/xs6ILUGZQlg/YJ3dJvu8YHOKOua5RzldlBiOgJB+wT+6S\nfN8xOMQdc12jnK/KDEZAST5gn9wl+b5jcIg75rpGOV+VGSgzGJS7JN93DA5xp1nXvr67v8i9e0o6\nX5UZKDNQZjACh7jTrGtJ50xJ3F3RthiEtolK8gGVGcT0zmvlhrLOmZK4h0LoYlCSD6jMIKZ3Xis3\nlHXOlMQ9GNpcPuTQlBkoMxD3ONwl3aMnyvmqzECZgTKDETjE3Z671nOmJO6uaFsMQttEJfmAygxi\neuelcdd6zpTEPRR6KQZmdoWZuZlta55/j5ndYWb3NP9+99Rrdzb9R8zsBjOzPo5hPZTkAyoziOmd\nl8Zd6zlTEvdgaHP5MK8B24EPA58EtjV93wKc1jw+Bzg29frbgQsAA24BXt6GR5lBtzF93bNokTE5\nckTnXsTPj3CPnlq5u4CxMgPgZuB84OhXisFx2w34PHAicCpw39S2y4B9bXiUGWx+Xzn62mNwROfO\n8f08Bket3F3RthhsyiYys93N//rvnvOyS4E73f1/gNOBB6e2Pdj0zdr/HjNbNbPVtbW1zsdXkg84\nBneOvvYYHNG5c3w/j8FRK/dQ2LAYmNkBMzu8TtsNXA1cO2fs2cBvAa9Z5ODcfb+7L7v78tLSUufx\nJfmAY3Dn6GuPwRGdO8f38xgctXIPhjaXD+s14FzgISb20FHgy8ADwLOb7WcAnwC+fWrMqDaRe1k+\n4BjcOXrqY3Dkyj30PXpS+9oRzpkcubuAsX9nwFRmADwDuBv4oXVed3yAfEmb/SszUGYQjTvCeyr6\nOZMjd1e0LQZD/c5gL/A84Fozu6tpX9dsey3wDuAIcH9TEAZBST6gMoOyfXvduyfmOZMj91DorRi4\n+w53f7h5/Bvu/nR3f8FUe6jZturu57j7me6+t6lcg6AkH1CZQdm+ve7dE/OcyZF7MLS5fMihKTNQ\nZpALt+7dkw9HrdxdgO5NVJYPqMygjMwgx3XNkTu6vtRz2wVti8FQmUEWKMkHHIM7R099DI4+uXNc\n1xy5o+tLPbdDIHQxKMkHVGZQRmaQ47rmyB1dX+q5HQRtLh9yaItmBgcPul9//ZMvsbr297mvlNyH\n9h30v9x1/ROsknn9i4zJkaNP7hzXNUfu6PpSz21b0NIm2jpwrREywxfOvZDbPnchF52b+kjGxTzd\ns7bVOldCpWhTMXJoCpCH44geINe6rqlDzsj6Us9tF6AAuaxQKCV39AC51nVNHXJG1pd6bodA6GJQ\nUiiUkjt6gFzruqYOOSPrSz23Q2DLysrKsAw9Yf/+/St79uzpNGb7dlhagi99CV73OnjFKxbr73Nf\nOXI/a+d2Dh5Z4vPHvsTRH3odL3r94+SztvXVv9GYP/3v7+bvP/tc/vWnf43v+pULF+qvdV1TckfX\nl3puu+C666779MrKyv4NX9jGS8qhKTOoLzOIMLe1ckfXl3puuwBlBmX5gCm5c80MIsxtrdzR9aWe\n2yEQuhiU5AMqMxjuB161rmtqXzuyvtRzOwSUGWTiA0bKDLr4+fO2RZjbWrmj60s9t12gzMDL8gGj\nZAY56hO3MoNI3F2BMoOyfMCU3Clv5BZ9bmvljq4v9dwOgdDFoCQfMEpmkKM+cSsziMQ9FJQZZOID\nps4Muvr5ffn80ee2Vu7o+lLPbRcoM/CyfMBauaPrq5U7ur7Uc9sFKDMoyweslTu6vlq5o+tLPbdD\nIHQxKMkHrJU7ur5auaPrSz23Q0CZQSY+YK3c0fXVyh1dX+q57QJlBl6WD1grd3R9tXJH15d6brsA\nZQZl+YC1ckfXVyt3dH2p53YIhC4GJfmAtXJH11crd3R9qed2CCgzyMQHrJU7ur5auaPrSz23XaDM\nwMvyAWvljq6vVu7o+lLPbRegzKAsH7BW7uj6auWOri/13A6B0MWgJB+wVu7o+mrljq4v9dwOAWUG\nmfiAtXJH11crd3R9qee2C5QZeFk+YK3c0fXVyh1dX+q57QLGzAzM7AozczPbdlz/15vZF8zsF6f6\ndprZPWZ2xMxuMDPr4xjWQ0k+YK3c0fXVyh1dX+q5HQKbLgZmth3YBTywzubfAW45ru9G4NXAWU27\neLPHMAsl+YC1ckfXVyt3dH2p53YQtLl8mNeAm4HzgaPAtqn+VwBvAFaAX2z6TgXum3rNZcC+NjyL\n2ETu7vv2ue/aNfl3M/197kvcdemrlTu6vtRz2xa0tIk2Wwh2A29qHv9/MQC+Grit+Xe6GCwDB6bG\nvxj4UBsuZQYxuaPrq5U7ur7Uc9sFbYvBhjaRmR0ws8PrtN3A1cC16wxbAd7o7l9oeYEyi3uPma2a\n2era2lrn8SX5gLVyR9dXK3d0fanndhC0qRjrNeBc4CEmVwRHgS8zyQ2eDfzNVP+/A58H9jKyTVRS\nta+VO7q+Wrmj60s9t13AGDbRE3Z0XGYw1b9CYxM1z28HLgCMSbh8SZv9KzOIyx1dX63c0fWlntu2\nyLkYLAOHgfuBtwDWZv+6MojJHV1frdzR9aWe2y5oWwx6+Z1BYzftcPeH1+lfcfffnnq+6u7nuPuZ\n7r63OdhBUJIPWCt3dH21ckfXl3puh0BvxSBHlPTd4Vq5o+urlTu6vtRzOwR0b6JM7jdSK3d0fbVy\nR9eXem67QPcm8rJ8wFq5o+urlTu6vtRz2wXo7xmU5QPWyh1dX63c0fWlntshELoYlOQD1sodXV+t\n3NH1pZ7bIaDMIBMfsFbu6Ppq5Y6uL/XcdoEyAy/LB6yVO7q+Wrmj60s9t12AMoOyfMBauaPrq5U7\nur7UczsEQheDknzAWrmj66uVO7q+1HM7BJQZZOID1sodXV+t3NH1pZ7bLlBm4GX5gLVyR9dXK3d0\nfanntgtQZlCWD1grd3R9tXJH15d6bodA6GJQkg9YK3d0fbVyR9eXem6HgDKDTHzAWrmj66uVO7q+\n1HPbBcoMvCwfsFbu6Ppq5Y6uL/XcdgHKDMryAWvljq6vVu7o+lLP7RAIXQxK8gFr5Y6ur1bu6PpS\nz+0QsMlVRP5YXl721dXVzuNuu21SUS+6CC68cPH+Pvcl7rr01codXV/quW0LM7vD3Zc3fF30YiAI\nglAz2haDp4xxMIIgCELeUDEQBEEQVAwEQRAEFQNBEAQBFQNBEAQBFQNBEASBgr5aamZrwCdTH8cC\n2AY8nPogEkC664J054tvcPeljV5UTDEoFWa22uY7vtEg3XVBusuHbCJBEARBxUAQBEFQMRgDG99H\nPCakuy5Id+FQZiAIgiDoykAQBEFQMegVZvYMM7vZzO4zs3vN7EIzO9nM/tzM/qn595mpj7NvmNkv\nmNlHzeywmb3fzL4qqm4z+10ze8jMDk/1zdRqZleZ2REz+7iZvSzNUW8eM3S/oXmvHzKzD5rZM6a2\nhdU9te0KM3Mz2zbVV6xuFYN+8Sbgz9z9m4HzgXuBXwH+wt3PAv6ieR4GZnY68HPAsrufA2wBXklc\n3e8GLj6ub12tZvZ8JnNxdjPmbWa2ZbxD7RXv5sm6/xw4x93PAz4BXAVV6MbMtgO7gAem+orWrWLQ\nE8zsa4HvBN4J4O6PuPu/A7uB9zQvew+w4J+1zhpbgZPMbCvwNOBfCarb3f8a+Pxx3bO07gZ+z93/\nx93/BTgCvGiUA+0Z6+l294+4+5ebp38HnNE8Dq27wRuBK4Hp0LVo3SoG/eE5wBrwLjP7RzN7h5k9\nHXiWu3+6ec1ngGclO8IB4O7HgN9m8j+kTwP/4e4fIbju4zBL6+nAp6Ze92DTFxE/BdzSPA6t28x2\nA8fc/e7jNhWtW8WgP2wFXgjc6O7fAnyR46wRn3x1K9TXtxp/fDeTYnga8HQz+7Hp10TUPQs1af0K\nzOwa4MvA+1Ify9Aws6cBVwPXpj6WvqFi0B8eBB50979vnt/MpDh81sxOBWj+fSjR8Q2FlwL/4u5r\n7v4o8AfAtxFf9zRmaT0GbJ963RlNXxiY2U8A3wdc7o9/Tz2y7jOZ/MfnbjM7ykTbnWb2bArXrWLQ\nE9z9M8CnzOybmq6XAB8D/hj48abvx4E/SnB4Q+IB4AIze5qZGRPd9xJf9zRmaf1j4JVmdqKZPQc4\nC7g9wfENAjO7mIlv/gPu/l9Tm8Lqdvd73P3r3H2Hu+9g8p/AFzbnf9m63V2tpwa8AFgFDgF/CDwT\nOIXJN0z+CTgAnJz6OAfQfR1wH3AYeC9wYlTdwPuZZCOPMvkgeNU8rcA1wP3Ax4GXpz7+nnUfYeKR\n39W0t9eg+7jtR4FtEXTrF8iCIAiCbCJBEARBxUAQBEFAxUAQBEFAxUAQBEFAxUAQBEFAxUAQBEFA\nxUAQBEFAxUAQBEEA/g8qgSed2L7D/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111ebf6a0>"
      ]
     },
     "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",
    "\n",
    "miny = -450\n",
    "maxy = -350\n",
    "minx = 50\n",
    "maxx = 150\n",
    "x = []\n",
    "y = []\n",
    "import desimodel.focalplane\n",
    "import desimodel.footprint\n",
    "\n",
    "while miny < maxy:\n",
    "    startx = minx\n",
    "    while startx < maxx:\n",
    "        x.append(startx)\n",
    "        y.append(miny)\n",
    "        startx += 2\n",
    "    miny += 2\n",
    "\n",
    "#import desimodel.io\n",
    "#tiles = desimodel.io.load_tiles()\n",
    "x = np.asarray(x)\n",
    "y = np.asarray(y)\n",
    "\n",
    "# TYPE HERE TO TEST THE CODE FOR CHANGING GFA BOUNDARIES USING DEGREES_TOLERANCE\n",
    "targetindices = is_on_gfa_test(x, y, .1)\n",
    "#print(targetindices)\n",
    "plt.plot(x, y, 'b.')\n",
    "for gfaid in targetindices:\n",
    "    for indices in gfaid:\n",
    "        plt.plot(x[indices], y[indices], 'r.')\n",
    "\n",
    "\n",
    "        \n",
    "targetindices = is_on_gfa_test(x, y, 0)\n",
    "for gfaid in targetindices:\n",
    "    for indices in gfaid:\n",
    "        plt.plot(x[indices], y[indices], 'g.')\n",
    "        \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 396,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 469.80372211]\n"
     ]
    }
   ],
   "source": [
    "#6.644525362152656, -9.055425745149217 GUARANTEED TO BE IN GFA (RA, DEC)\n",
    "def retrieve_minimum_boundary(x_tolerance, y_tolerance):\n",
    "    import desimodel.footprint\n",
    "    import desimodel.focalplane\n",
    "    \n",
    "    targetx = 116.279135121\n",
    "    targety = -372.885546514\n",
    "    #x, y = desimodel.focalplane.radec2xy(7.11, -10.53, targetx, targety)\n",
    "    #print(x, y)\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",
    "    #print(targetx)\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",
    "        #gfaarea = a1 * a2\n",
    "        #gfaarea2 = a2 * a3\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",
    "        #print(targetarea)\n",
    "        return targetarea\n",
    "        \n",
    "    \"\"\"targetindices is an array 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 targetindices array.\"\"\"\n",
    "    #print(gfatable['X'], gfatable['Y'])\n",
    "    #return targetindices\n",
    "\n",
    "area = retrieve_minimum_boundary(0, 0)\n",
    "print(area)\n",
    "# CODE IS QUITE REDUNDANT AND SHOULD BE CLEANED UP BUT IS USED FOR TESTING RIGHT NOW"
   ]
  }
 ],
 "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
}