ellipse_and_stick.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import linalg\n",
    "from scipy.special import gamma, factorial, binom, comb,i0,iv,ive\n",
    "\n",
    "## Quadrature ellipse given mean vector, and squeezing parameter d=np.exp(2*np.abs(z)), thet is argument of z\n",
    "\n",
    "#Mean vector\n",
    "m=np.array([-.1,.8])\n",
    "\n",
    "#Consider covariance matrix defined by \n",
    "# gg=(d/2)+(1/(2*d))\n",
    "# hh=(d/2)-(1/(2*d))\n",
    "# xx=(1/2)*np.array([[gg-(hh*np.cos(thet)),-np.sin(thet)*hh],[-np.sin(thet)*hh[0],gg+(hh[0]*np.cos(thet))]])\n",
    "\n",
    "#Can solve for d and thet to proceed to plot. C XXVII p.104. If given pure state covariance\n",
    "#matrix [[a,b],[b,((b**2)+(1/4))/a]], get thet = np.arcsin(np.abs(b)/((1/2)*gg)) and\n",
    "# d=(1/(2*a))*\n",
    "#   ((2*(np.abs(b)**2))+(2*(a**2))+(1/2)+np.sqrt( (((2*(np.abs(b)**2)) + (1/2) +(2*(a**2)))**2)-(4*(a**2)) ))\n",
    "\n",
    "d=3\n",
    "## Negative correlation goes with positive thet in our convention\n",
    "thet=.8\n",
    "\n",
    "dats=[]\n",
    "dats.append([m,d,thet])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.1,  0.8])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dats[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAJDCAYAAADuENvfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXBc9d3n+8/3nN7UrW7t8iJDIOxrIHHCviUsBsxiCAnkgYQsmDz3SdUzVTM1M7mpmltTU1N3qp6quXOr5rkPmD0JIQtewGwGAgRICMGAIRCzBQhYXrTYVmuX+vTv/qFjYoMtdUu9SXq/qlyFrHOOvuki+O3TfX4/c84JAAAAklftAQAAAGoFYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAoZKFkZn5ZvaqmT1UqmsCAABUUinvGP2zpM0lvB4AAEBFlSSMzGyJpEsl3V6K6wEAAFRDqe4Y/S9J/1FSvkTXAwAAqLjITC9gZssldTnnXjazcyc5bqWklZKUSqW+dPTRR8/0RwMAABTs5Zdf7nHOtU12jM10rzQz+78l3SApJykhKSNpjXPu+gOds3TpUrdx48YZ/VwAAIBimNnLzrmlkx0z47fSnHM/ds4tcc4dIulaSU9NFkUAAAC1inWMAAAAQjP+jNHenHPPSHqmlNcEAACoFO4YAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQmnEYmVnCzP5kZq+Z2Ztm9l9LMRgAAEClRUpwjVFJX3XODZhZVNLzZvaoc+6PJbg2AABAxcw4jJxzTtJA+GU0/OVmel0AAIBKK8lnjMzMN7NNkrokPeGce7EU1wUAAKikkoSRcy5wzp0kaYmkr5jZ8Z8+xsxWmtlGM9vY3d1dih8LAABQUiV9Ks05t1vSM5KW7ed7q5xzS51zS9va2kr5YwEAAEqiFE+ltZlZY/jPdZLOl/TWTK8LAABQaaV4Km2RpHvMzNdEaP3aOfdQCa4LAABQUaV4Ku11SSeXYBYAAICqYuVrAACAEGEEAAAQKsVnjABJ0trBp9WV26XlqTPVEWmv9jgAABSNMELJJBTTy2Ob9UFuq85KfEEXJE9VzKLVHgsAgILxVhpK5py6L2mJv0ANXkrPjWzSbdm12hn0VXssAAAKRhihZJJeQtfWX6jd+X4t8JvVHezS/9t3n/4y+n61RwMAoCCEEUrqyNjBuiR5hrYHvWrxGpSyOt098JAeHnxOY2682uMBADApwggld1biZH0hdqS2B71Kegl1+G16bmST7squ10B+qNrjAQBwQIQRSs4zTytS56nFb9TOICvfPHX4bdqS26F/67tfXcHOao8IAMB+EUYoizovrhvSlyhQoMH8sMxMCyLNGnGj+te+3+j98c5qjwgAwGcQRiibNr9J304v1+78gEbcmCSpyc8oaXHdll2rjSN/kXOuylMCAPB3hBHK6rDoEl1Xf5G6g10adzlJUsqrU6vfqN8MPqkNQy8ocEGVpwQAYAJhhLI7MX6ELkuerW1BzycRFLeoFvutenpko+4b2KCR/GiVpwQAgDBChZyR+IK+mviKtgY9yodvn/nma4nfrr+MfaDb+9dpd9Bf5SkBAPMdYYSKMDNdmDxFS+PHalvQ/clni8xMiyOt6gp26//L/kadue4qTwoAmM8II1SMZ56uTJ2rI6Of0/agd5/vtfuNcnK6pe9+bR79oEoTAgDmO8IIFRW1iK6rv0gdkXZ15Xbt870Gr15pL6l7Bh7Sc8Ov8sQaAKDiCCNUXCJc4yjjp9T7qU1mk15CC/xmrR96Tg8NPccTawCAiiKMUBVpL6Xvpi+XL199+YF9vhe1iDr8Nj0/8pp+NfAke6wBACqGMELVtPgN+l7mco26cQ3mh/f5nm+elvht+vPYu7o7+9Bnvg8AQDkQRqiqxZE23Zi+bJ/VsfcwMy32W/VRsF2rsmu0i8f5AQBlRhih6j4f7dC36pepO9il0U+9bWZmWug3a3d+QP+WvV/bc70HuAoAADNHGKEmnBA/XNekLtCOYKfGwq1D9tbmNyrvAv1b9jf6YHxrFSYEAMwHhBFqxtLEMboyeY62Bz3K7edptEY/rTqL67bsGr0x+l4VJgQAzHWEEWrKaYkTdXHydG3da1+1vdV7STV7Gf1s4BG9MPI6ax0BAEqKMEJNMTOdm1iqrya+rM6gR4HLf+aYhBfXAr9Fawef0YahPyq/n2MAAJgOwgg1Z8++amcmTtLWoPuTTWf3FrOIOvxWPT3yklYPPqXx/XwuCQCAYhFGqEmeebo0eYaWxo/V1r02nd2bb746/Ha9PLpZP+9/RMP50SpMCgCYSwgj1CzffK1InacTY4dra9Cz3zjyzNTht+m93Bbd0b9O2fxgFSYFAMwVhBFqWsR8fb3+fB0ZPVjbgt79xpGZaZHfoq7cLt3Sd7+6g137uRIAAFMjjFDzYhbVt9LL9LnIIu3I7zzgce2RJo26Mf1b3/36OLejghMCAOYKwgizQtxiuiF9iRb6LdqRO3AcNfkZRczXLX2r9fbYh5UbEAAwJxBGmDWSXkI3pi9Ts59RV7D7gMdlvJQavJTu6l+vV0fequCEAIDZjjDCrFLvJfW9zBVKW516g74DHpf0EmrzG3Xf4ON6bvhVFoIEABSEMMKs0+DV6/uZKxS1qHYF/Qc8Lm4xLfJbtH7oWT3OQpAAgAIQRpiVmv0G/SBzpSSpLz9wwOOiFlGH36anRl7SusFn9rsHGwAAexBGmLXa/SZ9P3OFxlxu0vWLJhaCbNOLo2/qlwMbNObGKzglAGA2IYwwqy2OtOkHmSs17EY1kB864HGeeVrit+nNsb/q7uxDGsqPVHBKAMBsQRhh1jsoskDfS1+h/vywBvPDBzzOzLTYb9NHwTbdnl036VtwAID5iTDCnHBodLG+l7lcffnBSe8GmZkW+i3qye/Wqr416pnksX8AwPxDGGHOOCy6RDeml2t3vn/KDWXb/SYNu1Hd0ne/tua6KzQhAKDWEUaYU46MfU431F+qnfmsRtzYpMc2+xmZmW7NrtaH49sqNCEAoJYRRphzjokfqm/VL1NPsFujU8RRg1evOotrVXYNW4gAAAgjzE0nxA/XtfUXqjvYrdEpHs+v95Jq9Op1V/96bRp9u0ITAgBqUaTaAwDlclL8KOVcoN8M/lYL/CbFLHrAY5NeQp55um9gg4bdqE6NnyAzq+C0AIBawB0jzGlLE8dqRepc7Qh2atzlJj02YTEt8Fu0bvAZPTX8EvurAcA8RBhhzjs1cYIuS56tbUHPlHEUs4gW+a3aMPxHPTT0nAK2EAGAeYUwwrxwRuILujR5prYFvVPulxYJtxB5fmSTVg8+NWVMAQDmDsII84KZ6ezEF3Vx8jRtDXqmjCPfPHX47Xp19C39YuCxKZ9uAwDMDYQR5g0z07mJpbqw7hRtDXqmfJvMC7cQeXvsQ92dfWjS7UYAAHMDYYR5xcz0tbqv6PzEl9UZ9Chw+SmPX+S36uNgu27PrlM2P1ihSQEA1UAYYd4xM12QPFXnJb6kzlx3QXG00G9Rb9CnW/tWqzfoq9CkAIBKI4wwL5mZliVP1zl1X1RnMHUcSVJ7pElDbkS3ZFdre663AlMCACqNMMK8ZWa6OHm6zkqcpK1Bt/IFxFGL3yDn8role78+ym2vwJQAgEoijDCveebp0uSZOj1xojoLjKNGP62YRbWqb43eHfuoAlMCACqFMMK855mn5cmzdEr8eHUGPcoXsOJ1xksp7aV0Z/+D+vPoexWYEgBQCYQRIMk3X1ekztHS+DHqDLoKiqOUl1CL36B7Bx7Vn0beZAsRAJgDCCMg5JuvFanz9MX40doadBcUOgmLqd1v0urB3+p3I68QRwAwyxFGwF4i5uvq1Nd0YuxIdRYYRzGLapHfqkeGfq9Hh/5Q0OeUAAC1acZhZGYHmdnTZrbZzN40s38uxWBAtUTM1zX1X9OJscMLjqM9+6v9buRlrR18ZsotRwAAtakUd4xykv69c+4YSadK+iczO7YE1wWqJmoRXVN/gY6LfV5bg56C4mhif7U2vTT6pn45sEFjbrwCkwIASmnGYeSc2+aceyX8535JmyV1zPS6QLVFLaJv1l+oo2OHFBxHXhhHb469r5/2P6Th/GgFJgUAlEpJP2NkZodIOlnSi6W8LlAtMYvquvqLdGT04ILjyMy02G/VB7mtuqv/QTafBYBZpGRhZGb1klZL+nfOuex+vr/SzDaa2cbu7u5S/Vig7GIW1bfSy3R49CBtC3oLjqNFfqu2BT1sPgsAs0hJwsjMopqIonudc2v2d4xzbpVzbqlzbmlbW1spfixQMXGL6fr0xTo0sljbC4wjSVrgN6s36NOq7BrtCvrLPCUAYKZK8VSaSbpD0mbn3P+c+UhAbYpbTDekL9HBkUXani88jtojTRrID+vW7Gr1BLvLPCUAYCZKccfoDEk3SPqqmW0Kf11SgusCNSfhxfXt9KXq8BdoR35nwee1+g0ac+O6pW+1tud6yzghAGAmSvFU2vPOOXPOneicOyn89UgphgNqUZ0X143p5Vrkt2pHUHgcNfsZSU63ZlerM9dVvgEBANPGytfANCS9hG5MX6Z2v1k7coXHUaOfVtQiujW7Rh+ObyvjhACA6SCMgGlKeXX6XvpytUeairpzlPFSSlpct2XX6t2xj8o4IQCgWIQRMAMTcXSFFvjNRcVRvZdUxkvpzv4H9ZfR98s4IQCgGIQRMEMpr07fTV+uBX6zthcRRykvoWYvo58OPKxNo2+XcUIAQKEII6AE9sTRwiLjqM6Lq9Vv1H0DG/SnkTfLOCEAoBCEEVAie+Josd+q7UHhj+QnLKYFfrPuH/ytnht+teD1kQAApUcYASWU8ur0nfRyLfbbilohO2ZRLfJbtH7oWT01/BJxBABVQhgBJZby6nRj+jIt9tu1I7+z4MiJWkSL/TZtGH5Bjw79QXmXL/OkAIBPI4yAMphY52h5GEeF3zmKmK8Ov12/G3lZDw49q8AFZZ4UALA3wggokz1xtMRfUNTear556vDb9MLIn7V68CmNu1yZJwUA7EEYAWWU9BL6TvoyHeQv0PYi3lbzwjh6ZfQt/WrgcY258TJPCgCQCCOg7Oq8uL6TvkwH+wuLjCNTh9+mN8fe18/7H9VIfrTMkwIACCOgAuq8uL6dvnQijop4Ws3MtNhv1XvjH+ue/oc1lB8p86QAML8RRkCFTNw5ulSfiywuOo4W+S36ONimO/sf0EB+qMyTAsD8RRgBFZTw4vp2+pJpxdFCv1U7cjt1W3at+vIDZZ4UAOYnwgiosD1xdGh0sbYVEUeStCDSrF35fq3qW6OdQV8ZpwSA+YkwAqog4cV1ff0l+vw04qjdb9KQG9Et2TXqCnaVcUoAmH8II6BKEl5cN6Qv1WHRJdoW9BQVRy1+gwKX0619q7Ut11PGKQFgfiGMgCqKW0zXpy/WYdGDtLXIOGryMzIz3ZpdrY9zO8o4JQDMH4QRUGV74uiI6MFFx1GjV6+YRbUqu0YfjHeWcUoAmB8II6AGzCSOMl5KKavTbdl1emfsb2WcEgDmPsIIqBExi+r69MU6chpxVO/VqdGr11396/XG6HtlnBIA5jbCCKghMYvqH9IX66jY54qOo6SXUIuX0c8HHtUrI5vLOCUAzF2EEVBjYhbVt+qX6ejYIeoMuouKo4QXV5vfqF8NPqEXRv5c1LkAAMIIqEkxi+q6+ot0TOzQou8cxS2mBX6z1g48rWdHXiGOAKAIhBFQo2YSRzGLalGkRQ8PPa8nh18kjgCgQIQRUMNiFtW19Rfq2NihRb+tFrWIOvw2PTH8Jz089LzyLl/GSQFgbiCMgBoXs6i+WX+hjot9XluLjCPffHX4bXpu5FWtG3xGOReUcVIAmP0II2AWmLhzdJGOix1W9J0j3zx1+G16cfQN3T/wpMZdroyTAsDsRhgBs0TUIvpm/YU6IXZ40XHkmaclfrs2jb2j+wY2aMyNl3FSAJi9CCNgFolaRN+ov0Anxg4v+m01M1OH36bNYx/oZ/2PaCQ/WsZJAWB2IoyAWSZqEV1Tf4FOiB1R9J0jM9Niv1Xvj2/Rnf3rNZgfLuOkADD7EEbALDQRR+frxNiR04qjhX6LtgbduiP7gPrzg2WcFABmF8IImKUm3lY7X1+Ydhw1qzvYpduya7Ur6C/jpAAwexBGwCwWMV/X1J+vk8I4yhe5kOOCSLP68oNalV2t3qCvTFMCwOxBGAGzXMR8fb3+fJ0cO2pacdTmN2rEjenW7Gp1BTvLNCUAzA6EETAHRMzX1fVf05fiR08rjlr8BgUu0C19q9WZ6y7TlABQ+wgjYI6ImK+rUl8N46ir6Dhq8jPyzdOq7Gp9lNtepikBoLYRRsAcsieOlsaPnVYcNXj1Slhct/Wt1V/Ht5RpSgCoXYQRMMdMxNF5+nL8uGnFUdpLqt6r0x3ZB/TW6IflGRIAahRhBMxBvvlakTp3rzjKF3V+yqtTo1evewYe0uuj75ZpSgCoPZFqDwCgPPbEkcn0p9E31OG3ybPC/y6U9BLyzNO9A49qzI1raeLYMk4LALWBO0bAHOabrytT5+iU+PHh02rF3TlKWEztfrN+PfiEnh/eVNQikgAwGxFGwBznm68rZhBHcYtqod+qB4d+p2dGNhJHAOY0wgiYB/bE0WmJE6cVRzGLaLHfpkeHXtCGoT8SRwDmLMIImCd883V58uxpx1HEfHX4rXpq+CU9MvR80ecDwGxAGAHziGeeLk+erdPDOAqKjBvffHVE2vTsyKtaP/ScAheUaVIAqA7CCJhnPPN0WfJsnZE4SVunFUeeFvtt+sPIa1o3+DvliCMAcwhhBMxDnnlanjxzRnHU4bfpxdE3tHrgt8QRgDmDMALmqT1xdOY048gzT0v8dr069rZ+PfC4xtx4mSYFgMohjIB5zDNPlybP1FmJk6f1mSPPTB1+m14fe0+/JI4AzAGEETDPeebpkuQZOjtxsjpzXUXHkYVxtHnsA/28/1GNurEyTQoA5UcYAfgkjs6p+5I6c91FP21mZlrst+q98Y/00+zDGs6PlmlSACgvwgiApL/H0Xl1XwrfVis+jhb5rfowt1X39K/XUH6kTJMCQPkQRgA+YWZaljxdF9Sdos6gp+inzcxMC/0WbQm6dGf/AxrMD5dpUgAoD8IIwD7MTOfXnaJlydO0NejRuMsVff5Cv0Xbc726PbtO/fnBMk0KAKVHGAH4DDPTeYmlujR5hrZNI44kaWGkRT3Bbt2WXae+/EAZpgSA0iOMAOyXmensxBd1efIcbQt6p/Uo/oJIs3bn+3Vbdq12Bf1lmBIASoswAnBAZqYz607SitR52j7NOGr3mzSQH9aq7Gr1Bn1lmBIASqckYWRmd5pZl5m9UYrrAagtpyVO0DWp87Uj2KnRacRRq9+gETemW7Nr1B3sKsOEAFAapbpjdLekZSW6FoAa9OXEcfpm6gJ1BTs1Mo1FHFv8BgUup1v71mhHrrcMEwLAzJUkjJxzz0raWYprAahdX0wco+vqL1JPsFsj01jEscnPyMnp1uwabc11l2FCAJgZPmMEoCgnxY/S9fUXqzefndYK101+Wp55WpVdo85cVxkmBIDpq1gYmdlKM9toZhu7u/mbIjCbHR8/XN+uv1Q789lprXDd6NUrZlHdml2jv41vK8OEADA9FQsj59wq59xS59zStra2Sv1YAGVyTPxQ3Zhert35fg1OI44yXkp1Ftft/ev0wXhnGSYEgOLxVhqAaTsqdoi+m75cffmBaW3/kfaSSlmdbs8+oHfHPirDhABQnFI9rn+fpBckHWVmW8zs+6W4LoDad0TsYN2UuVL9+SH154eKPr/eq1PGS+mu/gf19tiHpR8QAIpQqqfSrnPOLXLORZ1zS5xzd5TiugBmh0OjHfpBZoWG3aiy09gbLeUl1OildXf/Q3pz9K9lmBAACsNbaQBK4pDoIt2UWaFRNz6tOEp6CbV4Gf1s4BG9NvpOGSYEgKkRRgBK5qDIAq3MrNCYy01r49iEF1er36hfDDyml0c2l2FCAJgcYQSgpDoi7VqZWaHA5bV7GhvHJiymdr9Zvx58Qi+OsMsQgMoijACU3OJIm1ZmrpKTtCvIFn1+3KJa4Ddr9eBTen54k5xzpR8SAPaDMAJQFgsjLbq54Wp55mvnNOIoZlEt8lv04NDv9OzIK8QRgIogjACUTbvfpJWZFYpYRL1BX9HnRy2iRX6rHh76vZ4e2UgcASg7wghAWbX5Tbo5c5USFlPPNONosd+qx4Ze0JPDLxJHAMqKMAJQdi1+g1ZmrlLKEuoOdk967MiINNBv+/xexHx1+K16YvhFbRh6QXmXL+e4AOYxwghARTT5Gd2UWaGMl1JXsOuAx/3y53V6YkP8M7/vm68Ov11PD2/UI0O/J44AlAVhBKBiGv20bsqsUJOXVldu/3HU0OCU7bP9fs83T4sjbXpu5FWtH3pOgQvKOS6AeYgwAlBRGS+lH2SuVIvfoO253s9+vyGvvj5PB/ookW+eFvtt+sPIa3pg8HfEEYCSIowAVFw6jKOFkRZtD3r3+UB1Q6PT2KhpdOTA5/vmqcNv04ujb2jN4NPKEUcASoQwAlAVKa9O301frsV+u7bn/x5HmYaJzw719U3+nyfPPHX47do4uln3DzypcZcr+8wA5j7CCEDVTMTRZTrYX/TJnaOGholAOtDnjPbmmWmJ36ZNY+/oNwNPEEcAZowwAlBVdV5c30lfqs9Hl2hr0KN0ZuJtsanuGO1hZurw2/T62Hv65cAGjbnxco4LYI4jjABUXcKL64b0JTo6doh2RnsUix34ybT92RNHb469r3v7HyWOAEwbYQSgJsQsquvqL9Lx8c+rLjOmvt2Fh5H09zh6Z/wj/az/EY3kR8s0KYC5jDACUDNiFtU36y/U4qaEevryyhe5/YeZabHfqr+Ob9FP+x/RMHEEoEiEEYCaErWIvtx6kEazcW3JdRW9wrWZaZHfor/ltuqe/vUayk/y3D8AfAphBKDmHNyc0sio9AV3vDqDbgXTiKOFfou2BF26u3+9BvPDZZoUwFxDGAGoOUua6iRJJ46drDMTJ6kz1zXtONoa9OjO/gc0kB8qx6gA5hjCCEDNWdKUlCR17hrVpckzdV7dUnUGXdNa4Xqh36yu3C7dnn1A2fxgqUcFMMcQRgBqzp47Rlt2DckzT8uSp+uCulO0NeiZVhwtiDRrZ7Bbt2fXqS8/UOpxAcwhhBGAmtNQF1V9PKItuyY+G2RmOr/uFF2cPE1bg55prXDdHmlWX75ft2XXalfQX+qRAcwRhBGAmmNmWtJU90kY7fm9cxNLdVnyLG0LejU2jThq85vUnx/S7dm12hn0lXJkAHMEYQSgJk2E0b4fmDYznVV3slakztP2oHdaK1y3+Y0aciNalV2rXuIIwKcQRgBq0pKmpDp3DcvtZ5HH0xIn6JrU17Qj2KlRN1b0tVv8Bo26Md2aXa3uYFcpxgUwRxBGAGrSkqY69Y/mlB3e/1tmX04cp2tTF6or2D2t7T9a/AYFLtCtfWvUFeyc6bgA5gjCCEBN2vNk2se7Drz+0MmJo3V9/cXqzWentf1Hk5+Rk9MtfWu0Ldcz7VkBzB2EEYCatGcto70/gL0/J8QP17frL9WufHZa2380+Wl5klZl16oz1z2dUQHMIYQRgJq091pGUzkmfqhuTF+m3fmBaW3/0einFTFPq7JrtCW3o+jzAcwdhBGAmvTptYymcmTsc/pB5kr154fUP43tPxq8esUtqtuya/VRbnvR5wOYGwgjADVpf2sZTeXz0Q79ILNCw250Wtt/ZLyUEhbXbdm1+nB8W9HnA5j9CCMANWt/axlN5ZDoIt2UWaFRNzat7T/SXlIpS+iO7Dp9ML616PMBzG6EEYCaNdlaRpM5KLJAKzNXK+fy2j2N7T/qvaRSXl0YR51Fnw9g9iKMANSsqdYymkxHpE03Z66Sk7QzyBZ9fr1Xp3ovqTuyDxBHwDxCGAGoWYWsZTSZhZEW3dxwtSLmTzuO0l5St2cf0PvEETAvEEYAalahaxlNpt1v0srMVYpZVD3B7qLPT3l1yoR3jogjYO4jjADUrGLWMppMq9+olZkVSlmdumcYR38d3zKjWQDUNsIIQM0qdi2jyTT7Dbops0IZL6WuXPEbx+6JozuzD+rdsY9mPA+A2kQYAahZ01nLaDKNflo3ZVaoyc9oR674jWMn4iilu/rXE0fAHEUYAahp01nLaDIZL6WbMldqQaRZ24OeopcCSHkJNXj1uqt/vd4hjoA5hzACUNOmu5bRZOq9pL6XvkIH+Yu0bQZxdHf/er0z9reSzQWg+ggjADVtJmsZTSbpJXRjZrmOiB6srdOMo0biCJhzCCMANW2maxlNJm4x/UP6Yh0fO0ydQbfyRcZR0kuo0Uvrrv71envsw5LPB6DyCCMANa0UaxlNJmZRfbP+Qn0pfow6gy7lXb6o85NeQk1eWvf0P6S3Rj8sy4wAKocwAlDTSrWW0WQi5uuq1Hk6PXGitgTdCqYRR41eWvcMEEfAbEcYAahppVzLaDK++bosebbOS3xJnUGXci4o6vykl1Czl9E9Aw9p8+gHZZoSQLkRRgBqWqnXMpqMZ56WJU/XsuTp2hr0aNwV94HvOi+uZi+jnw48rL+Mvl+mKQGUE2EEoOaVei2jyZiZzkss1WXJs7Qt6NWYG34Jj14AABomSURBVC/q/L3j6E3iCJh1CCMANa8caxlNxsx0Vt3Juib1Ne0IdmrEjRV1fp0XV4vXoJ8PPKI3R/9apikBlANhBKDmlWsto6l8OXGcrqu/SD3Bbg3nR4s6dyKOMvr5wKN6Y/S9Mk0IoNQIIwA1r5xrGU3lpPhR+k79cu3KZzWYL+5zTgniCJh1CCMANa/caxlN5Zj4ofpu+nL154fUny8uzhJeXK1eg34+8Kj+TBwBNY8wAlDzKrGW0VSOiB2smzIrNOJG1ZcfKOrchBdXq9+oewce1euj75ZpQgClQBgBqHmVWstoKp+LLtLKzNXKubx2Bf1FnZuwmFr9Rv1i4DHiCKhhhBGAmlfJtYym0hFp0w8zV8kzTzuDbFHn7omjewce06bRt8s0IYCZIIwAzAqVXMtoKgsiLbo5c5XiFlV3sLuocxMWU5vfqF8OPE4cATWIMAIwK1R6LaOptPqNWpm5Wmkvqa7czqLO3XPn6L6BDXp15K0yTQhgOkoSRma2zMzeNrP3zOw/l+KaALC3aq1lNJkmP62VmavU4jdqR7CzqGhLWEztfpN+Ofg4cQTUkBmHkZn5kv5V0sWSjpV0nZkdO9PrAsDeqrmW0WQyXko/yFypRX6rtueLi6N4GEe/GnxCr4xsLuOUAApVijtGX5H0nnPufefcmKRfSrqiBNcFgE9Uey2jyaS8On03fbkOjSzS1qCn6DhqC+PoZeIIqLpShFGHpI/3+npL+HsAUDK1sJbRZOq8uL6dXq6jY4eoM+guMo6iaveb9evBJ7SROAKqqhRhZPv5vc/8F8HMVprZRjPb2N3dXYIfC2A+aaiL6pRDm5Wpi1Z7lAOKWVTfql+mE2NHakvQrfw04ug3g09o48hfyjglgMlESnCNLZIO2uvrJZK2fvog59wqSaskaenSpbXxWAmAWcPM9KubT6v2GFOKWkTfqD9f8cGo/jT6hhb7bfKtsL+Dxi2qBX6zfjP4pJycvpw4rszTAvi0UtwxeknSEWZ2qJnFJF0r6cESXBcAZqWI+VqROldnJk7S1ly3AhcUfG4sjKP7B5/Sn0beLOOUAPZnxneMnHM5M/uRpA2SfEl3Ouf4fzOAec0zT8uTZymmqH478pIW+62KmF/QuRNx1KT7B38rJ6dTEseXeVoAe5TirTQ55x6R9EgprgUAc4WZ6cLkqYp7MT0y9LwW+q2KWWH/2Y1ZVAv9Zq0efEpOTqcmTijztAAkVr4GgLIyM52T+KJWpM7T9qBXo2684HP3xNGawaf1wsifyzglgD0IIwAoMzPTaYkT9c3UBeoKdmrEjRV87p44Wjv4tF4Yeb2MUwKQCCMAqJgvJY7R9fWXqDfo01B+pODzJuKoResGn9Efhl8r44QACCMAqKAT4ofrxvRy7c4PaDBf+CreMYtogd+iB4Z+RxwBZUQYAUCFHRU7RD/IXKkBN6z+fOEree8dR88PbypqdW0AhSGMAKAKPh/t0E2ZFRp1Y+rLDxR8XswiWui36MHBZ/X8CHEElBphBABVcnBkoVZmrlLgnHYF2YLPi1pEiyItemjoOeIIKDHCCACqaHGkTT9suFq++eoN+go+LxreOVo/xJ0joJQIIwCosna/STdnrlbSEuoOdhd8XtQiWuS3av3Qs3pu5FXiCCgBwggAakCL36CbMivU4NWrK7er4PP2xNFDQ8/r2ZFXiCNghggjAKgRjX5aKzMr1BZp0vagt+DImYijFj0y9Hv9jjgCZoQwAoAaUu8l9f30FVrit08jjlqJI2CGCCMAqDFJL6Eb05fpsOhB2hr0FBw5EfO1OIyjZ0Y2EkfANBBGAFCDEl5c16cv1rGxQ9UZdCtfZBw9OvQHPU0cAUUjjACgRsUsqmvrL9LJsaPUGXQp7/IFnTcRR216bOgFPTX8EnEEFIEwAoAaFrWIvl7/NZ0aP0GdQbeCouKoVRuG/0gcAUUgjACgxvnm64rUOTo78UV1Bt3KuaCg8/bE0ePEEVAwwggAZgHPPF2SPEMX1J2ibUGPxl2uoPMi5muR36oNwy/omZGXiSNgCoQRAMwSZqbz676i5cmztC3o1ZgbL+i8v3/m6A+skA1MgTACgFnEzHRW3cm6KnWedgQ7NerGCjovYr4W+mw8C0yFMAKAWejUxAm6NnWhuoPdGsmPFnTO37cPeU4vjLxe5gmB2YkwAoBZ6uTE0bqh/hL15rMayo8UdE7UIlrgt2jd0O/0x5E/l3lCYPYhjABgFjsufpi+m75cfflBDeSHCjonZhEt9Ju1ZvBpvTTyZpknBGYXwggAZrkjYwfrB5krNeRG1ZcfKOicmEW10G/W/YO/1csjm8s8ITB7EEYAMAccGl2smzNXK3BOu4L+gs6JWVTtfrN+PfiENo2+XeYJgdmBMAKAOaIj0qYfNlytqEXUHewu6Jy4RdXuN+m+gcf1+ui7ZZ4QqH2EEQDMIe1+k27OXK0Gr147cjsLOiduMbX5jfrFwGN6Y/S9Mk8I1DbCCADmmCY/rZWZFVoUadW2oLegNYsSFlOr16B7Bx7TX0bfr8CUQG0ijABgDqr3kvpu+nJ9PtKhzqC7sDjy4mr2MvrZwCN6e+zD8g8J1CDCCADmqDovrhvSl+jE2OHaEnQr7/IFndPkpXVP/0N6Z+yjCkwJ1BbCCADmsJhF9Y36C3VK/Dh1Bt0KCoijpJdQo5fW3f3r9d7YxxWYEqgdhBEAzHER83Vl6lydl1iqzqBb4y435TlJL6EGL6W7+tfrg/HOCkwJ1AbCCADmAc88XZQ8TcuTZ2p70KsxNz7lOSmvThkvqTuyD+jD8W0VmBKoPsIIAOYJM9PZdV/U11Nf045gZ0Gbz6a8OtV7dbqjf50+ym2vwJRAdRFGADDPfDlxnG6ov0Q781kN5oenPL7eSyppCd2eXaePczsqMCFQPYQRAMxDx8cP1/fSV2jADSubH5zy+LSXVMJiuj27Vp25rgpMCFQHYQQA89ThsYO0MnOVxl2uoP3VMl5KcYvptuxabc11V2BCoPIIIwCYxw6KLNAPM1fLN0+9Qd+Ux2e8lKIW0e3ZddqW66nAhEBlEUYAMM8tiLToh5mrlfLq1FXA/moNXr08mW7LrlNXUNh+bMBsQRgBANTsN2hl5iq1RZq0vYD91Rr9tExOt/WtVXewq0JTAuVHGAEAJE28Tfb99JU62F+krUHPlHHU5GcUKK/bsusKehsOmA0IIwDAJ5JeQt/JXKpjYoeqM+hWfoo4avYzGnfjuj27tqAPcAO1jjACAOwjbjFdV3+RvhQ/Rp1B15T7q7X4DRp2Y7qjf512E0eY5QgjAMBnRC2iq1Ln6azEyeoMupRzwaTHt/oNyuaHdGf/gwWtiwTUKsIIALBfvvm6NHmmliVP19Zcz5T7q7X7jdodZHVX9kEN5IcqNCVQWoQRAOCAzEznJZbq6/WF7a/WHmlWT7Bbd/evL2i7EaDWEEYAgEmZmb5SxP5qCyLN2h7s1D39D2u4gI1qgVpCGAEACnJ8/HB9P3NlQfurLfCa1Bns0E/7H57yLhNQSwgjAEDBDosu0Q8zVyvnAu0Ksgc8zsy00GvRR7lt+sXAhik/nwTUCsIIAFCUjki7ftjwdUUtqq5g9wGPMzMt9Fv07vhH+vXAExp3uQpOCUwPYQQAKFq736QfZq5Wk5fWjmDnAVfJNjMt9lv1xthftXbwaQVTPPYPVBthBACYlkY/rZWZFerw27Rtki1EJuKoTS+Pbtb6oeeUn2LBSKCaCCMAwLSlvDp9N3O5joodMukWIp6ZOvw2vTDyuh4bemHKfdiAaiGMAAAzEreYvlW/bK8tRPb/dplnnhb7bXpmZKOeGn6JOEJNIowAADMWtYiuTn1V5yW+rM6gR2MH+KC1H8bRhuE/6vmRTRWeEpgaYQQAKAnPPF2UPFUrUueqa5JVsiPma5HfovVDz+mlkTcrPCUwOcIIAFAyZqbTEifq+vqLtTOfPeCeaVGLaKHfrNWDT+m10XcqPCVwYIQRAKDkjo8frpsyV2nYjWp30L/fY2IWVavfqPsGNmjz6AcVnhDYP8IIAFAWh0YX6x8z18gzT90HWAgyYTG1eA362cDD+uv4lgpPCHzWjMLIzK4xszfNLG9mS0s1FABgblgYadE/Zq5Ro1ev7QdYCLLOi6vBq9fd2fX6KLe9ClMCfzfTO0ZvSLpK0rMlmAUAMAc1+WmtzFylg/2F2nqAhSBTXp2SXkJ3ZNdpa667ClMCE2YURs65zc65t0s1DABgbkp5dboxs1wnxI7QlqBLwX5Wv057ScUsqtuz69QV7KzClACfMQIAVEjMovpm/fk6M3GSOoOu/W4q2+DVyyTdnl2n3qCv8kNi3psyjMzsSTN7Yz+/rijmB5nZSjPbaGYbu7u5TQoA85FvvpYnz9Ly5FnaFvRq1I195pgmP6Nxl9Od2QcO+EQbUC5WiiXZzewZSf/BObexkOOXLl3qNm4s6FAAwBz16shb+vXgk2r00kp5ic98vyvYrSZvYqPalFdXhQkx15jZy865SR8W4600AEBVnJw4Wt9NX64BN6S+/MBnvt/uN6o3v1v39j+mMTdehQkxH830cf0VZrZF0mmSHjazDaUZCwAwHxwZO1g/zFytvHP7/UzRAq9ZH+Q6df/AkwfcnBYopZk+lbbWObfEORd3zi1wzl1UqsEAAPNDR6Rd/0fDNUp5ddqR23etIzPTYr9Vr429q4eHfr/fR/2BUuKtNABA1bX4Dbo5c5UWR9q0NehR/jNx1Kbfj2zS0yMbiSOUFWEEAKgJaS+l72eu0EmxI9UZdCm311tnvnla5Ldqw9AL2jj6lypOibmOMAIA1IyYRXVN/fm6oO4UbQ16NLLX4/wR87XAb9bqwafYdBZlQxgBAGqKZ57OT56ib9VfpJ1Bn/rzQ598L2bRcNPZR/S38W1VnBJzFWEEAKhJJ8WP0s2ZqzXucvs8sVbnxZX2krqr/0HtyPVWcULMRYQRAKBmfS66SP/U8A2lvZS2Bb2ffPA67SXly9cd/Q9qF6tjo4QIIwBATdvzxNphkSXqDLqVDzegbfLTGnVjurv/QQ3mh6s8JeYKwggAUPNSXp1uSF+iUxMnqDPo/mQD2ja/Ub35rH4x8Nh+N6UFikUYAQBmhahFdEXyHC1PnqXtQa+G8iOSpAVek/46vkUbhl5gjSPMGGEEAJg1zExn1Z2sG+svU78b0u6gX2amRX6rnh15RZtG3672iJjlCCMAwKxzTPxQ/WPm64pbXNuCXplM7X6zfjP4W3Xmuqo9HmYxwggAMCstjrTpRw3f0PHRw7Ql2CGTKeUl9NP+h9WfH6z2eJilCCMAwKxV58V1bf2F+nrqfPUGfXLOaSg/ojWDT/N5I0wLYQQAmNXMTF9JHKcfNXxDcYspp5xeG31H745/VO3RMAsRRgCAOWHirbVv6qt1X1HEInp97N1qj4RZKFLtAQAAKJU6L64LkqfotMQJ8uVXexzMQoQRAGDOqfeS1R4BsxRvpQEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgNCMwsjM/sXM3jKz181srZk1lmowAACASpvpHaMnJB3vnDtR0juSfjzzkQAAAKpjRmHknHvcOZcLv/yjpCUzHwkAAKA6SvkZo+9JerSE1wMAAKioyFQHmNmTkhbu51s/cc49EB7zE0k5SfdOcp2VklZK0sEHHzytYQEAAMppyjByzp0/2ffN7DuSlkv6mnPOTXKdVZJWSdLSpUsPeBwAAEC1TBlGkzGzZZL+k6RznHNDpRkJAACgOmb6GaP/LSkt6Qkz22Rmt5RgJgAAgKqY0R0j59zhpRoEAACg2lj5GgAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACA0IzCyMz+m5m9bmabzOxxM1tcqsEAAAAqbaZ3jP7FOXeic+4kSQ9J+i8lmAkAAKAqZhRGzrnsXl+mJLmZjQMAAFA9kZlewMz+u6RvS+qTdN6MJwIAAKiSKe8YmdmTZvbGfn5dIUnOuZ845w6SdK+kH01ynZVmttHMNnZ3d5fufwEAAECJmHOleffLzD4n6WHn3PFTHbt06VK3cePGkvxcAACAQpjZy865pZMdM9On0o7Y68vLJb01k+sBAABU00w/Y/Q/zOwoSXlJf5P0w5mPBAAAUB0zCiPn3NWlGgQAAKDaWPkaAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIBQScLIzP6DmTkzay3F9QAAAKphxmFkZgdJukDSRzMfBwAAoHpKccfo/5H0HyW5ElwLAACgamYURmZ2uaRO59xrJZoHAACgaiJTHWBmT0pauJ9v/UTS/ynpwkJ+kJmtlLQy/HLUzN4odMh5oFVST7WHqCG8Hvvi9dgXr8e+eD32xeuxL16PfR011QHm3PTeATOzEyT9VtJQ+FtLJG2V9BXn3PYpzt3onFs6rR88B/F67IvXY1+8Hvvi9dgXr8e+eD32xeuxr0JejynvGB2Ic+7Pktr3+mEfSlrqnKNMAQDArMQ6RgAAAKFp3zH6NOfcIUUcvqpUP3eO4PXYF6/Hvng99sXrsS9ej33xeuyL12NfU74e0/6MEQAAwFzDW2kAAAChqocR24lMMLP/Zmavm9kmM3vczBZXe6ZqMrN/MbO3wtdkrZk1VnumajKza8zsTTPLm9m8fMLEzJaZ2dtm9p6Z/edqz1NtZnanmXWx9MnEDgxm9rSZbQ7/f/LP1Z6pmswsYWZ/MrPXwtfjv1Z7plpgZr6ZvWpmD012XFXDiO1E9vEvzrkTnXMnSXpI0n+p9kBV9oSk451zJ0p6R9KPqzxPtb0h6SpJz1Z7kGowM1/Sv0q6WNKxkq4zs2OrO1XV3S1pWbWHqBE5Sf/eOXeMpFMl/dM8//djVNJXnXNfkHSSpGVmdmqVZ6oF/yxp81QHVfuOEduJhJxz2b2+TGmevybOucedc7nwyz9qYp2secs5t9k593a156iir0h6zzn3vnNuTNIvJV1R5Zmqyjn3rKSd1Z6jFjjntjnnXgn/uV8Tf/h1VHeq6nETBsIvo+Gvef1nipktkXSppNunOrZqYcR2Ip9lZv/dzD6W9A/ijtHevifp0WoPgarqkPTxXl9v0Tz+gw8HZmaHSDpZ0ovVnaS6wreNNknqkvSEc25evx6S/pcmbsTkpzqwZI/r70+pthOZKyZ7PZxzDzjnfiLpJ2b2Y0k/kvR/VXTACpvq9QiP+YkmbpPfW8nZqqGQ12Mes/383rz+GzA+y8zqJa2W9O8+dRd+3nHOBZJOCj+fudbMjnfOzcvPo5nZckldzrmXzezcqY4vaxg5587f3++H24kcKuk1M5Mm3iZ5xcym3E5kNjvQ67Efv5D0sOZ4GE31epjZdyQtl/Q1Nw/WlSji34/5aIukg/b6es8WRIAkycyimoiie51za6o9T61wzu02s2c08Xm0eRlGks6QdLmZXSIpISljZj93zl2/v4Or8laac+7Pzrl259wh4cKQWyR9cS5H0VTM7Ii9vrxc0lvVmqUWmNkySf9J0uXOuaGpjsec95KkI8zsUDOLSbpW0oNVngk1wib+hn2HpM3Ouf9Z7Xmqzcza9jzJa2Z1ks7XPP4zxTn3Y+fckrA3rpX01IGiSKr+h6/xd//DzN4ws9c18RbjvH7cVNL/lpSW9ES4hMEt1R6omsxshZltkXSapIfNbEO1Z6qk8IP4P5K0QRMfrP21c+7N6k5VXWZ2n6QXJB1lZlvM7PvVnqmKzpB0g6Svhv+92BTeHZivFkl6Ovzz5CVNfMZo0kfU8XesfA0AABDijhEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgj/P/t1rEAAAAAwCB/6z3DKIoAgIkRAMDECABgAev1k112TUraAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ell=[]\n",
    "for k in range(len(dats)):\n",
    "    mean = dats[k][0]\n",
    "    width = 1/(2*dats[k][1])\n",
    "    height = 1/width\n",
    "    rot = 180*dats[k][2]/np.pi\n",
    "    ell.append(mpl.patches.Ellipse(xy=mean, width=width, height=height, angle = rot))\n",
    "    \n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "for j in range(len(ell)):\n",
    "    ell[j].set_alpha(0.5)\n",
    "    col=np.random.rand(3)\n",
    "    ell[j].set_color(col)\n",
    "    ax.add_artist(ell[j])\n",
    "    \n",
    "    x_values = [0,m[0]]\n",
    "    y_values = [0,m[1]]\n",
    "    plt.plot(x_values, y_values)\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlim(-4, 4)\n",
    "ax.set_ylim(-4, 4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we rotate the state\n",
    "phi=np.pi/4\n",
    "oo=np.array([[np.cos(phi),-np.sin(phi)],[np.sin(phi),np.cos(phi)]])\n",
    "\n",
    "m=m@oo\n",
    "\n",
    "#Consider covariance matrix defined by \n",
    "# gg=(d/2)+(1/(2*d))\n",
    "# hh=(d/2)+(1/(2*d))\n",
    "# xx=(1/2)*np.array([[gg-(hh*np.cos(thet)),-np.sin(thet)*hh],[-np.sin(thet)*hh[0],gg+(hh[0]*np.cos(thet))]])\n",
    "\n",
    "d=3\n",
    "## Negative correlation goes with positive thet in our convention\n",
    "thet=.8-(2*phi)\n",
    "\n",
    "dats=[]\n",
    "dats.append([m,d,thet])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAJDCAYAAADuENvfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZRcdb3v/c+3uquqM89jDwkIRGbBOKIoHlQUDjgriqACEUhIGD3X433us+5z1133rMsYwmQYRAERjogoiAwKggcUAZknEUx6SNLd6XR6rml/nz9qRxLS6a7qru7qqnq/1spadLpS/aNStfc7v733b5u7CwAAAFKk2AMAAACYKAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgVLIzMrMrM/mpm9xTqOQEAAMZTIWeM1kh6pYDPBwAAMK4KEkZmVifpWEnXF+L5AAAAiqFQM0aXS/qepKBAzwcAADDuqkf7BGZ2nKRWd3/azD4+xONWSFohSVOmTHnvu9/97tH+aAAAgJw9/fTT7e4+b6jH2GjvlWZm/0fSNyWlJdVImi7pF+5+0p7+zPLly/2pp54a1c8FAADIh5k97e7Lh3rMqA+lufv33b3O3ZdK+pqk3w8VRQAAABMV6xgBAACERn2O0c7c/RFJjxTyOQEAAMYLM0YAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIRGHUZmVmNmT5rZc2b2kpn9z0IMDAAAYLxVF+A5EpI+4e49ZhaV9Eczu8/d/1SA5wYAABg3ow4jd3dJPeGX0fCXj/Z5AQAAxltBzjEysyoze1ZSq6QH3f3PhXheAACA8VSQMHL3jLu/R1KdpPeb2UHvfIyZrTCzp8zsqba2tkL8WAAAgIIq6FVp7t4p6RFJxwzyvfXuvtzdl8+bN6+QPxYAAKAgCnFV2jwzmxn+9yRJR0t6dbTPCwAAMN4KcVXaIkk/NrMqZUPrDne/pwDPCwAAMK4KcVXa85IOK8BYAAAAioqVrwEAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAGGPuriDwYg8DOSjELUEAAMAeNDcFuv2WpPr7pB/8fzXFHg6GQRgBADAGgsD1yENp/eKOlLZtdX3ys+xySwF/SwAAFFh/v+v2W5L68+NpzZsXUcSkY0+IFntYyAFhBABAAbW1BrruqoRaml0NSyLauCHQV0+Kad58TustBYQRAAAF8spLaV1/dVJmUn1DRM1NgQ5fXq0jj2J3Wyr4mwIAYJTcXY89ktbtN6c0e45p6jTTtg7XjBmmE0+JKRKxYg8ROSKMAAAYhSBw/eZXKd17d1qLa02xmCmRcHV3u77332s0dSpRVEoIIwAARiiddt1xS0qP/SGtunpTdbUpCFwtzYG+9s2YluzFeUWlhjACAGAEBgZcN12X1PN/zahhiSkSMbm7mhpdH/5olPOKShR/awAA5Kmry3XdlQm99WZGDUsiMsseLtu8ybXPfhF97ZtRzisqUYQRAAB52Noe6KrLEtra7qqrfzuKtrYHmj7DdOqZccViRFGpIowAAMhRc2OgKy9LKJV0La59+/yh7m5XOi2duSau6dOJolLGWWEAAOTgjdczuvj/DCgIfJfFGhMJ17YO13fPjmvRYnarpY4ZIwAAhvHsM2ndcHVSM2aapu00I5TJuFqaXSefGtWy/auKOEIUCmEEAMAe7Fi48Wc/SWn+AtOkybbL95o2uj59bLU+9BF2p+WCv0kAAAbh7rrv1yn9+q60Fi02xeO7njvU3OQ6bHmVjv9C9J8nYKP0EUYAALxDELh++fOUHrwvpbr6iKqrdw2fLZsD1dabTvpOTFVVRFE5IYwAANhJJuO646cp/eF3aTUsiewWPts6XPEa03dXxTVpElFUbjh9HgCAUCrluvnGpB77fVpLltpuUdTX6+rvd521Jq5Zs9mFliNmjAAAUPay+5uuS+q5ZzKqD2/x8c7vt7W6zjwnpvolRFG5IowAABWvv991/dUJvfZKoIYlttvJ1Ol09rL8r58S1cGHsussZ/ztAgAqWk+Pa/26hN58M1Bd/e5RFASuxo2uzx5frY9+nN1mueNvGABQsbq2u666bECbNrnq6naPIvdsFH3kY9U67nNcll8JCCMAQEXq2BroyksT6tjqqq0d/Jyh5ibXQYdE9NWTorudc4TyRBgBACpO65ZA6y5OqK/P93h/s82bAtXVR/Tt78YVjRJFlYIwAgBUlJbmQOsuSSidcs1fMHgUtbcFmjbddMZq1iqqNIQRAKBiNG4IdMUlAzKT5s4bPIq2d7rMpJXnxjVjJlFUaViIAQBQEd58I6NL/2NA1VWmOXMG3/319rh6e10rz63RgoXsIisRM0YAgLL36stpXbM2qalTTdNnDD4LlBhwbW13rTwvpiV7EUWVijACAJS1559N64frkpo9xzR16uBRlEq5Wlpc3zo9qgMOYtdYyfjbBwCUraf+nNaN1yY1b4Fp8uTBoyiTcTU1Bvrcl2L64BHRcR4hJhrCCABQlp74r5Ruvj6lBQtNNXu4sszd1bTR9YlPRvXpY9klgjACAJShxx5J6dabUlq82BSv2XMUNW50Hba8Sl/4KqtaI4swAgCUDXfXI79L6/ZbUlpca4rH9xw7m1pc79onom+eGlN1NVGELMIIAFAW3F0P3Z/WnbenVFtnisX2HDutWwLNmWs6fVVcNXuYUUJlIowAACXP3fXbe9K6+86U6uptyFt4bOtwxeKmlefGNW0aUYRdsVADAKCkubvu+WVKd9+ZVH3D0FHU3eVKJF2rzo1rzlx2gdgdM0YAgJIVBK6770zp/ntTqm+IDHmuUF+vq6vLtebCuGrriSIMjjACAJSkIHDd+bOUfv9gWg1LIqqq2nMUJQZc7W2uM8+J6V37Vo3jKFFqCCMAQMnJZFy335rSYw+nVd9gQ0bRjlWtv70iqoMOYbeHofEOAQCUlHTaddtPUnrisbQalpgikT1HUTrtamp0ffnEqD7wYVa1xvAIIwBAyUinXbf8KKknH8+ofpgoCoLsAo7HHFetT3yK3R1ywzsFAFASUinXzTck9dSTw0fRjlWtP/rxah3/BVa1Ru4IIwDAhJdKuX58fVLP/CWjhiU2ZOjsiKJDD6/SV0+KDhlQwDsRRgCACS2ZdP1ofVLP/3X4KJKklmbXPvtFdMpp3OoD+SOMAAATViLhuvHapF58PqP6huGjaMumQAsXmU5fya0+MDKEEQBgQhoYcF1/dUKvvhzkFEXtrYGmTjOddU6Npk4lijAyhBEAYMLp789G0WuvBKqrHz6KtnW4ItXSqvPjmjmLKMLIsSY6AGBC6etz/XBdQq+/mlsUdW13JZOu1efXaP4CdmsYHd5BAIAJo7c3G0V//1ug2rrho6i3x9Xd7Vp1Hvc/Q2GM+l1kZvVm9rCZvWJmL5nZmkIMDABQWXp7XdesTeitv+cWRf39rvZ21xmrY9p7H+5/hsIoxDlGaUnnu/szZjZN0tNm9qC7v1yA5wYAVIDu7mwUNW0MtLh2+ChKJFybN7lWrIzpgIM4XRaFM+oZI3ff5O7PhP/dLekVSbWjfV4AQGXo6nJdeemAmhsD1dZFho2iVMrV0uz6+ilRHf4+ogiFVdB3lJktlXSYpD8X8nkBAOWpa3s2ilq3uBbXDv9v9UzG1dQY6IQvxvTRjxNFKLyCnalmZlMl3SnpHHfvGuT7K8zsKTN7qq2trVA/FgBQorZ3uq64eECtra5Fi4ffHQWBa+MG1yc/E9Uxx1Vz/zOMiYKEkZlFlY2iW939F4M9xt3Xu/tyd18+b968QvxYAECJ6tzmWnvRgLa2uxYtGn5X5O5q3OA64shqff7L3BQWY2fU85CWfXfeIOkVd7909EMCAJSzjq2Brrg4oa7troW5RtFG16HvrdKJJ3NTWIytQswYHSHpm5I+YWbPhr8+W4DnBQCUma3tgdZelFB3l2vBwtx2QS3Nrn2XRfSt07kpLMbeqGeM3P2PkninAgCG1N6WjaL+Ps95herNmwItqs3eFDYeZ1eDsccp/QCAMde6JdAVFyU0MOCaNz+3KGrbEmjGTNNZa2o0ZQpRhPFBGAEAxtSWzdmZolQy9yjaujVQNG5adV5cM2YSRRg/3FgGADBmNm8KdPn/TSidcs2dl9suZ3unK8hIZ58fz/nPAIXCOw4AMCayUTSgTNo1Z25uu5uebldvr2vVeTU5LfgIFBrvOgBAwe04fBZklHMU9fe5tm1znbUmpqV7s3tCcfDOAwAU1JbN2cNn+cwUJRKu1i2uU8+MadkBnP6K4iGMAAAF07olO1OUTuUeRTtuCnvSd6I67L1EEYqLMAIAFETrluxMUSqZ+4nW6bSrqdH1xa9F9eGPRsd4hMDwCCMAwKi1teYfRUGQvdXHMcdV6+hPM1OEiYEwAgCMSltroLV5RtGOm8IeeVS1jv8CN4XFxEEYAQBGrK01e05RYiC/KNq4wfW+D1bpK9/gprCYWAgjAMCI7Lj32UC/a26OK1q7Zw+fHXJYlU76DjeFxcRDGAEA8tbelj18NtCf+20+3F1NTa5l+0f07RUxRaNEESYewggAkJet7dmZov6+3KNIklpaXEv3iui0s+KKx4kiTEyEEQAgZ1vbszNFfb2ueQty34Vsagm0aJHpjNVxTZ5MFGHiIowAADnp2JqdKerrc83PI4paNweaPce08twaTZ1KFGFiI4wAAMPaEUW9PfkdPmtvDTR5imnVeXFNn0EUYeIjjAAAQ9rWEeiKixLq6c5vpmhre6DqmGn1hXHNnsPuBqWBdyoAYI+2dWRnirq78ouibR0ul7T6gnheM0xAsfFuBQAM6p9RtN01f2Huu4uu7a5E0rX6/BotWsxuBqWFdywAYDed21xrL0qoK88o6u529fS4zj4vrvol7GJQenjXAgB20bnNdcXFA+rqdC3II4p6e12d21xnnRPX3vtUjeEIgbFDGAEA/ml7ZzaKOre5FizKfRfR3+/a2u46Y3VMy/YnilC6CCMAgKRsFK29KIyiPGaKEgOuLZtcp54R00GHVI/hCIGxRxgBALS907XukgFt68gvipJJV0uL6+TTojr8fUQRSh9hBAAVrmt7Noq2trsW5nH4LJVyNTe5vvqNqD70kegYjhAYP4QRAFSwru3Zc4ryjaJ02tXUGOhzX4rq40czU4TyQRgBQIXq6nJdeemA2vOMokzG1bgx0DHHRfXpY6tlxq0+UD4IIwCoQF1drisvGVDrFteiPKIoCFyNG1yf+GRUx38hShSh7BBGAFBhdomiPFamdndt3OD68JHV+uLXoopEiCKUH8IIACpId7frqstGHkXvfX+VTjw5qqoqogjliTACgArR3Z09p2hzS/5R1NToOuiQiE4+NabqaqII5YswAoAK0N3tuvqyhDZvci2uzW/T39zk2me/iL5zRlyxGFGE8kYYAUCZ6+nJRlFLc6DFed7tvqU5UP2SiE5fGVdNDVGE8kcYAUAZ6+lxXX15Qs3NQd4zRZs3BZq/wHTmmrimTCGKUBkIIwAoUzuiqKkxUG2eUdS6JdCMmaZV59Vo2jSiCJWDMAKAMtTb67rm8oSaNuYfRe1tgSZNNp19flwzZhJFqCyEEQCUmd5e1zVrE2rcGKi2Lr/NfMfWQJEq6ezz45ozl10EKg/vegAoI729rmvXJrTxH4EW1+Y329O5zZXOSKsvqNGCheweUJl45wNAmejrc/3wioT+8VY2ivK5XUd3l6t/wLX6/Jq8Z5mAcsK7HwDKQF9fdqborTcD1dblF0U9Pa6uLteqc+Nashe7BVQ2PgEAUOL6+lw/XJfQW2/lH0V9fa6Ora4z18S0z35VYzhKoDQQRgBQwvr7s1H05huBavM8fDbQ72pvda1YFdP+B1aP4SiB0kEYAUCJ2hFFf38j/5miRMK1ebPrWytiOvQwogjYgTACgBK0I4re+FugujyjKJl0tTS7vvGtqN73QaII2BlhBAAlZmDAtf7KhP4+gihKp13NTa4vnxjVRz4WHcNRAqWJMAKAEjIwkJ0p+ttr+R8+y2RcjRsDHf+Fan3iU8wUAYMhjACgRAwMuK67KqG/vZZRXX1+URQEro0bXJ/6bFSf+ddoXn8WqCSEEQCUgIEB13VXJ/TaKxnV1UdGFEUf+5dqfe5LRBEwFMIIACa4RMJ1/dUJvfZykHcUubsaN7g+dES1vvL1qCIRoggYCmEEABPYjih65aUg78Nn7q7Gja7DllfpxFOiqqoiioDhEEYAMEElEq4brkno5RcD1TfkH0VNja79D4zo5NNiikaJIiAXhBEATECJhOuGa5MjiiJJaml2vWvfiE49M654nCgCckUYAcAEk0y6brw2qZeez//qM0na1BKoti6iFavimjSJKALyQRgBwASyI4pefD4zopmiLZsCzZlrOvOcuKZMIYqAfBFGADBBpFKum9Yn9cJzI4uiti2Bps0wrTovrunTiSJgJAgjAJgA0mnXzTck9ewzI4uire2BYjWms8+Pa9ZsNu3ASPHpAYAiS6ddt/woqaeezKhhSf5RtK3DJZNWXxDX3Hls1oHR4BMEAEWUybh+dnNKTz4+spmi7Z2uZMq15oIaLVzEJh0YrYJ8iszsRjNrNbMXC/F8AFAJgsB1+60p/dejadUvsbxXpe7ucvX1uVafX6PaeqIIKIRCfZJuknRMgZ4LAMpeELh+fltKjz2cVsMIoqi3x9XZ6Vp5blxL9yaKgEIpyKfJ3R+V1FGI5wKAcufu+uXPU3r4obTqG/KPov4+V3u768w1Me27rGqMRglUJv6ZAQDjyN31q1+k9MBvUqpvsLzvXzYw4NqyxXX6ypgOPLh6jEYJVK5xCyMzW2FmT5nZU21tbeP1YwFgwnB33furlO77dUr1DZG8oyiRcG1qcX3r9KgOey9RBIyFcQsjd1/v7svdffm8efPG68cCwITg7vrtPWndc1da9Q0RVVfnF0XJpKul2XXyqVF94MPRMRolAA6lAcAYc3c9dH9ad9+ZUl295R1FqZSrucl14slRffijRBEwlgp1uf5tkp6QtMzMmszs1EI8LwCUg0ceSuvO27NRFI3mH0VNja4vnRjVkUdx+AwYawX5lLn7iYV4HgAoN48+nNLtt6ZUW5d/FKXT2Sg64YtR/cunqvNe/BFA/jiUBgBj5PHHUvrpj1NaXGuKxfKLmkzG1bgx0HGfq9YxxxFFwHghjABgDPz58ZRuvjEbRfH4SKLI9eljozr2hChRBIwjwggACuzpJ9P68XUpLVyUfxQFgWvjBtdRR1frhC8SRcB4I4wAoID++nRaN16b1IKFppqakUXRkUdV60snRvNeERvA6BFGAFAgLzyX1vVXJzVvvqlmUn5R4549fPahI6r11ZOIIqBYCCMAKICXX0zrh+uSmjvXNGly/lG0cYNr+fur9PVvRfNeERtA4RBGADBKr76c1jWXJzV7tmnylJHNFB16eJW+eWos78UfARQWYQQAo/C31zK6+rKkZs4yTZmafxQ1NboOPDiib6+I5b3OEYDCI4wAYITefCOjKy9NaPoM09RpI4iiJte+yyL6zhnxvNc5AjA2CCMAGIF/vBnoiksSmjrVNG16/lHT0uzae++IVqyK5331GoCxQxgBQJ42/iPQFRcPaPIk0/QZI4miQHUNEX13dVyT8rx6DcDYIowAIA9NGwOtvXhA8bhpxsz8o2ZTS6CFi0xnrolrSp4nagMYe4QRAOSopTk7UxStNs2clX/UbNkcaO4808pzazQtz3OSAIwPwggAcrB5U6C1Fw3ITJo1O/+oad0SaPoM06rzakZ0+A3A+CCMAGAYrVsCrb0oIQ+k2XPy32y2tQaaMsW0+oL4iGaaAIwfwggAhtDeFuiKixJKp1xz5ua/yWxvCxSPm1ZfGB9RVAEYX3xKAWAPOrZmo2hgwDV3Xv6by63tgaqqpdUXxkf05wGMPz6pADCIzm2uKy5OqLfXNW9+/pvKbR0umbTmwhotWMimFigVfFoB4B26truuvHRAXdtd8xfsvpns88SQf75zmyuVdq25oEaLFrOZBUoJn1gA2El3t+uqywbU3u6DzvS8mt6kn/Y/ofage9A/37XdNZDIRlFtPZtYoNTwqQWAUG+v65q1CW3e5Fq0aPfNY3vQrUeTr2peZLpm25Tdvt/d5ertda0+P66GpWxegVLEJxcAJPX3u354RUJNGwMtrt1905jwlO5PvKAai+qT8QMVsV0f09Pt6upyrTovrr3eVTVewwZQYIQRgIo3MOBaf2VCb74ZaHHt7usMubt+l3hZPT6gT8UO1mSL7/L93l5X5zbXynNi2mc/oggoZYQRgIqWSLhuvDahv70WqK7OZLZ7GD2T/oc2BO36cHRfLayascv3+vpcW9tdZ6yJadkB1eM1bABjhDACULFSKddN1yX10guB6uoHj6LGzFY9mXpT+1Yt0EHVdbt8r7/f1d7qWrEqpgMPJoqAckAYAahI6bTr5huSev6ZjOobBo+i7mBADyVe0myboo/F9t/lMQMDrtYtrtPOiunQw4gioFwQRgAqTibj+ulNKT31ZEb1SwaPoowHeiD5ggIF+nT8EEXt7XOHEgOuLZtc314R02HLiSKgnBBGACpKELhuvzmlP/1Xeo8zRZL0x9Trag26dFTsAM2MTP7n7ycSrk2bXCefFtX7PkgUAeWGMAJQMYLA9fPbUnr0kbTql5gikcGj6NX0Jr2cbtZ7qpdo7+r5//z9ZNLV0uw66dtRffCI6HgNG8A4IowAVAR31913pvTwQ2k1DBFFOxZxXByZpQ9E9/7n7yeTruZm14knR3XEkUQRUK4IIwBlz93167tSuv/elOobTFVVg0fRnhZxTKVczU2ur5wY1ZFHcfgMKGeEEYCy5u767T1p/eZXKdU3RPYYRXtaxDGddjU1uj7/laiO+mT1Hs9JAlAeCCMAZcvd9bsH0rr7zpTq6iOqrt5z1Ay2iGM67WrcGOj4L1TrU58hioBKQBgBKFuPPpzWz29Lqa7eFI3uOWoGW8Qxk8lG0bEnRPWZf40SRUCFIIwAlKXHH0vptp+kVFs3dBQNtohjJuPauCHQp4+N6rjPEUVAJSGMAJSdJ59I6ZYbU1pca4rF9hw1gy3iGASujRtcRx8T1QlfJIqASkMYASgrz/wlrZuuS2nBIlM8PnTUvHMRxx1R9PGjq/WFr0T3eEk/gPJFGAEoG88/m9YN1yQ1f4GppmboqHnnIo5B4Grc4Prox6r15ROJIqBSEUYAysLLL6a1/sqk5s43TZo0dNS8cxFH9+xM0Qc/Uq2vfjO6x0v6AZQ/wghAyXvtlYyuuTyp2XNMkycPHTXvXMTRZNq4wfX+D1Xp66cQRUClI4wAlLQ3Xs/oqssSmjnLNGXK0FHzzkUcJymmjRtch7+vSid9OzbkOkcAKgNhBKBk/ePNQFdemtD06aap04aPmp0XcVwQma7Gja5DDqvSyafGhrykH0DlIIwAlKTGDYGuuHhAkyebpk0fPmp2XsTxwKpaNTW6Djw4om+viA15ST+AykIYASg5zY2B1l40oHjcNGPm8FGz8yKOR0bfrZZmadn+EX3njPiwl/QDqCyEEYCSsnlToLUXD6i62jRz1vBRs/Mijp+KHazWlojetW9Ep50VH/aSfgCVhzACUDLaWgNdcXFCcmnW7NyiZudFHPs2TdLSpRGtWBUf9pJ+AJWJMAJQErZ1BFp3SULJhGvO3Nw2Xa/ttIhjzZa5qmuI6Iw18WEv6QdQuQgjABNe13bXuksS6u12zZuf22arPejWH8JFHBval2rhYtOZa+LDXtIPoLIRRgAmtJ4e19WXJ9Sx1TV/YW6brJ0XcTy08wAtmF+llefWaFoOl/QDqGyEEYAJq7/ftX5dQi3NgRYtzm1ztfMijsu7D9SiWXGtOq9G03O4pB8ACCMAE1Ii4br+moTeejPQ4trco2bHIo4HD+yjvabN1OoLanK6pB8AJMIIwASUSrl+fF1Sr74UqLbOZJZb2OxYxLE+tUCHxOq0+oJ4Tpf0A8AOhBGACSWTcd16U1LPPpNRfUPuUbRjEcfpmSk6IvJurbmwRrPnsIkDkB+2GgAmjCBw3fHTlJ58PL8o2rGIY9oDfSw4WOdfOFlz57F5A5A/thwAJgR31y9/ntIffpdW/RJTJJL7IbAdizh+KHOA/p/vzcr5kn4AeCe2HgCKzt11769SeuA3KTXkGUU7FnHcP9OgSy6s0/wFbNYAjFx1sQcAAL9/IK177kqrviGiqqrco6g96NYfEq9qgc/U9efuTxQBGDW2IgCK6rFHUvrP21KqqzdVV+ceRQlP6Td9zyumqG4643DVLubfeQBGjzACUDRPPpHST3+c0uJaUzSaexS5u+7reUn9ltAVXz1M++89aQxHCaCSFCSMzOwYM3vNzN4ws/9WiOcEUN6efSatm65LaeEiUzye31pDf9z+ljZVbdW5R+2vTx4+Z4xGCKASjTqMzKxK0lWSPiPpAEknmtkBo31eAOXrlZfSuv7qpObNN9XU5BdFL3W068XoW/rUskVa9emlYzNAABWrEDNG75f0hru/6e5JST+TdEIBnhdAGXrj9YyuuTypWbNNkyfnF0UbO/r0p/jL2nvuVF3+jUNyXucIAHJViDCqldS409dN4e8BwC42vBXoyssSmj7DNHVqflHTtjWtJ2IvKpFrtXMAABf5SURBVBp3Xf+t92pyjJOtARReIbYsg23dfLcHma2QtEKSGhoaCvBjAZSS5qZA6y4Z0ORJpml53um+Y2uglDLaqyGm0z++n/aeN3WMRgmg0hUijJok1e/0dZ2klnc+yN3XS1ovScuXL98tnACUry2bA11x8YCqqpX3ne47OgK5pO9dOE2Laz/A4TMAY6oQh9L+ImlfM9vLzGKSvibpVwV4XgBlYGt7oHUXJ+SBNHt2fpucjo5AHkhrLqxRbV2EKAIw5kY9Y+TuaTNbJel+SVWSbnT3l0Y9MgAlr3Oba90lCQ0MeN73L9vW4QoC6ZwwigBgPBTk7EV3/42k3xTiuQCUh+5u15WXDqhru2vBwvyjKBN4NorqiSIA44ctDoCC6+11XbM2obbWEUZRxrXmAqIIwPhjqwOgoAYGXNddlVDjhkCLa/PbxHRuc6XTrjUX1qiugc0TgPHHlgdAwSSTrhuvTeiN1zOqrcvvROnOba5UiigCUFxsfQAURDrt+skNSb30QqC6+vyuIOvc5koks1FUv4TNEoDiYQsEYNQyGddtP0npmb9kVN9gI4qic4giABMAWyEAoxIErp//LKXHH0vnHUXbO12JRDaKGpayOQJQfGyJAIyYu+tXv0jpkYeyURSJ5BdFAwPZw2dEEYCJgq0RgBFxd/32nrR+e09K9Q2mqqr8oqi/PxtFS/ZiMwRg4mCLBGBEHnkorbvvTKmuPpJXFHVtd/UPEEUAJia2SgDy9vhjKd1+a0q1daZoNL8o6uvLLt64dG82PwAmHrZMAPLy9JNp3XJjSotrTbFYflHU2+dafUGcKAIwYbF1ApCzF55L68Zrk1qw0BSP5xlFva7V58e117uqxnCEADA6hBGAnLz2clo/XJfU3Hmmmkm5R1F3VxhFF8S19z5EEYCJjTACMKw338jo6rVJzZplmjwlvyjq6XGdfT5RBKA0EEYAhtS4IdC6SxOaOtU0dVp+UdTdnY2id+1LFAEoDYQRgD1qaQ50xSUDqombps8gigCUP8IIwKDaWgOtuyQhM2nmrDyiqNvV1eVadV5c++xHFAEoLYQRgN10bA10xUUJpZOuOXNy30x0d7u6tmdnivZdRhQBKD2EEYBddG13XXlpQn19rrnzc99E9HS7ujqzM0VEEYBSRRgB+KeeHtdVlw1oW4dr/oL8oqiz07XyvLj2ezdRBKB0EUYAJEl9fa5rr0ho0ybXwkV5RFGPa/t216pzY1q2P1EEoLQRRgCUSLiuuyqhDW8FWrw49xOte3pcndtcZ50T07IDqsdwhAAwPggjoMKlUq4f/TCp118NVFtnMsstjHp7XJ0drpXnxvRuoghAmSCMgAqWTrtuuTGpF57NqL4hvyjq6HCdRRQBKDOEEVChgsB1xy0p/eVPGdUvySOKesMoOiem/Q8kigCUF8IIqEDurl/ckdKjj6Tzj6Kt2Sg64CCiCED5IYyACuPuuueXKT3025Qalpgikfyi6Mw1RBGA8kUYARXE3fXgb9P6za/Sqm+IqKoqtyjqC6PojNUxHXgwUQSgfBFGQAV59OG0fnF7SrV1purq3KOovT0bRQcdQhQBKG+EEVAh/vx4Sj+7ORtF0WiOUdTnamsjigBUDsIIqADP/TWtH1+f0qJFplgs9yhqD6Po4EOJIgCVgTACytzrr2Z03VVJzZtvitfkFkX9YRStWBXTIe8higBUDsIIKGMb/xHo6rUJzZxlmjw59yhqa3WdvjKmQw8jigBUFsIIKFObWgKtu3RAk2pM06blF0WnrYzpPYcTRQAqD2EElKGt7YHWXZKQSZo5K8co6n87ig57L1EEoDIRRkCZ6druuvLShBIDrjlzc/uI9/e7Wre4TjuLKAJQ2QgjoIz09bmuWZvQtg7X/AV5RtGZMR22nCgCUNkII6BMJBKu665KqKkx0KLFOUZRn6t1czaKDn8fUQQAbAmBMpBOu35yfVJ/ey1QXX3uK1q3tblWnM2J1gCwA1tDoMQFgev2W1L669MZNSwxmQ0fRjvfEJbFGwHgbWwRgRLm7rr7zpQeeyStJUtzi6KeHte2jmwUcUNYANgVW0WgRLm7HrgvrfvvTalhSUSRSA5R1O3q7HStPDem/Q/k4w8A78SWEShRjz+W1l3/mVJdfURVVcNHUXeXq2u7a9V5cS3bv2ocRggApYcwAkrQX59K69YfpbR4sSkazS2KurtdZ18Q177LiCIA2BPCCCgxr76c1vXXJDV/oSkeHz6Kura7entdqy+I6137EkUAMBTCCCghb/09o2vWJjV7jmnSpOGjaHunq3/AtebCuPZ6F1EEAMMhjIAS0dwU6KrLEpo61TR16vBR1LnNlUi6zrmwRkv2Yi1XAMgFYQSUgLbWQFdemlBVtTR9xvBRtK3DlU67zv1ejeqXEEUAkCu2mMAEt70ze1PYdNI1e/bwH9mOjkCZwHUOUQQAeWOrCUxgPT2uqy4bUHeXa+784T+uW7cGcpfO+V6N6hr4eANAvthyAhPUwED2prCbN7kWLMwhitoDmUnn/luNauv4aAPASHCOETABpVKum9Yn9fe/5XZT2Pa2QNVR05oL41q4iCgCgJEijIAJJpNx3fbjlF54LqP6huHvf9beGigaz0ZRLjNLAIA9I4yACcTdddd/pvTEf6XVsGT4KGrbEmjSZNPqC+Oal8M5SACAoRFGwATh7nrgN2n97v6U6huGvyls65ZAU6aaVl8Q19x5RBEAFAJhBEwQT/wxrbt+nlJd3fA3hW3dHGja9OxM0ew5RBEAFAphBEwAzz6T1i053hR2y+ZAM2ZmZ4pm5bCuEQAgd4QRUGSvv5rR9VcnNW/+8DeF3bwp0Ow5prPPr9HMWcNfrQYAyA9hBBRR44ZAV69NaOYs0+TJQ4dOS3Og+QtMq86r0YyZRBEAjIVRzcOb2ZfN7CUzC8xseaEGBVSCLZsDrbt0QJNqTNOm7Tl03F3NzYEW10Z09gVEEQCMpdGeoPCipC9IerQAYwEqRue27P3PPNCQh8TcXc1NroaGiFaeF9f06UQRAIylUR1Kc/dXJA271gqAt/X0uK6+fEC93a4FQ6xSvSOK9t4nou+eHR/2UBsAYPS4pAUYR4nETvc/GyaKGje69nt3RGesJooAYLwMO2NkZg9JWjjIt37g7nfn+oPMbIWkFZLU0NCQ8wCBcpHJuG6+Mak3Xg9U3zD04bPGja4DD47oO2fEVVNDFAHAeBk2jNz96EL8IHdfL2m9JC1fvtwL8ZxAqXB33fmzlJ55MqOGpXu+1Ye7q3GD65DDq/St02PDXr4PACgsLtcHxsED96X1+wdTalgS2WMUBUF2pui976/SyafGhl3oEQBQeKO9XP/zZtYk6UOS7jWz+wszLKB8/PnxlO66I6W6+j3f6iMIsjNFH/hQtU45jSgCgGIZ7VVpd0m6q0BjAcrOKy+l9ZMbUlo0xK0+MpnsTNERR1brxJOjw94nDQAwdjiUBoyRDW8F+uG6pGbPsT2eQJ1Ou5oaAx11dFRf/BpRBADFRhgBY6B1S6CrLh/QpEmmqVMHj51UytXU6Prs8VH96+ejrAcGABMAYQQUWNd211WXhatazxk8dpLJ7OKNn/9KVJ/6TDVRBAATBGEEFFB/v+vadQlt73QtWjz4tQ2JAVdLi+vEk6M68iiiCAAmEsIIKJBUynXT+qQ2/iNQfcPgUdTf52ptdX17RVQf+HB0nEcIABgOYQQUQBC47rg1pReey6hhyeAzQL09rq1bXaevjOmw9/LRA4CJiK0zMErurt/8KqXHHklryR5Wte7ucnV1uVaeG9MBB/GxA4CJii00MEp//ENa9/wyrfoGUySyexRt73T19bnOPj+ufZdVFWGEAIBcEUbAKDz317Ru+3FKi2tN1dW7R9G2Dlc67TrnezVauveoFpoHAIwDwggYoTffyOj6q5Oat8AGvdnr1vZAMuncf6tRbT1RBAClgDACRmBTS6CrLkto+nTT5Mm7R1HrlkCTJpnOviCuBQuJIgAoFYQRkKeu7a6rL0+oqkqaPmP3KNqyOdD0Gaazz49rzlyiCABKCVttIA+JhOv6qxPq7vJBo6elOdDsOaZzvldDFAFACWLGCMhRELh+dnNSf38jUH3DrjNF7tlbfNQviejMNXFNm8Zq1gBQiggjIEf335vWE3/M7LZWkXv2ZrD77BfRilXxQc85AgCUBsIIyMHTT6Z1952p3dYqCgJX40bXoYdX6ZTTYqqpIYoAoJQRRsAw3nwjox+tT2rhol3XKkqnXY0bAx15VFRf+UZ00HWMAAClhTAChtDWGujaKxKaMcM0adLb4ZNIuFqaXcd9PqrP/mt00BWvAQClhzAC9qCnx3XN2oSCYNfL8vv6XG1bXN/4VlQf+Vj1oPdGAwCUJsIIGEQq5bppfULtba7aurcvu+/ucm3vdK04O6b3HM7HBwDKDVt24B3cXXfentLLLwZqWPL2bFBHR6BkUlrzvbj22Y+bwQJAOSKMgHd45KG0HnkopSVLI/88TNa6JVAsbrrg3+O7zCABAMoLYQTs5OUX0/rPn6ZUWxdRJGJyd21qcc2Za1p5Lrf4AIByRxgBoY6tgW64Nqk580yxmIULNwbaa+8qrTib1awBoBIQRoDCk62vS8oDaepUUyaTXc360MOrWbgRACoIYQRIuu/XKb3xevZk64GB7OGzoz5ZrS9+lYUbAaCSEEaoeC+9kNZ9v06rrt60vVPq6XZ96/SoPvBh1igCgEpDGKGidW5z/eiHSc2da2prdU2eYjr/32u0dG9OsgaASkQYoaLd9+uUenuzq1nvtyyib383rhkzmSUCgEpFGKFidW5zPfxQStVVps8eX61PHxtVNEoUAUAlI4xQsaZOk045LaZ9l1Vp/gIOnQEACCNUsOpq0xFHRos9DADABMI/kwEAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABAijAAAAEKEEQAAQIgwAgAACI0qjMzsIjN71cyeN7O7zGxmoQYGAAAw3kY7Y/SgpIPc/RBJr0v6/uiHBAAAUByjCiN3f8Dd0+GXf5JUN/ohAQAAFEchzzH6jqT7Cvh8AAAA46p6uAeY2UOSFg7yrR+4+93hY34gKS3p1iGeZ4WkFZLU0NAwosECAACMpWHDyN2PHur7ZnaKpOMk/Yu7+xDPs17Seklavnz5Hh8HAABQLMOG0VDM7BhJ/ybpY+7eV5ghAQAAFMdozzG6UtI0SQ+a2bNmdm0BxgQAAFAUo5oxcvd9CjUQAACAYmPlawAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQoQRAABAiDACAAAIEUYAAAAhwggAACBEGAEAAIQIIwAAgBBhBAAAECKMAAAAQqMKIzP7X2b2vJk9a2YPmNniQg0MAABgvI12xugidz/E3d8j6R5J/6MAYwIAACiKUYWRu3ft9OUUST664QAAABRP9WifwMz+t6STJW2XdNSoRwQAAFAkw84YmdlDZvbiIL9OkCR3/4G710u6VdKqIZ5nhZk9ZWZPtbW1Fe7/AAAAoEDMvTBHv8xsiaR73f2g4R67fPlyf+qppwrycwEAAHJhZk+7+/KhHjPaq9L23enL4yW9OprnAwAAKKbRnmP0H2a2TFIgaYOkM0Y/JAAAgOIYVRi5+xcLNRAAAIBiY+VrAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABChBEAAECIMAIAAAgRRgAAACHCCAAAIEQYAQAAhAgjAACAEGEEAAAQIowAAABCBQkjM7vAzNzM5hbi+QAAAIph1GFkZvWSPilp4+iHAwAAUDyFmDG6TNL3JHkBngsAAKBoRhVGZna8pGZ3f65A4wEAACia6uEeYGYPSVo4yLd+IOnfJX0qlx9kZiskrQi/TJjZi7kOsgLMldRe7EFMILweu+L12BWvx654PXbF67ErXo9dLRvuAeY+siNgZnawpN9J6gt/q05Si6T3u/vmYf7sU+6+fEQ/uAzxeuyK12NXvB674vXYFa/Hrng9dsXrsatcXo9hZ4z2xN1fkDR/px/2D0nL3Z0yBQAAJYl1jAAAAEIjnjF6J3dfmsfD1xfq55YJXo9d8XrsitdjV7weu+L12BWvx654PXY17Osx4nOMAAAAyg2H0gAAAEJFDyNuJ5JlZv/LzJ43s2fN7AEzW1zsMRWTmV1kZq+Gr8ldZjaz2GMqJjP7spm9ZGaBmVXkFSZmdoyZvWZmb5jZfyv2eIrNzG40s1aWPsnegcHMHjazV8LPyZpij6mYzKzGzJ40s+fC1+N/FntME4GZVZnZX83snqEeV9Qw4nYiu7jI3Q9x9/dIukfS/yj2gIrsQUkHufshkl6X9P0ij6fYXpT0BUmPFnsgxWBmVZKukvQZSQdIOtHMDijuqIruJknHFHsQE0Ra0vnuvr+kD0paWeHvj4SkT7j7oZLeI+kYM/tgkcc0EayR9MpwDyr2jBG3Ewm5e9dOX05Rhb8m7v6Au6fDL/+k7DpZFcvdX3H314o9jiJ6v6Q33P1Nd09K+pmkE4o8pqJy90cldRR7HBOBu29y92fC/+5WdudXW9xRFY9n9YRfRsNfFb1PMbM6ScdKun64xxYtjLidyO7M7H+bWaOkb4gZo519R9J9xR4EiqpWUuNOXzepgnd82DMzWyrpMEl/Lu5Iiis8bPSspFZJD7p7Rb8eki5XdiImGO6BBbtcfzCFup1IuRjq9XD3u939B5J+YGbfl7RK0v87rgMcZ8O9HuFjfqDsNPmt4zm2Ysjl9ahgNsjvVfS/gLE7M5sq6U5J57xjFr7iuHtG0nvC8zPvMrOD3L0iz0czs+Mktbr702b28eEeP6Zh5O5HD/b74e1E9pL0nJlJ2cMkz5jZsLcTKWV7ej0G8VNJ96rMw2i418PMTpF0nKR/8QpYVyKP90clapJUv9PXO25BBEiSzCyqbBTd6u6/KPZ4Jgp37zSzR5Q9H60iw0jSEZKON7PPSqqRNN3MbnH3kwZ7cFEOpbn7C+4+392XhgtDNkk6vJyjaDhmtu9OXx4v6dVijWUiMLNjJP2bpOPdvW+4x6Ps/UXSvma2l5nFJH1N0q+KPCZMEJb9F/YNkl5x90uLPZ5iM7N5O67kNbNJko5WBe9T3P377l4X9sbXJP1+T1EkFf/ka7ztP8zsRTN7XtlDjBV9uamkKyVNk/RguITBtcUeUDGZ2efNrEnShyTda2b3F3tM4yk8EX+VpPuVPbH2Dnd/qbijKi4zu03SE5KWmVmTmZ1a7DEV0RGSvinpE+H24tlwdqBSLZL0cLg/+Yuy5xgNeYk63sbK1wAAACFmjAAAAEKEEQAAQIgwAgAACBFGAAAAIcIIAAAgRBgBAACECCMAAIAQYQQAABD6/wFG9O/RdGJmmQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ell=[]\n",
    "for k in range(len(dats)):\n",
    "    mean = dats[k][0]\n",
    "    width = 1/(2*dats[k][1])\n",
    "    height = 1/width\n",
    "    rot = 180*dats[k][2]/np.pi\n",
    "    ell.append(mpl.patches.Ellipse(xy=mean, width=width, height=height, angle = rot))\n",
    "    \n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "for j in range(len(ell)):\n",
    "    ell[j].set_alpha(0.5)\n",
    "    col=np.random.rand(3)\n",
    "    ell[j].set_color(col)\n",
    "    ax.add_artist(ell[j])\n",
    "    \n",
    "    x_values = [0,m[0]]\n",
    "    y_values = [0,m[1]]\n",
    "    plt.plot(x_values, y_values)\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlim(-4, 4)\n",
    "ax.set_ylim(-4, 4)\n",
    "plt.show()"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}