ZWO_capture_2_element_polaris_240124.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4ccebea4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import cv2\n",
    "import argparse\n",
    "import os\n",
    "import sys\n",
    "import time\n",
    "import zwoasi as asi\n",
    "from time import time\n",
    "\n",
    "from scipy import optimize\n",
    "import laserbeamsize as lbs\n",
    "import matplotlib.animation as animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "47e48c8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import pearsonr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6b046e5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.animation as animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9bb2eb39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_dict_to_h5(ddict,filename='star_images.h5', initialdir=\"C:\\\\Starshot\\\\Polaris\"):\n",
    "    #use group with attributes to mock up dictionary, save to h5 file 'h5f'\n",
    "    h5file=initialdir+'\\\\'+filename\n",
    "    keys=ddict.keys()\n",
    "    with h5py.File(h5file,'w') as h5f:\n",
    "        for key in keys:\n",
    "            h5f.create_dataset(key,data=np.array(ddict[key]))\n",
    "        h5f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f7b5ee4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_dict_from_h5(filename='star_images.h5', initialdir=\"C:\\\\Starshot\\\\Polaris\"):\n",
    "    #if there is a group called 'dictionsary', read its attributes as key/value, otherwise\n",
    "    #read directly from an h5 file with no groups, datasets convert to key/value dictionary\n",
    "    h5file=initialdir+'\\\\'+filename\n",
    "    with h5py.File (h5file,'r') as h5f:\n",
    "        h5keys=list(h5f.keys())\n",
    "        if 'dictionary' in h5keys:\n",
    "            grp=h5f.get('dictionary')\n",
    "            out_dict = {}\n",
    "            for key in grp.attrs.keys():\n",
    "                out_dict[key]= grp.attrs[key]\n",
    "        else:\n",
    "            out_dict={}\n",
    "            for key in h5keys:\n",
    "                out_dict[key]=h5f[key][:]\n",
    "        h5f.close()\n",
    "    return out_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "161e578c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def smooth(y, box_pts):\n",
    "    box = np.ones(box_pts)/box_pts\n",
    "    y_smooth = np.convolve(y, box, mode='same')\n",
    "    return y_smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "71b541a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def star_centroid(image):\n",
    "    #quick and dirty centroid, only viable for 1 star in field\n",
    "    gray_image=image.copy()\n",
    "    ret,thresh = cv2.threshold(gray_image,247,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "    # calculate moments of binary image\n",
    "    M = cv2.moments(thresh)\n",
    "    # calculate x,y coordinate of center\n",
    "    cX = (M[\"m10\"] / M[\"m00\"])\n",
    "    cY = (M[\"m01\"] / M[\"m00\"])\n",
    "    return cX,cY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b72b7948",
   "metadata": {},
   "outputs": [
    {
     "ename": "OSError",
     "evalue": "[WinError 126] The specified module could not be found",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mOSError\u001b[0m                                   Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-6c4b5d7ed8e7>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0menv_filename\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetenv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'ZWO_ASI_LIB'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0masi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'C:\\\\Users\\\\peter\\\\Downloads\\\\ASI_Windows_SDK_V3.24.0.0\\\\ASI SDK\\\\lib\\\\x64\\\\ASICamera2.dll'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\zwoasi\\__init__.py\u001b[0m in \u001b[0;36minit\u001b[1;34m(library_file)\u001b[0m\n\u001b[0;32m    819\u001b[0m             \u001b[1;32mraise\u001b[0m \u001b[0mZWO_Error\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Require filename of the ASI SDK library'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    820\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 821\u001b[1;33m     \u001b[0mzwolib\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcdll\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mLoadLibrary\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlibrary_file\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    822\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    823\u001b[0m     \u001b[0mzwolib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mASIGetNumOfConnectedCameras\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0margtypes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\ctypes\\__init__.py\u001b[0m in \u001b[0;36mLoadLibrary\u001b[1;34m(self, name)\u001b[0m\n\u001b[0;32m    440\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    441\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mLoadLibrary\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 442\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_dlltype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    443\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    444\u001b[0m \u001b[0mcdll\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mLibraryLoader\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mCDLL\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\ctypes\\__init__.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, name, mode, handle, use_errno, use_last_error)\u001b[0m\n\u001b[0;32m    362\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    363\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mhandle\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 364\u001b[1;33m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_handle\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_dlopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_name\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    365\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    366\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_handle\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mhandle\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mOSError\u001b[0m: [WinError 126] The specified module could not be found"
     ]
    }
   ],
   "source": [
    "\n",
    "env_filename=os.getenv('ZWO_ASI_LIB')\n",
    "asi.init('C:\\\\Users\\\\peter\\\\Downloads\\\\ASI_Windows_SDK_V3.24.0.0\\\\ASI SDK\\\\lib\\\\x64\\\\ASICamera2.dll')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d20890bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "asi.get_num_cameras()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fb365e16",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0  Cameras\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "Invalid id",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-55925923991b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mnum_cameras\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0masi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_num_cameras\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnum_cameras\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m' Cameras'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mcamera\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0masi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mCamera\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0mcamera\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_control_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0masi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mASI_HIGH_SPEED_MODE\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mcamera\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_roi\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\zwoasi\\__init__.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, id_)\u001b[0m\n\u001b[0;32m    377\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mid_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mint\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    378\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mid_\u001b[0m \u001b[1;33m>=\u001b[0m \u001b[0mget_num_cameras\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mid_\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 379\u001b[1;33m                 \u001b[1;32mraise\u001b[0m \u001b[0mIndexError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Invalid id'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    380\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mid_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msix\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstring_types\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    381\u001b[0m             \u001b[1;31m# Find first matching camera model\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: Invalid id"
     ]
    }
   ],
   "source": [
    "num_cameras = asi.get_num_cameras()\n",
    "print(num_cameras,' Cameras')\n",
    "camera=asi.Camera(0)\n",
    "camera.set_control_value(asi.ASI_HIGH_SPEED_MODE, 1)\n",
    "camera.set_roi()\n",
    "if num_cameras==2:\n",
    "    camera2=asi.Camera(1)\n",
    "    camera2.set_control_value(asi.ASI_HIGH_SPEED_MODE, 1)\n",
    "    camera2.set_roi() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b90177eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\magics\\pylab.py:160: UserWarning: pylab import has clobbered these variables: ['time']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n",
      "  \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n"
     ]
    }
   ],
   "source": [
    "%pylab auto\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "74036d38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7989492266666666\n",
      "remote_microns_per_local_pixel 0.7733333333333333\n"
     ]
    }
   ],
   "source": [
    "pix_size=2.9  #microns/pixel\n",
    "pix_scale = 206.6248*pix_size/750    #in arcsec/pixel\n",
    "print(pix_scale)\n",
    "micron_scale_remote = 206.6248/200  # in arcsec/micron_remote\n",
    "remote_microns_per_local_pixel = pix_scale/micron_scale_remote\n",
    "print('remote_microns_per_local_pixel',remote_microns_per_local_pixel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b5caade6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def smooth(y, box_pts):\n",
    "    box = np.ones(box_pts)/box_pts\n",
    "    y_smooth = np.convolve(y, box, mode='same')\n",
    "    return y_smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "be9e8ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_beacon_and_set_roi_max(camera=camera,n8s=1,exposure=10000,maxpix=None):\n",
    "    camera.set_roi()\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    test=camera.capture()\n",
    "    figure()\n",
    "    imshow(test)\n",
    "    if maxpix==None:\n",
    "        maxpix = np.unravel_index(np.argmax(test, axis=None), test.shape)\n",
    "    start_x = int(maxpix[1]) - n8s*8\n",
    "    start_y = int(maxpix[0]) - n8s*8\n",
    "    print(maxpix)\n",
    "    camera.set_roi(start_x=start_x,start_y=start_y,width=2*n8s*8,height=2*n8s*8)\n",
    "    #test =  camera.capture()\n",
    "    #figure(1)\n",
    "    #imshow(test)\n",
    "    return(start_x,start_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3e9ce5d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import butter, filtfilt\n",
    "import numpy as np\n",
    "\n",
    "def butter_highpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = cutoff / nyq\n",
    "    b, a = butter(order, normal_cutoff, btype='high', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_highpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_highpass(cutoff, fs, order=order)\n",
    "    y = filtfilt(b, a, data)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ccb85ac2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vidcapture(camera=camera,nframes=100,exposure=3000,gain=500,n8s=2,maxpix=[500,500]):\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    camera.set_control_value(asi.ASI_GAIN, gain)\n",
    "    offsetxy = find_beacon_and_set_roi_max(camera=camera,exposure=10000,maxpix=maxpix,n8s=n8s)\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    figure()\n",
    "    imshow(camera.capture())\n",
    "    print(offsetxy)\n",
    "    ims=[]\n",
    "    camera.start_video_capture()\n",
    "    t0=time()\n",
    "    x=[]\n",
    "    y=[]\n",
    "    for i in range(nframes):\n",
    "        test=camera.capture_video_frame()\n",
    "        xypos = star_centroid(test)\n",
    "        x.append(xypos[0])\n",
    "        y.append(xypos[1])\n",
    "        ims.append(test) \n",
    "    camera.stop_video_capture()\n",
    "    dt=time()-t0\n",
    "    outdict={}\n",
    "    outdict['xcen']=np.array(x)\n",
    "    outdict['ycen']=np.array(y)\n",
    "    outdict['images']=ims\n",
    "    print('dt = ',dt)\n",
    "    return outdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "515f15e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vidcapture_star(camera=camera,nframes=100,exposure=3000,gain=500,arcsec_per_pix=0.7989492266666666):\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    camera.set_control_value(asi.ASI_GAIN, gain)\n",
    "    offsetxy = find_beacon_and_set_roi_max(camera=camera,exposure=10000)\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    figure()\n",
    "    imshow(camera.capture())\n",
    "    print(offsetxy)\n",
    "    ims=[]\n",
    "    camera.start_video_capture()\n",
    "    t0=time()\n",
    "    x=[]\n",
    "    y=[]\n",
    "    for i in range(nframes):\n",
    "        test=camera.capture_video_frame()\n",
    "        xypos = star_centroid(test)\n",
    "        x.append(xypos[0])\n",
    "        y.append(xypos[1])\n",
    "        ims.append(test) \n",
    "    camera.stop_video_capture()\n",
    "    dt=time()-t0\n",
    "    outdict={}\n",
    "    outdict['xcen_arcsec']=np.array(x)*arcsec_per_pix\n",
    "    outdict['ycen_arcsec']=np.array(y)*arcsec_per_pix\n",
    "    outdict['images']=ims\n",
    "    print('dt = ',dt)\n",
    "    return outdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b28a8960",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vidcapture_2_star(camera=camera,camera2=camera2,nframes=100,exposure=3000,gain=500,arcsec_per_pix=0.7989492266666666,n8s=2):\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    camera.set_control_value(asi.ASI_GAIN, gain)\n",
    "    camera2.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    camera2.set_control_value(asi.ASI_GAIN, gain)\n",
    "    offsetxy = find_beacon_and_set_roi_max(camera=camera,exposure=10000,n8s=n8s)\n",
    "    camera.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "    offsetxy2 = find_beacon_and_set_roi_max(camera=camera2,exposure=10000,n8s=n8s)\n",
    "    camera2.set_control_value(asi.ASI_EXPOSURE, exposure)\n",
    "\n",
    "    figure()\n",
    "    imshow(camera.capture())\n",
    "    print(offsetxy)\n",
    "    figure()\n",
    "    imshow(camera2.capture())\n",
    "    print(offsetxy2)\n",
    "\n",
    "    ims=[]\n",
    "    ims2=[]\n",
    "    camera.start_video_capture()\n",
    "    camera2.start_video_capture()\n",
    "    t0=time()\n",
    "    x=[]\n",
    "    y=[]\n",
    "    x2=[]\n",
    "    y2=[]\n",
    "\n",
    "    for i in range(nframes):\n",
    "        test=camera.capture_video_frame()\n",
    "        xypos = star_centroid(test)\n",
    "        x.append(xypos[0])\n",
    "        y.append(xypos[1])\n",
    "        ims.append(test) \n",
    "        test2=camera2.capture_video_frame()\n",
    "        xypos2 = star_centroid(test2)\n",
    "        x2.append(xypos2[0])\n",
    "        y2.append(xypos2[1])\n",
    "        ims2.append(test2)\n",
    "    camera.stop_video_capture()\n",
    "    camera2.stop_video_capture()\n",
    "    dt=time()-t0\n",
    "    outdict={}\n",
    "    outdict['xcen_arcsec']=np.array(x)*arcsec_per_pix\n",
    "    outdict['ycen_arcsec']=np.array(y)*arcsec_per_pix\n",
    "    outdict['images']=ims\n",
    "    outdict['xcen2_arcsec']=np.array(x2)*arcsec_per_pix\n",
    "    outdict['ycen2_arcsec']=np.array(y2)*arcsec_per_pix\n",
    "    outdict['images2']=ims2\n",
    "\n",
    "    print('dt = ',dt)\n",
    "    return outdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "1e23eb22",
   "metadata": {},
   "outputs": [],
   "source": [
    "camera.stop_video_capture()\n",
    "camera2.stop_video_capture()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "04237c0b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(646, 1251)\n",
      "(484, 1348)\n",
      "(1219, 614)\n",
      "(1316, 452)\n",
      "dt =  300.4722671508789\n"
     ]
    }
   ],
   "source": [
    "polaris_2x_240124_3ms_13=vidcapture_2_star(camera=camera,camera2=camera2,exposure=3000,nframes=100000,n8s=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "809802da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1cf967bcf48>"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure()\n",
    "imshow(polaris_2x_240124_3ms_13['images'][4])\n",
    "figure()\n",
    "imshow(polaris_2x_240124_3ms_13['images2'][4])\n",
    "figure()\n",
    "imshow(polaris_2x_240124_3ms_13['images'][98000])\n",
    "figure()\n",
    "imshow(polaris_2x_240124_3ms_13['images2'][98400])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "26ac6562",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "save_dict_to_h5(polaris_2x_240124_2ms_12,filename='polaris_2x_240124_2ms_20sec_capture_12.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "d421bfba",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddict=polaris_2x_240124_2ms_9.copy()\n",
    "figure(figsize=[8,6])\n",
    "xlabel('Sample [2ms]')\n",
    "ylabel('Y position [arcsec]')\n",
    "xsdev=np.std(ddict['xcen_arcsec'])\n",
    "ysdev=np.std(ddict['ycen_arcsec'])\n",
    "\n",
    "xsdev2=np.std(ddict['xcen2_arcsec'])\n",
    "ysdev2=np.std(ddict['ycen2_arcsec'])\n",
    "\n",
    "title('Polaris 2 camera positions at 500 Hz, 1/24/24(9), 20sec\\n  xstdev:'+str(xsdev)+' ystdev:'+str(ysdev))\n",
    "plot(ddict['ycen_arcsec'],label='Camera 1')\n",
    "plot(ddict['ycen2_arcsec'],label='Camera 2')\n",
    "legend()\n",
    "savefig('Polaris/polaris_2x_240124_2ms_20sec_9_yy.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "4de1635f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddict=polaris_2x_240124_2ms_9.copy()\n",
    "figure(figsize=[8,6])\n",
    "xlabel('Sample [2ms]')\n",
    "ylabel('X position [arcsec]')\n",
    "xsdev=np.std(ddict['xcen_arcsec'])\n",
    "ysdev=np.std(ddict['ycen_arcsec'])\n",
    "\n",
    "xsdev2=np.std(ddict['xcen2_arcsec'])\n",
    "ysdev2=np.std(ddict['ycen2_arcsec'])\n",
    "\n",
    "title('Polaris 2 camera positions at 500 Hz, 1/24/24(9), 20sec\\n  xstdev:'+str(xsdev)+' ystdev:'+str(ysdev))\n",
    "plot(ddict['xcen_arcsec'],label='Camera 1')\n",
    "plot(ddict['xcen2_arcsec'],label='Camera 2')\n",
    "legend()\n",
    "savefig('Polaris/polaris_2x_240124_2ms_20sec_9_xx.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "id": "725e5e82",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddict=polaris_2x_240124_2ms_9.copy()\n",
    "figure(figsize=[8,6])\n",
    "xlabel('X position [arcsec]')\n",
    "ylabel('Y position [arcsec]')\n",
    "xsdev=np.std(ddict['xcen_arcsec'])\n",
    "ysdev=np.std(ddict['ycen_arcsec'])\n",
    "\n",
    "xsdev2=np.std(ddict['xcen2_arcsec'])\n",
    "ysdev2=np.std(ddict['ycen2_arcsec'])\n",
    "\n",
    "title('Polaris 2 camera positions at 500 Hz, 1/24/24(9), 20sec.\\n  xstdev:'+str(xsdev)+' ystdev:'+str(ysdev)+ '\\n x2stdev:'+str(xsdev2)+' y2stdev:'+str(ysdev2))\n",
    "plot(ddict['xcen_arcsec'],ddict['ycen_arcsec'],'x',label='Camera 1')\n",
    "plot(ddict['xcen2_arcsec'],ddict['ycen2_arcsec'],'x',label='Camera 2')\n",
    "legend()\n",
    "savefig('Polaris/polaris_2x_240124_2ms_9_20sec_xy.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "d61d2081",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1cfe26de5c8>"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nsmth=10 #smoothing bin width\n",
    "hp=0.1  #highpass cutoff in hz\n",
    "fs=500  #sample rate in hz\n",
    "ddict=polaris_2x_240124_2ms_9.copy()\n",
    "figure(figsize=[8,6])\n",
    "xlabel('Sample [3ms]')\n",
    "#ylabel('Y position [arcsec]')\n",
    "#xlabel('Cam 1 X position [arcsec]')\n",
    "ylabel('Cam 1, 2 Y position [arcsec]')\n",
    "\n",
    "x1 = smooth(butter_highpass_filter(ddict['ycen_arcsec'],hp,fs),box_pts=nsmth)\n",
    "x2 = smooth(butter_highpass_filter(ddict['ycen2_arcsec'],hp,fs),box_pts=nsmth)\n",
    "title('Polaris 2 camera positions at 500 Hz, 1/24/24(9), 200sec \\n Highpass %s Hz smoothed %s moving average' %(str(hp), str(nsmth)))\n",
    "#plot(x1,x2,'x',label='Camera 1')\n",
    "plot(x1,label='Camera 1')\n",
    "plot(x2,label='Camera 2')\n",
    "legend()\n",
    "#savefig('Polaris/polaris_2x_11_240110_2ms_100sec_yy.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "3d1e0a97",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'polaris_240109_8_2ms' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-123-dea7d1eb5a56>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mddict\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mpolaris_240109_8_2ms\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mfigure\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mxlabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'X position [arcsec]'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mylabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Y position [arcsec]'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mxsdev\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstd\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mddict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'xcen_arcsec'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'polaris_240109_8_2ms' is not defined"
     ]
    }
   ],
   "source": [
    "ddict=polaris_240109_8_2ms.copy()\n",
    "figure(figsize=[8,6])\n",
    "xlabel('X position [arcsec]')\n",
    "ylabel('Y position [arcsec]')\n",
    "xsdev=np.std(ddict['xcen_arcsec'])\n",
    "ysdev=np.std(ddict['ycen_arcsec'])\n",
    "title('Polaris_8 position at 500 Hz, 1/9/24, 100sec\\nUpgraded centroids:  xstdev:'+str(xsdev)+' ystdev:'+str(ysdev))\n",
    "plot(ddict['xcen_arcsec'],ddict['ycen_arcsec'],'x')\n",
    "#savefig('Polaris/polaris_240109_8_2ms_100sec_xy.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "a5cefa81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1ef22d6a6c8>"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure()\n",
    "imshow(polaris_2x_5_240110_2ms['images2'][100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "b0935011",
   "metadata": {},
   "outputs": [],
   "source": [
    "camera.stop_video_capture()\n",
    "camera2.stop_video_capture()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "c99e5d5d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64.85607709750568\n",
      "67.56395691609977\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x177ab3d2fc8>"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arcsec_per_pix=0.7989492266666666\n",
    "print(polaris_240109_4_5ms['xcen_arcsec'][0]/arcsec_per_pix)\n",
    "print(polaris_240109_4_5ms['ycen_arcsec'][0]/arcsec_per_pix)\n",
    "figure()\n",
    "imshow(polaris_240109_4_5ms['images'][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "id": "1bdeae35",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = polaris_240109_5_2ms['images'] # some array of images\n",
    "frames = [] # for storing the generated images\n",
    "fig = plt.figure()\n",
    "for i in range(len(img)):\n",
    "    frames.append([plt.imshow(img[i], cmap=cm.Greys_r,animated=True)])\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, frames, interval=50, blit=True,\n",
    "                                repeat_delay=1000)\n",
    "ani.save('Polaris_240109_5_500Hz_10sec.mp4')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "1119fc1a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(polaris_2x_240124_2ms_9['images'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "9fe9444c",
   "metadata": {},
   "outputs": [],
   "source": [
    "img1 = polaris_2x_240124_2ms_9['images'][:500] # some array of images\n",
    "img2 = polaris_2x_240124_2ms_9['images2'][:500] # some array of images\n",
    "frames = [] # for storing the generated images\n",
    "fig,(ax1,ax2) = plt.subplots(1,2)\n",
    "plt.title('Polaris 240124 2ms')\n",
    "for i in range(len(img1)):\n",
    "    frames.append([ax1.imshow(img1[i], cmap=cm.Greys_r,animated=True),ax2.imshow(img2[i], cmap=cm.Greys_r,animated=True)])\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, frames, interval=50, blit=True,\n",
    "                                repeat_delay=1000)\n",
    "ani.save('Polaris_2x_240124_2ms_9.mp4')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "bb552a60",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x175eb815688>"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure()\n",
    "imshow(ddict['images'][1300])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "3e13c188",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dt =  11.810209035873413\n"
     ]
    }
   ],
   "source": [
    "x1,y1,ims1,x2,y2,ims=vidcapture2(camera1=camera,camera2=camera2,nframes=10000,offsetxy1=[512,512],n8s=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "19bf35f3",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'camera2' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-465e280a3305>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mdef\u001b[0m \u001b[0mvidcapture2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcamera1\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcamera\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcamera2\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcamera2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnframes\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mexposure1\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m32\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgain1\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m400\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mexposure2\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m32\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgain2\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0moffsetxy1\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mn8s\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m     \u001b[0mcamera1\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_control_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0masi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mASI_EXPOSURE\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mexposure1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m     \u001b[0mcamera1\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_control_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0masi\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mASI_GAIN\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgain1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0moffsetxy1\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m         \u001b[0moffsetxy1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfind_beacon_and_set_roi_max\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcamera\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcamera1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mn8s\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mn8s\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'camera2' is not defined"
     ]
    }
   ],
   "source": [
    "def vidcapture2(camera1=camera,camera2=camera2,nframes=10,exposure1=32,gain1=400,exposure2=32,gain2=0,offsetxy1=None,n8s=8):\n",
    "    camera1.set_control_value(asi.ASI_EXPOSURE, exposure1)\n",
    "    camera1.set_control_value(asi.ASI_GAIN, gain1)\n",
    "    if offsetxy1 == None:\n",
    "        offsetxy1 = find_beacon_and_set_roi_max(camera=camera1,n8s=n8s)\n",
    "    else:\n",
    "        start_x = offsetxy1[0] - 4*n8s*8\n",
    "        start_y = offsetxy1[1] - n8s*8\n",
    "        camera1.set_roi(start_x=start_x,start_y=start_y,width=2*n8s*8,height=2*n8s*8)\n",
    "        camera2.set_roi(start_x=start_x,start_y=start_y,width=2*n8s*8,height=2*n8s*8)\n",
    "    camera2.set_control_value(asi.ASI_EXPOSURE, exposure2)\n",
    "    camera2.set_control_value(asi.ASI_GAIN, gain2)\n",
    "    offsetxy2 = offsetxy1  #find_beacon_and_set_roi_max(camera=camera2)\n",
    "\n",
    "    ims1=[]\n",
    "    ims2=[]\n",
    "    camera1.start_video_capture()\n",
    "    camera2.start_video_capture()\n",
    "    t0=time()\n",
    "    x1=[]\n",
    "    y1=[]\n",
    "    x2=[]\n",
    "    y2=[]\n",
    "    for i in range(nframes):\n",
    "        test1=camera1.capture_video_frame()\n",
    "        beamparam1 = lbs.basic_beam_size(test1)\n",
    "        x1.append(beamparam1[0])\n",
    "        y1.append(beamparam1[1])\n",
    "        ims1.append(test1) \n",
    "        test2=camera2.capture_video_frame()\n",
    "        beamparam2 = lbs.basic_beam_size(test2)\n",
    "        x2.append(beamparam2[0])\n",
    "        y2.append(beamparam2[1])\n",
    "        ims2.append(test2) \n",
    "    camera1.stop_video_capture()\n",
    "    camera2.stop_video_capture()\n",
    "    dt=time()-t0\n",
    "    print('dt = ',dt)\n",
    "    return x1,y1,ims1,x2,y2,ims2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "0742de4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "camera2.stop_video_capture()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "a9325aba",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "freqs1=[]\n",
    "xamplitudes1=[]\n",
    "yamplitudes1=[]\n",
    "for freq in linspace(1,500,100):\n",
    "    freq=int(freq)\n",
    "    print('Frequency',freq)\n",
    "    set_freq_volts(sdg1025,freq,2.0)\n",
    "    x,y,ims = vidcapture(camera,nframes=1000,exposure=32,gain=0)\n",
    "    xamp = max(x)-min(x)\n",
    "    yamp = max(y)-min(y)\n",
    "    xamplitudes1.append(xamp)\n",
    "    yamplitudes1.append(yamp)\n",
    "    freqs1.append(freq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "b856f936",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13d82167b88>]"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure()\n",
    "title('Beacon fiber motion from image at ~1kHz framerate \\n 750mm/200mm Fl telescope, 4V PP on FTA')\n",
    "xlabel('Frequency [Hz]')\n",
    "ylabel('X p-p amplitude [$\\mu$]')\n",
    "legend()\n",
    "plot(freqs1,np.array(xamplitudes1)*remote_microns_per_local_pixel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ebadb2c1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\IPython\\core\\magics\\pylab.py:160: UserWarning: pylab import has clobbered these variables: ['time']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n",
      "  \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n"
     ]
    }
   ],
   "source": [
    "%pylab auto\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "f81b2264",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1cfe19b12c8>"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure(2)\n",
    "camera.set_roi()\n",
    "camera.set_control_value(asi.ASI_EXPOSURE, 10000)\n",
    "camera.set_control_value(asi.ASI_GAIN, 550)\n",
    "test=camera.capture()\n",
    "imshow(test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "d6f9862d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1cfa75b6d88>"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure(3)\n",
    "camera2.set_roi()\n",
    "camera2.set_control_value(asi.ASI_EXPOSURE, 10000)\n",
    "camera2.set_control_value(asi.ASI_GAIN, 550)\n",
    "test=camera2.capture()\n",
    "imshow(test)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "415c5327",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_lagged_cc(aa,bb,lags=arange(-10,11,1)):\n",
    "    #calculate 1d cross correlation over various lags, input by the user as a list, or assumed +- 10 samples, return a dict\n",
    "    #just ignore the ragged edges from shifting timelines\n",
    "    out_cc={}\n",
    "    cc = []\n",
    "    a=aa-mean(aa)\n",
    "    b=bb-mean(bb)\n",
    "    for lag in lags:\n",
    "        if lag == 0:\n",
    "            cc.append(np.correlate(a[lag:],b[lag:]))\n",
    "        if lag < 0:\n",
    "            cc.append(np.correlate(a[-1*lag:],b[:lag]))\n",
    "        if lag > 0:\n",
    "            cc.append(np.correlate(a[:-lag],b[lag:]))\n",
    "    out_cc['cc'] = np.array(cc)\n",
    "    out_cc['lags'] = np.array(lags)\n",
    "    return out_cc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "id": "7b8d2784",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_lagged_pearson(aa,bb,lags=arange(-10,11,1)):\n",
    "    #calculate 1d cross correlation over various lags, input by the user as a list, or assumed +- 10 samples, return a dict\n",
    "    #just ignore the ragged edges from shifting timelines\n",
    "    out_cc={}\n",
    "    cc = []\n",
    "    a=aa-mean(aa)\n",
    "    b=bb-mean(bb)\n",
    "    for lag in lags:\n",
    "        if lag == 0:\n",
    "            cc.append(pearsonr(a[lag:],b[lag:])[0])\n",
    "        if lag < 0:\n",
    "            cc.append(pearsonr(a[-1*lag:],b[:lag])[0])\n",
    "        if lag > 0:\n",
    "            cc.append(pearsonr(a[:-lag],b[lag:])[0])\n",
    "    out_cc['cc'] = np.array(cc)\n",
    "    out_cc['lags'] = np.array(lags)\n",
    "    return out_cc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "id": "49f6905d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x1_hpr1=butter_highpass_filter(polaris_2x_240124_2ms_9['xcen_arcsec'],.1,500)\n",
    "x2_hpr1=butter_highpass_filter(polaris_2x_240124_2ms_9['xcen2_arcsec'],.1,500)\n",
    "y1_hpr1=butter_highpass_filter(polaris_2x_240124_2ms_9['ycen_arcsec'],.1,500)\n",
    "y2_hpr1=butter_highpass_filter(polaris_2x_240124_2ms_9['ycen2_arcsec'],.1,500)\n",
    "sccd9_hp_x = get_lagged_pearson(x1_hpr1,x2_hpr1,lags=arange(-1000,1001,1))\n",
    "sccd9_hp_y = get_lagged_pearson(y1_hpr1,y2_hpr1,lags=arange(-1000,1001,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "6de67549",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1d029534b08>"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure(1)\n",
    "title('Lagged pearson cross correlation 2 camera Polaris measurements \\n  2msec sampling, 240124')\n",
    "xlabel('Lag time [msec]')\n",
    "ylabel('Cross correlation [$arcsec^2$]')\n",
    "plot((sccd9_x['lags']),sccd9_x['cc'],label='X position')\n",
    "plot((sccd9_y['lags']),sccd9_y['cc'],label='Y position')\n",
    "legend()\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "id": "54a231df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1d02bcc6588>"
      ]
     },
     "execution_count": 204,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ddict_list=[polaris_2x_240124_2ms_2,polaris_2x_240124_2ms_3,polaris_2x_240124_2ms_4,polaris_2x_240124_2ms_5,polaris_2x_240124_2ms_6,polaris_2x_240124_2ms_7,polaris_2x_240124_2ms_8,polaris_2x_240124_2ms_9,polaris_2x_240124_2ms_10]\n",
    "nlist=['2','3','4','5','6','7','8','9','10']\n",
    "hp=1\n",
    "fs=500\n",
    "figure()\n",
    "title('Lagged Pearson cross correlation 2 camera Polaris measurements \\n  X position: 2msec sampling, Highpass 1 Hz 240124')\n",
    "xlabel('Lag time [msec]')\n",
    "ylabel('Cross correlation [$arcsec^2$]')\n",
    "for ddict,n in zip(ddict_list,nlist): \n",
    "    print(n)\n",
    "    x1_hp=butter_highpass_filter(ddict['xcen_arcsec'],hp,fs)\n",
    "    y1_hp=butter_highpass_filter(ddict['ycen_arcsec'],hp,fs)\n",
    "    x2_hp=butter_highpass_filter(ddict['xcen2_arcsec'],hp,fs)\n",
    "    y2_hp=butter_highpass_filter(ddict['ycen2_arcsec'],hp,fs)\n",
    "    cc_hp_x = get_lagged_pearson(x1_hp,x2_hp,lags=arange(-1000,1001,1))\n",
    "    cc_hp_y = get_lagged_pearson(y1_hp,y2_hp,lags=arange(-1000,1001,1))\n",
    "    plot((cc_hp_x['lags']),cc_hp_x['cc'],label='Scan '+n)\n",
    "legend()\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "c7057db1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1d029534b08>"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure(1)\n",
    "title('Lagged pearson cross correlation 2 camera Polaris measurements \\n  2msec sampling, 240124')\n",
    "xlabel('Lag time [msec]')\n",
    "ylabel('Cross correlation [$arcsec^2$]')\n",
    "plot((sccd9_x['lags']),sccd9_x['cc'],label='X position')\n",
    "plot((sccd9_y['lags']),sccd9_y['cc'],label='Y position')\n",
    "legend()\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "c6be82ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1513885.57372974])"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.correlate(polaris_2x_240124_2ms_9['xcen_arcsec'],polaris_2x_240124_2ms_9['xcen_arcsec'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "2cf1446c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[14.395612429575756,\n",
       " 14.381086079999998,\n",
       " 14.473481568798183,\n",
       " 14.407717720888888,\n",
       " 15.443766119352748,\n",
       " 13.40124268880503,\n",
       " 12.856775055438595,\n",
       " 13.82467501142857,\n",
       " 15.05712004102564]"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b[slc]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "23d3ed0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cc=np.correlate(polaris_2x_240124_2ms_9['xcen_arcsec'],polaris_2x_240124_2ms_9['xcen2_arcsec'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "ecd0afca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1cfe418e808>]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "figure()\n",
    "plot(cc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "c1581525",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([18521.71589814])"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b726e8c",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}