amarvutha/gaussian_beams.ipynb

## gaussian_beams.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian beam propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beam waist radius, $w_0$\n",
    "\n",
    "Rayleigh range, $z_R = \\frac{\\pi w_0^2}{\\lambda}$\n",
    "\n",
    "Beam radius, $w(z) = w_0 \\sqrt{1 + \\left(\\frac{z}{z_R}\\right)^2}$\n",
    "\n",
    "Complex beam parameter, $q = (z-z_0) + j z_R$ for a waist located at $z_0$\n",
    "\n",
    "Radius of curvature, $R(z) = \\mathrm{Re}\\left(\\frac{1}{q}\\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All lengths in mm\n",
    "\n",
    "from __future__ import division, print_function\n",
    "import numpy as np\n",
    "from numpy import sin,cos,tan,arctan,pi,sqrt,abs\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import optimize\n",
    "from functools import partial\n",
    "\n",
    "speed_of_light = 3e11   # mm/s\n",
    "I = np.mat(np.eye(2))\n",
    "lamda = 0.684e-3        # mm\n",
    "\n",
    "class OpticalElement:\n",
    "    def __init__(self,T = np.matrix(np.identity(2)),label=''):\n",
    "        self.T = I\n",
    "        self.label = label\n",
    "\n",
    "class FreeSpace(OpticalElement):\n",
    "    def __init__(self,d,label=''):\n",
    "        self.T = np.matrix( [[1,d],[0,1]] )\n",
    "        self.label = label\n",
    "\n",
    "class ThinLens(OpticalElement):\n",
    "    def __init__(self,f,label=''):\n",
    "        self.T = np.matrix( [[1,0],[-1/f,1]] )\n",
    "\n",
    "class Cascade(OpticalElement):\n",
    "    def __init__(self,cascade_list,label=''):\n",
    "        \"\"\" optical elements in order from left to right, inputs at left \"\"\"\n",
    "        self.assembly = cascade_list\n",
    "        self.label = label\n",
    "        self.T = I\n",
    "        for element_i in cascade_list: self.T = element_i.T * self.T \n",
    "        \n",
    "def propagate(T,q):\n",
    "    return (T[0,0]*q + T[0,1])/(T[1,0]*q + T[1,1])\n",
    "\n",
    "def z_R(w0,lamda=1e-3):\n",
    "    return pi*w0**2/lamda\n",
    "\n",
    "def q(z,w0):\n",
    "    return z + 1j*z_R(w0)\n",
    "\n",
    "def R(q):\n",
    "    return 1/((1/q).real)\n",
    "\n",
    "def w(q,lamda=1e-3):\n",
    "    return sqrt(lamda*q.imag/pi) * abs(q)/q.imag\n",
    " \n",
    "def eigenmode(T):\n",
    "    a,b,c,d = np.array(T.flatten())[0]\n",
    "    return np.roots([c,d-a,-b])\n",
    "\n",
    "def make_mode_matching_system(parameter_array):\n",
    "    f1,f2,d2 = parameter_array\n",
    "    d1 = f1+ f2 + 5\n",
    "    return Cascade( [ThinLens(f1),FreeSpace(d1),\n",
    "                    ThinLens(f2),FreeSpace(d2)] )\n",
    "\n",
    "def mode_matching_function(q_in,q_target,parameter_array):    \n",
    "    system = make_mode_matching_system(parameter_array)\n",
    "    q_out = propagate(system.T,q_in)\n",
    "    return (q_target.real - q_out.real)**2 + 3*(q_target.imag - q_out.imag)**2 \n",
    "    \n",
    "\n",
    "def mode_matching_nagourney(w_laser,q_target,distance_from_cavity_input=50):\n",
    "    d0 = (-q_target.real) + distance_from_cavity_input\n",
    "    print(d0)                              # distance of lens2 from cav waist\n",
    "    w0 = sqrt( lamda*q_target.imag/pi )    # cavity waist\n",
    "    w2 = w0 * sqrt(1 + (d0/q_target.imag)**2)      # beam size at exit of lens2\n",
    "    print(w2)\n",
    "    f_ratio = w_laser/w2\n",
    "    return d0, f_ratio\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian beam propagation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEOCAYAAACuOOGFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4lHeW4PvvUU4oo5zJQSQDBhOMs40D7uTUzgH3hN29d+7snZ71fWZ3em/f8czshJ57Z6fXbXvancbtdoONAds4G4yNCQKRk4SQShlllKVz/3hLQjAySFClqpLO53n0UPW+r946SFCnfun8RFUxxhhjrlWQrwMwxhgzPlhCMcYY4xGWUIwxxniEJRRjjDEeYQnFGGOMR1hCMcYY4xGWUIwxxniEJRRjjDEeYQnFGGOMR1hCMcYY4xEhvg7A25KTkzUvL8/XYRhjTEDZu3dvvapOHs33jPuEkpeXx549e3wdhjHGBBQRKRvt91iXlzHGGI+whGKMMcYjLKEYY4zxCEsoxhhjPMISijHGGI+whGKMMcYjxn1CqWzuoKKx3ddhGGPMuDfuE8q5tm5u/NtP+U+vF3G4stnX4RhjzLg17hc2xkaEoKq8vb+St/dXsmpaMutXF7ByajIi4uvwjDFm3Bj3LZS0uEh+9thi7pufQURoENtP1vPYK19z9z/t4O39Lnr7+n0dojHGjAuiqr6OYZCI3An8BAgGXlbVFy85/wPgj4A+oA1Yr6pHLnfPwgWL9DebPwGgrbOXrYeqeKe4kqb2HgAy4yN5ZmU+Dy7JJjp83DfYjDFmRERkr6ouHtX3+EtCEZFg4ARwG1AB7AYeHpowRCRWVVvcj+8D/lBV77zcfYcmlAHdvf18cryWjUUuXE0dAMRFhvLYslyeuCGPyZPCPflXM8aYgHM1CcWfPpIvBU6pagmAiLwOrAMGE8pAMnGLBq4qG4aFBHHHnDRum53KrtIGNu6r4Gh1K//fJ6d4aXsJ31mUxXOr8imYHHMNfx1jjJlY/CmhZALlQ55XANdfepGI/BHwJ0AYcPO1vGCQCMsLklhekMSRqhY27KtgV2kD//b1WV7ffZbbZ6fy/I1TWJSTcC0vY4wxE4I/JZQRUdV/Bv5ZRB4B/i/giUuvEZH1wHqAjKzsEd13dnoss++eTXljO28Vufj4WC3vH67h/cM1LMlL4PnVU7h5ZgpBQTYzzBhjhuNPYyjLgf+mqne4n/85gKr+1TdcHwQ0qmrc5e473BjKSDSc72ZzcSVbD1VxvqsPgCmTo3l+9RTWLcwgPCR41Pc0xphAEeiD8iE4g/K3AC6cQflHVPXwkGumqepJ9+N7gf96pb/w1SaUAe3dvWw7UsPb+yupb+sCIGVSOE+uyOP71+cSFxl61fc2xhh/FdAJBUBE1gL/iDNt+FVV/bGI/AjYo6qbROQnwK1AD9AI/PHQhDOca00oA3r7+tl+qp4N+yo4c84p5RIdFszDS3N4emU+GfGR1/waxhjjLwI+oXiDpxLKAFWl6GwTG4oqOFDhlHIJCRLum5/B+hsLmJkW67HXMsYYX7GEMgxPJ5ShTtW2sbGogh2n6ul3/xhvnD6Z528sYHlBkpV2McYELEsow/BmQhlQ3dLJ2/tdfHCkhq5ep5RLYWYc61cXcNfcNEKCx32FG2PMOGMJZRhjkVAGtHT0sPVQFZuLq2jucEq7ZCdG8uzKAr63OIuosICbpW2MmaAsoQxjLBPKgK7ePj4+5pR2qWruBCAhKpTHlufxxPJckmKstIsxxr9ZQhmGLxLKgL5+5auSc2woquBETRsA4SFBfG9xFs+tKiA3KdoncRljzJVYQhmGLxPKAFXlcGULG4oq2H2mEYAggTvnpvH86inMz473aXzGGHOpQC8OOW6JCHMz45ibGUfZufNsLHLx2Yk6th6sZuvBaq7PT+QHN05hzYzJNjPMGBOwrIXiI+fauninuJJ3D1XT3u2UdpmeGsP61VO4b34GYSE2M8wY4zvW5TUMf00oA8539fL+4WrePlBJw/luANJiI3h6ZR4PL81hUoSVdjHGjD1LKMPw94QyoKevn89P1LGhyMXZBqe0y6TwEB5ZlsPTK/JJjY3wcYTGmInEEsowAiWhDFBV9pY18vt9FRyqdPYTCw0W1i3IZP3qAqanTvJxhMaYicASyjACLaEMdaKmlQ1FLr48faG0y80zU3h+dQFL8xNtAN8Y4zWWUIYRyAllQGVTB2/td/HR0Vq6+5zSLvOz4/nB6gJun5NGsG36ZYzxMEsowxgPCWVAc0cPm4sr2XKwitbOXgDykqJ4dlUB370ui4hQ2/TLGOMZllCGMZ4SyoDOnj4+PFrDW/td1LQ4m34lRYfx2PJcHltmpV2MMdfOEsowxmNCGdDXr+w8Xc+GfS5O1V0o7fLd67J4ZmU+BZNjfByhMSZQWUIZxnhOKANUlUOuZjYUudhT5pR2EYHbZqWyfnUB1+Um2AC+MWZUrPTKBCUiFGbFU5gVz9mGdt7a7+KTY7VsO1LDtiM1LMyJZ/0qG8A3xniXtVDGqcbz3Ww5WMXWg1W0djkD+DmJUTyzMt/2ZjHGXJF1eQ1joiaUAZ09fXx0tIa39ldS3eLszRIXGcpjy3J5/IZcUibZCnxjzL9nCWUYEz2hDBjYm2VjkYvjNa0AhAUH8a2FmTy7Kp9ptgLfGDOEJZRhWEL5945WtbCxyMVXJecY+O3fPDOF51YVsKzAVuAbY8bBoLyI3An8BAgGXlbVFy85/yfAs0AvUAc8raplYx5ogJuVHsus9NiLVuB/fMz5mpsZy3OrClhbmE5osJXQN8aMnN+0UEQkGDgB3AZUALuBh1X1yJBrbgJ2qWq7iPwBsEZVH7zcfa2FcmXNHT1sPVjFloNVNHf0AJAZH8lTK/J4aGkOMeF+9bnDGDMGrqaF4k8fQZcCp1S1RFW7gdeBdUMvUNVPVLXd/fQrIGuMYxyX4iJDeXhpDq88sZg/WjOVzPhIXE0d/N9bjrL8rz7ir949SnVzp6/DNMb4OX/66JkJlA95XgFcf5nrnwHe9WpEE0x4SDB3zk3j9jmp7DnTwIYiF4crW/hfn5XwyvZS7pufwbOrCpidEevrUI0xXlTirrwxWv6UUEZMRB4FFgM3fsP59cB6gIys7DGMbHwIEmFpfhJL85M4UdPKxiKXU+KlyMWGIherpiXz3KoCVk1LtgF8Y8YJVeXr0gZ+tr2Uj47VXNU9/CmhuICh7/5Z7mMXEZFbgReAG1W1a7gbqepLwEvgjKF4PtSJY3rqJP7szplUt3Syab+LD47WsP1kPdtP1jMzbRLPrSrg3vkZhIX4U++pMWakevr62Xqwild2lFJc0Qw4m/pdDX8alA/BGZS/BSeR7AYeUdXDQ65ZCLwJ3KmqJ0dyXxuU96y2zl7ePVzF5gNVNLR3A5AaG85TK/J5eGkOcZGhPo7QGDMSLZ09/Pbrcn6+8wyupg4AYiNCWFuYztrCdG6ckRLY61BEZC3wjzjThl9V1R+LyI+APaq6SUQ+BAqBKve3nFXV+y53T0so3tHT189nJ+rYWOTibIMzTyI6LJgHl+Tw9Mo8shKifByhMWY4rqYO/nVHKa/vLqfNXZYpMz6SdQsyuHlmCuEhzr5K87LjAzuheIMlFO9SVfadbWJjUQUH3M3l4CBhbWE6z63KZ15WvI8jNMYAHChv4mfbS3j3UDV97j3FCzPjuH9BJovzEgi6ZDz0ahKKP42hmAAkIlyXm8B1uQmcrmvjrSIXn5+s450DlbxzoJKl+Yk8uzKfW2alWqVjY8ZYX7/y0dEaXt5eytdnGgDnA9+a6ZNZtyCTqSme3TPJWijG4+pau3inuJL3DlXT0dMHQG5SFE+vyOe712URbQsljfGqju4+3txbzis7Sjlz7kKX9B1z0rhnXgaTJ115V1fr8hqGJRTfae/u5YMjNWw6UEltqzMhLzYihIevz+GJ5XlkxEf6OEJjxpfa1k5+sbOMX+0qo6ndqXqRMimc++ZncNvs1FFtW2EJZRiWUHxvoNLx2wcqOVrVAjjN7rsL03lmZT7zs22cxZhrcay6hZe3l7JpfyXdff0ATE+N4f4FmdwwJfmqupttDMX4peAgYcXUZFZMTeZETStv73ex41Q9mw5UsulAJUvyEnhmZQG3zbZxFmNGSlXZfrKen20vYfvJegAEWF6QxP0LM5mVNmnMFx5bC8X4RG1rJ5uLq9h2uJrz3c44S05iFE+tyON7i7OtIKUx36Czp4+3ily8+kUpJ2qcEinhIUHcNiuVe+dneKwr2bq8hmEJxb+1d/fy4dFa3jlwYUfJSREhPLw0hyduyCPTxlmMAaC2pZNfflXGr3edpeG8s6g4MSqMe+alc+fcNCZFeHZRsSWUYVhCCQx9/crXpc44y+HKC+Msd81N49lVBSywcRYzQR1yNfPqjlLeKa6kp895v546OYZ1CzJYMTXZa/sW2RiKCVjBQcLyKcksnzIwzlLJjlN1bC6uYnNxFdflJvDMynxun51KiG38Zca5vn7lgyM1vPpFKV+XOutHgsQZH1m3IIPZ6bF+WZjVWijGb9W1drHlYCXvHa7mfJczzpKVEMmTN+Tx4JJsjzfxjfG11s4e3thTwc93llLe4NTXigwN5vbZqdwzP4O02Igxi8W6vIZhCSXwdXT38dExZz1LlXujr5jwEB5aks0TN+SRnWh1w0xgK29o51+/OMMbey7U10qLjeDe+encOmt060c8xRLKMCyhjB99/cruMw28td81OM4SJHDn3DSevCGfJXkJftkNYMxwVJXdZxp5ZUcJHxypwV1ei7kZsdy3IJOleYk+nUZvYyhmXAsOEpYVJLGsIIlTtW28vd/F9lP1bD1YzdaD1czJiOXJG/K4d34GEaHBvg7XmGF19/az5WAlr+wo5ZDL+WAU4q6vde/8DI/X1xpL1kIxAe1cWxdbD1Xz3qEqWjqdroKk6DAeuT6HR5flkjqGfc7GXE7D+W5+s6uMX3xZdlEporsK01k7N53E6DAfR3gx6/IahiWUiaG7t3+wynFJ/XnA+dR3V2E6T63IY2F2vHWHGZ845GrmtZ1nePtAJd29TlmU3MQo7luQwY3TJw/uP+JvvNLlJSKJI7hPv6o2jeaFjfGksJAgbp2Vyi0zUzhS1cI7Byr5suTcYBn9+VlxPLkij7sLbbti4309ff28f7ian39xhj1ljYPHF+cmsG5BJvOz4sblB5wrtlBEpBOoxCkT802CVTXHk4F5irVQJq7a1k62Hqxm25FqWt3dYZMnhfP963N45PocUiZZd5jxrPq2Lv5t11l+vevsYOWHqLBgbp2Vyt2F6QFVYdsrXV4iUqSqC6/1Gl+xhGI6e/r47ITTHVbm3q44NFi4d14GT67Is10lzTUrrmji5zvPsPlA1WC136yESO6Zl8HNM1KIDPPPbq3L8VZCiVDVzmu9xlcsoZgBqspBVzPvFFeyq6SBgX/5i3LieXJFPnfNTfNaGQsz/nT39vPuoSpe23mGfWedHn8BluQlcs+8dBYE+LidV8ZQRpIo/DWZGDOUiDAvK555WfFUt3SypbiKD45Us+9sE/vOFpEyKZyHlzrdYTY7zHyT2tZOfuPu1qpzz9aKDgvmttmp3F2YQVrcxP23M+JZXiKyGHgByMVJRAKoqs7zXnjXzloo5nI6uvv49IRT7bi80Sl1ERwk3DEnlUeX5bK8ICmgP2Uaz1BV9p1t5JdflrHlYNVgkcbsxCjunZfOmumB2a11OV6dNiwix4H/DBwE+geOq2rZaF5wrFlCMSMx0B229WAVX5acG1y1PGVyNI8ty+Xb12URa7XDJpy2rl7eKnLxq6/KOFbdCjjVGZbmJ3LPvAzmZY7P2Vrg/YSyQ1VXXlVkIw1G5E7gJ0Aw8LKqvnjJ+dXAPwLzgIdU9c0r3dMSihmtc21dbDtSw3uHqwf3nYgMDeb+hZk8tiyX2RmxPo7QeNux6hZ+9VUZbxVVDtbWio0I4bbZadw5N21MizT6ircTyi3Aw8BHQNfAcVXdMJoXvMz9g4ETwG1ABbAbeFhVjwy5Jg+IBf4U2GQJxXhTb18/u0ob2HqwimJX8+Dx63ITeGxZLncVpvntojQzel29fbx3qJpffll20dqR2emx3DU3zat7j/gjb9fyegqYCYRyoctLAY8kFGApcEpVSwBE5HVgHTCYUFT1jPtc/3A3MMaTQoKDWDE1mRVTkylvaGfroSo+PlbL3rJG9pY18t83h/HAkmweWZpjFY8DWHlDO7/edZbf7Snn3JAW6U0zU7hrThp5ydE+jjBwjCahLFHVGV6LBDKB8iHPK4Drvfh6xoxYdmIUz6+ewuPL8vjsRB1bDlZy5lw7//Lpaf7l09OsmpbMg0uyuW12qrVaAkBPXz8fHa3l9d1n+exEHQMdNXlJUawtTOfG6ZN9UjI+0I3mJ7ZTRGYP7YLyVyKyHlgPkJGV7eNozHgSGRbMnXPTuGNOKseqW9l6qIovTtWz/aTzlRgdxrcXZvLQ0mympkzydbjmEqX153l991l+v9dFfZvTcx8SJKycnszauenMTJs0bgfZx8JoEsoy4ICIlOCMoXh62rALGPrun+U+Nmqq+hLwEjhjKNcemjEXExFmpccyKz2W51dN4ZPjtWw7Us2Zc+28vKOUl3eUsjg3gYeW5nB3Yfq4m1IaSDp7+nj3UBWvf13OLvd2ugDZCZHcPjuNm2amEBdpM/g8YTQJ5c5hjnnyzXo3ME1E8nESyUPAIx68vzFeERMRwr3zM7hnXjona9vYdriaz0/Ws6eskT1ljfzlpsOsW5jBQ0tymJsZ5+twJ4xDrmbe2FPOxiLXYC238JAgVk1L5vbZadYa8YJrWdgIgCcXNorIWpxpwcHAq6r6YxH5EbBHVTeJyBJgI5AAdALVqjrncve0WV7GFzq6+9h+qo5th2s4XtM6eHxm2iS+vSiTdQsybTW+F9S1drHpQCUb9lUM7uoJMC0lhttnp7F6erKNjYyQLWwchiUU42tn6s+z7Ug1nx6vo9W9piFIYOW0yXx7YSa3z/HNnuHjRWdPH9uO1LBxXwWfn6ynz70qNSY8hDUzJnP77FTykwN3F0RfCfiFjd5gCcX4i56+fvaUNfLJsVp2n2mg1/3GFx0WzF2F6Xx7YSbLCpII8uE+4oGiv1/ZVdrAxqIKth6sHlx8GBwkXJeTwM0zU1iSl2h731yDgF7Y6C2WUIw/aunoYcepej4+VntRl1h6XAR3zU3n7nnpLMyOt+QyhKqyv7yJLcVVbD1YRWXzhZq001JiuGlGCqunT7YBdg/xdkL5Fc7CxsMMWdioqk+PKsoxZgnF+LvKpg4+Pl7LJ8dqB/caBye5rC1MZ23hxE0uqsqBCqfG2pbiKlxNHYPnkmPCuWnGZG6amUJ2gi0s9TSvj6F4eWGjV1hCMYGiX5UTNa3sOFnPF6frqW/rHjyXERfBnXPTuXW205UznkuA9PT1s/tMAx8dreX9w9VUNF5IIonRYaycmszKqcnMSJtEkM3S8hpvJ5R/Bf42EBY2DmUJxQSiflVOVLey/VQ9Oy9JLpPCQ1g9YzK3zExhzYwUEqPDfBipZzS1d/Pp8To+OlbLp8drB6f5AiRGhbFiahIrp01mpiWRMePthHIUmAKU4p2FjV5hCcUEun5Vjle38lXJOXafaRjctwWc2WILsuNZMTWZ5VOSWJSTQESo/y+i7Ozpo+hsEztP17Pz9Dn2lzcNzs4CZ9HhkrxEluYnMis91pKID3g7oeQOd9ymDRsztqqaO9h9ppHdZxo45GoenC0GEBYSxKKceJYXJLM4L4HCrDi/2MelpbOH4vJm9pc38mXJOfacaaSr90KN1+AgYU56LEvyE1mal0hGfKQPozXg5YQSqCyhmPGsvbuXQ65miiuaKXY1U1p//t9dUzA5mgVZ8RRmxTEtZRJTU2JIjQ33yipxVaWutYsTNW2crG3lcGUL+8ubOF3XxqVvNXlJUczLimd+VhxzMuKIDre1OP7EK+XrRWSfqi661muMMZ4XFRbC0vwkluYnAc505EOVzRx0NXOypo3TdW2U1J2npO48G4oulMabFB5CQUoMWQmRpMVGkB4XQWpsBPFRocSEhxATHkJ0eAghQYICqqAo57v6aO3sobWzl5bOHmpauqhq6qCquRNXUwcldW20DBn/GBASJBSkRDM9dRKz02MpzIwjPirwx37MxUbykWCWiBRf5rwAVqDIGD8QGxnKDVOSuWFKMuDMmCo7186JmlZK6toob+ygvLGd1s5eDpQ3caC8yeMxxISHkJMYRXZiFLmJUcxIm0R+cvS4nplmHCNJKDNHcE3ftQZijPG80OAgpqbEMDXl4tIjzR09VDS2U9faxbnz3Zxrc/4839VLe3cfHT19tHf3MdAlLggIRIQEEe1uvUSFBZMQFUZyTDjJMWFMnhRORlwk8VGhVnRxgrpiQvH3QXdjzOjFRYYSF2kdC8azrA1qjDHGIyyhGGOM8YhRJxQRiRYR/185ZYwxZkxdMaGISJCIPCIiW0SkFjgGVInIERH5WxGZ6v0wjTHG+LuRtFA+wSm58udAmqpmq2oKsBL4CvhrEXnUizEaY4wJACOZNnyrqvaIyHdxdmsEQFUbgN8DvxcR39d2MMYY41NXbKGoao/74S+B3wwdPxGRpy65xhhjzAQ1mkH5Y8BnXNwi+Q+eD8kYY0wgGk1CUVX9KbAB2CQikThlV4wxxphRJZRGAFX9BfAKsAXw6L6bInKniBwXkVMi8sNhzoeLyG/d53eJSJ4nX98YY8zVG3FCUdVbhjx+E/h7IMlTgbjHZv4ZuAuYDTwsIrMvuewZoFFVpwL/APy1p17fGGPMtRnJOpRhu7VUdbOqJl/umlFaCpxS1RJV7QZeB9Zdcs064DX34zeBWzz02sYYY67RiNahiMh/EJGcoQdFJExEbhaR14AnPBBLJlA+5HmF+9iw16hqL9CMB1tJxhhjrt5I1qHcCTwN/JuI5ANNQAQQDGwD/lFVi7wX4uiJyHpgPUBGVraPozHGmIlhJOXrO4H/CfxP93ThZKBDVT29M48LGPrun+U+Ntw1FSISgrOx17lhYn4JeAmcLYA9HKcxxphhjHhQXkQiVLVHVau8kEwAdgPTRCRfRMKAh4BNl1yziQvda98FPla9dKdqY4wxvjCaacNfi8jfeasYpHtM5I+B94GjwBuqelhEfiQi97kvewVIEpFTwJ8A/25qsTHGGN8YyRjKgAXA3cA/iEgQ8C/AFk+2EFR1K7D1kmN/MeRxJ/A9T72eMcYYzxlNCyUeOAz8Jc5q+b8BSrwRlDHGmMAzmhZKPfAl8AXQijPo3eKNoIwxxgSe0bRQFgMngELgCPBPqvqqV6IyxhgTcEZTemWfqj4FPApMBT4Xkf/itciMMcYElBF3eYnIZ0A0FwpC9uNM3f1/vBCXMcaYADOaMZTHcVbJt9jaD2OMMZcacUJR1TJvBmKMMSawjWZQ3hhjjPlGllCMMcZ4xGjGUC4iIulAg6p2eTAeY4wPdPb00XC+m/NdvbT39NHR3UdHTx/9/c5wqbPrkBARGkR0eAjRYSFEhQWTEBVGZFiwT2M3/uOqEwrwS2CKiPxeVf/UUwEZY7yjtbOH8sYOyhvaqWhsp6Kxg/q2Lurbumnr6r3q+0aHBzM5JpzkmHAy4iPJTogiJymKnIQoYiKu5S3GBJqr/m2r6q3u3RIv3abXGONj7d29nK5t40RtGydrWjlR20Zd6zd3JoQFB5ESG058VCgx4SHEhIcQHR5CcJCAggL9qrR399Ha2UNrZy8tnT3UtHRxvquP813tnDnXDmWNF903KTqM6amTmJYaw4zUSUxNiSEqzJLMeDWadSjhwHeAvKHfp6o/8nxYxpjR6Oju40hVCwddTRRXNHO6ro3+Syb3R4UFM2VyDFNTnK8pk6PJSogiLS6CxKgwgoJGv5u2qtJwvpuq5k5cTR2U1J3nZG0rJ2vaOFXbxrnz3XxZco4vS5xtiwTInxzNvMx45mfHMSc9zrrMxpHRfFR4G2fL3b2AjZsY40Oqyum68+w+00DR2UZO1LbRNySDhAQJczNjmZcVx/yseBZkx1MwOcZpcXiQiJAUE05STDhzM+MuOtffr5SeO8/+s03sL3e+jla1UFJ3npK687y130VwkDA9JYb52fEsyUtkakoMQeLZGM3YkZGuURSRQ6o618vxeFzhgkX6m82f+DoMY65ZZ08fByqa2F3awO6yRhrOdw+eCxIozIpneUESy6cksTg3gehw/+ta6ujuY29ZIztP17Pz9DmKK5ouakklRIWyOC+RpXmJLMiOJyLUWi++Mi87fq+qLh7N94zmX9xOESlU1YOjjMsYc5Xau3vZVdrAjpP1FJU30tN34d03NTacm2emcvPMFK4vSCQ2ItSHkY5MZFgwK6cls3JaMgAtnT3sKmngsxO1fHy0lsrmTj44UsMHR2oICw5iXlYcK6cmc31BEjF+mCDNxa7YQhGRgzhjciHANJw9ULpwukNVVed5O8hrYS0UE2jau3v5urSBHafq2Xf24iQyPzueW2amcPPMFOZkxCLjqHtIVTla1cpHR2v48FgtB8ov7DQeEiQszIln5dTJXJ+f6Jetr/HmalooI0kouZc77+8lWSyhmEDQ1dvHrpIGPj9Zd1ESEYEleYncXZjOnXPTSI2N8HGkY6e2tZNth2vYerCKr0rODXaNhQQJi3ISWDUtmWUFSdYt5iVeSSiDF4r8tar+2ZWO+RtLKMZf9atyuLKFT47X8sWpetq7+4CLk8hdc9NImUBJ5JvUtXbx3uFqthRXsqu0gYG3rcjQYJZPSeLmmSkUZsbZgL4HeTuh7FPVRZccK7YuL2NGp6KxnU+O1/Hp8Vpqh6wNmZ8Vx/0LM7m7MN2SyGXUtnby7sFqNha52D+kWyw5Jow101O4aWYKOYlRl7mDGQlvdXn9AfCHQAFwesipScAXqvroaAMdS5ZQjD9o7ezh85P1fHKsluM1rYPHM+Ii+NaiTL61MIupKTE+jDAwldS1sbHIxcYiFxWNHYPHp06O4aaZKdw0YzKTAmCygj/yVkKJAxKAvwJ+OORUq6o2jDrKMWZYNLEbAAAZeklEQVQJxfiKqnLI1cy2IzV8cbp+cFwkOiyYtYXpfHtRFtfnJ17VgkJzsf5+ZU9ZIxv2VbDlYBWtnU4pmdBg4YYpydw+O5W51iU2Kl7t8vImEUkEfouzCv8M8ICqNg5z3XvAMmCHqt4zkntbQjFjrbG9m4+O1rLtSDVVzZ2Dx1dNS+a712Vx++w0Wx3uRZ09fXx4tIY39lSw/WTd4HhLelwEt85K5ZaZKSTFhPs2yADgrRbKn1zuvKr+/Whe8Bte429wKhe/KCI/BBKGG+wXkVtwtiB+3hKK8Sd9/UpReSPbDtfw9ZmGwVXrqbHhPLA4mwcWZ5Nt/fpjrqKxnd/tqeB3e8qpdCf3IIHFuYncPieVxbmJHq8eMF54a2HjJPefM4AlwCb383uBr0fzYpexDljjfvwa8Cnw7xKKqn4kImsuPW6Mr9S2dvLhkRo+OFpLfZszwB4cJNw2O5WHlmRz4/TJhATbtkO+kpUQxf9+23T+4y3T2H6yjt/uLueDI07S//pMA4lRYdw6O5U75qSSMskmQlyrKyYUVf1LABH5HFikqq3u5/8N2OKhOFJVtcr9uBpI9dB9jfG4flWKzjax9WAVe8oaBtdH5CRG8eCSbL57XdaEWi8SCIKDhDUzUlgzI4X6ti427Kvg9d3llNSd54095by5t5zFuYmsLUxnYU68jbVcpdEsN00Fuoc872YUb/wi8iGQNsypF4Y+UVUVkWsa2BGR9cB6gIys7Gu5lTGDWjp6+PBoDe8dvjA2EhosrJ2TxsNLc1hekGQD7AEgOSac9aun8NyqAnafaeRXX5Xx7qGqwVZLWmwEd85N49ZZqcRF2gyx0RhNQvkF8LWIbMQpu3I/TvfUiKjqrd90TkRqRCRdVavcO0HWjiKu4V7rJeAlcMZQruVexpyoaWXLwSq2n6wbnKmVGR/JI9fn8OCSbJJtgDcgiQhL8xNZmp9Ifdts3thTzq+/OourqYOf7zzDr74qY+XUZO4qTGdW2qRxVebGW0acUFT1xyLyLrDSfegJVd3voTg2AU8AL7r/fNtD9zXmqnT29LHjZD1bDlVxqrZt8PiN0yfz2LJcbpqZYoO540hyTDh/uGYqz6+ewmcnavnVV2f55Hgtn56o49MTdeQlRXHX3HTWzJhsG4Rdxkhmee1Q1ZUi0opTJHLo/yJV1dhrDkIkCXgDyAHKcKYNN4jIYuAHqvqs+7rtwEwgBjgHPKOq71/u3jbLy4xGZVMH7x6q4sOjtYPb4sZHhfLA4mweWZpDXnK0jyM0Y6W8oZ3ffH2WN3aXc869VUBkaDA3zUzhnsL0cT9rL2DXoXiTJRRzJX39yp6yBrYUV1E0pJTH/Ox4HluWyz3z0q0A4QTW1dvHe4eq+fVXZ/n6zIW13Auy47l3XjrXjdOpx17dD0VEfgl8DmxX1WOjDc4Yf9PW1cuHR2rYfLCSmhZnym94SBDrFmTw6LJc5mXF+zhC4w/CQ4JZtyCTdQsyOVbdwms7y9hYVDG4C2VabAR3F6Zz6+zUCb9ny2iKQ94ErHJ/TQGKgM9V9SfeC+/aWQvFXKq8sZ3NxVV8fKyGzp5+ALITI3l8WR7fW5xFfFSYjyM0/q6pvZs39pTziy/LBmuIhYcEcfPMFO6ZlzEuilN6vctLRIJxFjfeBPwA6FDVmaOKcoxZQjHgrB3ZW9bIOwcqL+rWWjE1iSdvyOdmG2Q3V6GvX/n4WC0/31nKF6fODR6flxXHvfMyWJIXuN1h3u7y+giIBr4EtgNLVPWapvca423t3b18eLSWzcWVg2tHIkKD+NbCLJ68IY8ZaZOucAdjvtlAVYTbZqdyoqaV13aeYcM+F8UVzRRXNJMyKZy7C9O5bXbqhKh6PJour38ArsPZ/vcLnPGUL1W147Lf6GPWQpmYKps62FxcyYdHa+nocTauyoyP5PHluTy4JNu6tYzXNHf08Dt3d9jZhnbA6Q5bM30y98zLCJiZgmMyy0tEJgFPAn8KpKmqX6/qsoQycfSrsv9sE+8UV7Kn7EKx6uvzE3lqRR63zkq1ulpmzPT1K58er+XnO8+w/WT94PH5WXHcNz+TxXkJfl3ixdtdXn+MMyB/HU6J+Vdxur6M8amO7j4+Pu50aw0MkIaFBHH/ggyeuCGPORlxPo7QTETBQcIts1K5ZVYqp2rb+MWXZ3hzbwUHKpo5UNFMRlwE983P4OaZqeNmO4PRdHn9KU4C2auqvV6NyoOshTJ+1bR0srm4kg+O1HDevR97WmwEjy3P5eGlOSRGW7eW8S/NHT38dvdZXttZhqvJ+fATHR7MHbPTuGdeBpMn+U+Hjy1sHIYllPHnWFULb+138WXJucFKv9flJvDUijzumJNGqHVrGT/X29fP+4drePWLUva6u2eDBFZMTea++RnMTLvmAiTXzKtdXsb4Ul+/svN0PW/vrxzckz0kSLhvfjpPr8y3RYgmoIQEB3H3vHTunpfO/vImXt1RytaDVWw/Wc/2k/XMSJ3EugUZ3DAlOaCmHVsLxfi18129bDtSzTvFVdS1OqvZ4yJDeeT6HJ5YnkdanO07YsaHquYOfvFlGb/ZdZbmjh7AKVp5z7x07pidRkzE2H7+ty6vYVhCCUzVzZ284x4fGZj2m58czdMr8vjOdVlW8dWMW+3dvWzY5+LVL0opqTsPONOOb5mVyn3zMshMiByTOCyhDMMSSuBQVY5Wt/JWkYtdpRfGR5YXJPHMSmc1u21gZSaK/n7ls5N1vLqj9KJpx4tzE7h/QSbzsuK8ukeLjaGYgNTb18/O0+d4a7+Lk+69R0KDhfvnZfD0ynzmZtq0XzPxBAUJN81I4aYZKRyvbuVfvyhlQ5GLPWWN7ClrJDcxivsWZLBmegphIf4xEcVaKMZn2rp62Xa4mneKK6lvc/abiI8K5fvX5/D48jzbl92YS5xr6+I3u87yi6/KBscU4yNDWVuYztrCdI9uWWxdXsOwhOJ/qpo72HSgkg+PXqj2W5AczdMr8/nOoqxxs8jLGG/p7u1nc3Elr+wo5XBlCwBhwUHcNDOF+xdkkJVw7dWOLaEMwxKKf1BVjrjXj+wqaWDgX92Kqc74yJrpNj5izGipKl+WnOPl7aV8fOxCrd4leQl8a0EmczOvfpzFxlCM3+nt62fHKWf9yKm6C+Mj6xZk8vSKfGZn+H4BlzGBSkS4YUoyN0xJ5lRtG6/sKGXDvgp2n2lk95lGpkyO5v4FmaycmjwmdeyshWK8oq2zl/cOV7O5uHJwP+6EqFAeXZbLY8tySbHxEWO84lxbF7/8qoxfflk2+H8vOSaMe+dlcPuctBHvKmldXsOwhDK2alo62XSgkm1HqgfHR6ZMjuaZlQV8e1Gm7c1uzBjp7OljY5GLl7eXcNq9niUyNJjbZqdy3/yMK056sYQyDEsoY+NETSsbi1zsPF0/uH5kxdQknl1VwI3TJtv4iDE+0t+vfHaijp9tL2HnaWdXySCB5VOS+daCzG/cZM7GUMyY6ldlz5kGNhS5BmeahAQJ6xZk8OyqfCsbb4wfCAoSbpqZwk0zUzjkauaVHaW8c6CSL07V88Wpemalx/KtBRkszU+65rphftFCEZFE4LdAHs5eKw+oauMl1ywA/gWIBfqAH6vqb690b2uheF53bz+fHK9lY5FrsAR3THgIj1yfw5M35JERPzalIYwxV6e6uZOf7zzDb3aV0dLp7EaS7t6f5dZZqUSEBgdul5eI/A3QoKovisgPgQRV/bNLrpkOqKqeFJEMYC8wS1WbLndvSyie09zRw7uHqthSXEWTu3hdRlwET6/M58El2RNiz2xjxpPzXb28saecV78opbzhwofDu+am8T8eWBCwXV7rgDXux68BnwIXJRRVPTHkcaWI1AKTgcsmFHPtKps6eGu/i4+P1dLV6wy0z8mIZf3qAtYWptv+I8YEqOjwEJ5akc/jy/PYdrian20vYd/ZJn63t+Kq7ucvCSVVVavcj6uB1MtdLCJLgTDgtLcDm8iOVrWwscjFVyXnBhcirpkxmfWrClg+JcmrhemMMWMnOEi4qzCduwrT2VvWyMvbS/jpVdxnzBKKiHwIpA1z6oWhT1RVReQb++FEJB34JfCEqvZ/wzXrgfUAGVnZVx3zRNTXr+wqPcfGIhfHqp2NrMKCg7h/YQbPripgeurwM0KMMePDdbkJXJd7HT99bPTfO2YJRVVv/aZzIlIjIumqWuVOGLXfcF0ssAV4QVW/usxrvQS8BM4YyrVFPjF09vTx0bFa3t7voqq5E3A2snp0mbORlS1ENMZcib90eW0CngBedP/59qUXiEgYsBH4haq+ObbhjV+N7d1sKa5i66EqWt2zPbITI3lmRT7fW5xN9AhX1RpjjL+8W7wIvCEizwBlwAMAIrIY+IGqPus+thpIEpEn3d/3pKru90G8Aa+8oZ239rv45HgtPX1OI25+VhzrV0/hjjmpY1L3xxgzvvhFQlHVc8AtwxzfAzzrfvwr4FdjHNq4oqocqmxhY5FTPA5ABG6bncpzqwpYkpdgA+3GmKvmFwnFeFdfv7LzdD0bilyccu+IGB4SxHeuy+KZlflMmRzj4wiNMeOBJZRxrL27lw+P1vD2/kpq3bu7JUaH8diyXB5bnktyTLiPIzTGjCeWUMahc21dvFNcxXuHqjjf3QdAfnI0z9iOiMYYL7KEMo6cqT/Pxv0uPj9RR6+75O+SvASeW1XArbNSreKvMcarLKEEOFXlQEUzG4sq2HfWqUITJLC2MI1nVxWwKCfBxxEaYyYKSygBqq9f2XGqng1FFZQM2TzngcVZPL0yn9ykaB9HaIyZaCyhBJjOnj62Hanh7f2uwYH25Jgwnliex6PLckmIDvNxhMaYicoSSoBoau9mc3EVWw5W0dblrGgvSI7m2VW2ta4xxj9YQvFzlU0dbCxy8dGxmsEV7Qtz4nl+9RRum516zTusGWOMp1hC8VPHqlvYsO/i0vG3zkrl+RsLWJxrK9qNMf7HEoofGW6P9rDgIL61MJPnVuczNcVKxxtj/JclFD/Q09fPp+492ssbnW04J0WE8OiyXJ66wUrHG2MCgyUUH2rr6uW9Q9W8c6CShvZuANLjInhmZT4PLc0hxkrHG2MCiL1j+UB9Wxdv76/k/cPVdPQ4pVFmpk1i/eoC7p2fYXu0G2MCkiWUMXSm/jwbi1x8drKOPndplBumJLF+dQE3Tp9sA+3GmIBmCcXLVJWDrmY2FLnYW+bsQRIkcM+8dJ5fPYXCrDgfR2iMMZ5hCcVLhtuDJCI0iAcXZ/PMygJykqJ8HKExxniWJRQP6+zp46OjNby1v5Lqlk7A2YPkieV5PLY8l0QrjWKMGacsoXhIc0cPW4or2XKwipZOpzRKblIUz64q4Lu2B4kxZgKwhHKNqpo7eGt/JR8draGrtx+A+VlxPH/jFO6Yk2alUYwxE4YllKt0oqaVDUUuvjxdj3vCFjfPTGH96gKuz0+0GVvGmAnHLxKKiCQCvwXygDPAA6raeMk1ucBGIAgIBf5fVf3pWMapquw928iGfS4OupoBCA0Wvr0gk/WrC5ieaqVRjDETl18kFOCHwEeq+qKI/ND9/M8uuaYKWK6qXSISAxwSkU2qWunt4Hr6+tl+so4N+1yUNbQDEBMewvevz+GpFfmkxVlpFGOM8ZeEsg5Y4378GvAplyQUVe0e8jQcp6XiVe3dvbx/uJq391dy7rzz8qmx4Ty9Ip+Hr88hNiLU2yEYY0zA8JeEkqqqVe7H1UDqcBeJSDawBZgK/GdvtU7OtXXxTnEl7x6qpr3bKY0yLSWG9asLWLcgk7AQK41ijDGXGrOEIiIfAmnDnHph6BNVVRHRYa5DVcuBeSKSAbwlIm+qas0wr7UeWA+QkZU94hjLG9rZWOTik+O19LpH2q/PT+T5GwtYMz2FIJuxZYwx32jMEoqq3vpN50SkRkTSVbVKRNKB2ivcq1JEDgGrgDeHOf8S8BJA4YJFwyanIddypMrZzOrrMw3ueGBtYRrrV09hQXb8Ff9uxhhj/KfLaxPwBPCi+8+3L71ARLKAc6raISIJwErgH672Bfv6lV2l59iwz8XxmlYAwkOC+N7iLJ5dWUBecvTV3toYYyYkf0koLwJviMgzQBnwAICILAZ+oKrPArOAv3N3hwnwP1T14GhfqKu3j4+P1fJWkYvKZqc0SnxUKI8vz+Px5bkkx4R76K9kjDETi18kFFU9B9wyzPE9wLPuxx8A8672NVo7e9h6sIp3iqto7ugBIDsxkmdXFvC9xVlEhfnFj8IYYwLWuH8X7enr56XPT7PtyIXSKIWZcaxfXcBdc9MIsc2sjDHGI8Z9QimtP887xc6M5BunT+b51QUsn5JkpVGMMcbDxn1CEeDbCzN5bnUBs9JjfR2OMcaMW+M+ocxIm8TfP7jA12EYY8y4N+4HEEJtjMQYY8aEvdsaY4zxCEsoxhhjPMISijHGGI+whGKMMcYjLKEYY4zxCEsoxhhjPMISijHGGI8Q1ctuFxLwRKQVOO7rOEYgGaj3dRAjYHF6lsXpOYEQIwROnDNUddJovmHcr5QHjqvqYl8HcSUissfi9ByL07MCIc5AiBECK87Rfo91eRljjPEISyjGGGM8YiIklJd8HcAIWZyeZXF6ViDEGQgxwjiOc9wPyhtjjBkbE6GFYowxZgyM+4QiIv+HiKiIJLufi4j8k4icEpFiEVnk4/j+uzuO/SKyTUQy/DTOvxWRY+5YNopI/JBzf+6O87iI3OHDGL8nIodFpF9EFl9yzi9iHBLPne5YTonID30dzwAReVVEakXk0JBjiSLygYicdP+Z4MsY3TFli8gnInLE/Tv/T/4Yq4hEiMjXInLAHedfuo/ni8gu9+//tyIS5ss43TEFi0iRiGy+6hhVddx+AdnA+0AZkOw+thZ4F2czx2XALh/HGDvk8X8Efuqncd4OhLgf/zXw1+7Hs4EDQDiQD5wGgn0U4yxgBvApsHjIcb+J0R1PsDuGAiDMHdtsX/5+h8S2GlgEHBpy7G+AH7of/3Dgd+/jONOBRe7Hk4AT7t+zX8Xq/v8b434cCuxy/39+A3jIffynwB/4wc/0T4DfAJvdz0cd43hvofwD8H8CQweK1gG/UMdXQLyIpPskOkBVW4Y8jeZCrP4W5zZV7XU//QrIcj9eB7yuql2qWgqcApb6KMajqjrcIla/idFtKXBKVUtUtRt43R2jz6nq50DDJYfXAa+5H78G3D+mQQ1DVatUdZ/7cStwFMjEz2J1//9tcz8NdX8pcDPwpvu4z+MUkSzgbuBl93PhKmIctwlFRNYBLlU9cMmpTKB8yPMK9zGfEZEfi0g58H3gL9yH/S7OIZ7GaT2Bf8c5wN9i9Ld4riRVVavcj6uBVF8GcykRyQMW4nz697tY3V1J+4Fa4AOc1mnTkA9o/vD7/0ecD9/97udJXEWMAb1SXkQ+BNKGOfUC8F9wuml87nJxqurbqvoC8IKI/Dnwx8B/HdMA3a4Up/uaF4Be4NdjGduAkcRovEdVVUT8ZmqoiMQAvwf+N1VtcT5YO/wlVlXtAxa4xx03AjN9HNJFROQeoFZV94rImmu5V0AnFFW9dbjjIlKI01d+wP0PLAvYJyJLARfO2MqALPexMY9zGL8GtuIkFL+LU0SeBO4BblF3xypjHOcofpZDjfnP8gr8LZ4rqRGRdFWtcne71vo6IAARCcVJJr9W1Q3uw34ZK4CqNonIJ8BynC7sEHcLwNe//xXAfSKyFogAYoGfXE2M47LLS1UPqmqKquapah5Oc22RqlYDm4DH3bOolgHNQ5rIY05Epg15ug445n7sb3HeidMkvk9V24ec2gQ8JCLhIpIPTAO+9kWMl+FvMe4Gprln0YQBD7lj9FebgCfcj58AfN4SdPfxvwIcVdW/H3LKr2IVkckDMyJFJBK4DWe85xPgu+7LfBqnqv65qma53ysfAj5W1e9zNTH6embBWHwBZ7gwy0uAf8bpxzzIkNlAPort98AhoBh4B8j00zhP4fT773d//XTIuRfccR4H7vJhjN/C+fDQBdQA7/tbjEPiWYszM+k0TnedT+MZEte/AVVAj/tn+QxOf/pHwEngQyDRD+JciTO4XTzk3+Raf4sVmAcUueM8BPyF+3gBzoeaU8DvgHBf/0zdca3hwiyvUcdoK+WNMcZ4xLjs8jLGGDP2LKEYY4zxCEsoxhhjPMISijHGGI+whGKMMcYjLKEYY4zxCEsoxhhjPMISijFeICJ5ItLhLgro6XtHirN/Tre49/kxxh9YQjHGe06r6gJP31RVO9z3rfT0vY25FpZQjLlKIvKxu6WwX0Q6ReSBy1ybJ86Olz8XkRMi8msRuVVEvnDvLrh0NNcZ448soRhzlVT1ZndL4X/hFCX8/RW+ZSrwdzjly2cCj+DUpPpTnO0WRnudMX4loMvXG+NrIvI4cBfwHXX2vbicUlU96P6+w8BHqqoichDIu4rrjPErllCMuUoi8j2cXTbXqWrPCL6la8jj/iHP+7n4/+JIrzPGr9g/TmOugnuXuz8E7lHVTl/HY4w/sDEUY67Oazi72H3hHpR/xtcBGeNrth+KMV4gInk4GxXN9eJrnMHZeK3eW69hzGhYC8UY7+gD4ry5sBEIxRlXMcYvWAvFGGOMR1gLxRhjjEdYQjHGGOMRllCMMcZ4hCUUY4wxHmEJxRhjjEdYQjHGGOMRllCMMcZ4hCUUY4wxHvH/A5+Wytxf/FUvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5294673490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = np.linspace(-40,40,1000)\n",
    "q_in = -10 + 1j*z_R(0.05)       # at z=0, q = -z0 + 1j*z_R\n",
    "\n",
    "w_array = np.array( [w( propagate( FreeSpace(zi).T, q_in ) ) for zi in z] )\n",
    "R_array = np.array( [R( propagate( FreeSpace(zi).T, q_in ) ) for zi in z] )\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "ax.plot(z,w_array,z,-w_array,color='C0',lw=2)\n",
    "ax.fill_between(z,w_array,-w_array,lw=2,alpha=0.2)\n",
    "ax.set_xlim(z.min(),z.max())\n",
    "ax.set_ylabel(\"width, $w(z)$ [mm]\")\n",
    "ax.set_xlabel(\"$z$ [mm]\")\n",
    "\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cavity eigenmode calculation: Half-confocal cavity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:47: RuntimeWarning: divide by zero encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEOCAYAAABfM7oIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XucXWV97/HPd++5JJncr8SEkCgRRESwA2KlrQIqWI/oqbXeClJ7co4Va2vbUyyvV632crS1XnqOWlNFwSpIRY95CZYqWvWoSIIgkZuEkEhiIOR+nWRm9u/8sZ6Z2XPNnpm1MrP3fN+v136ttZ7nWc/z7GxYv1lrPetZigjMzMyKVJroDpiZWeNzsDEzs8I52JiZWeEcbMzMrHAONmZmVjgHGzMzK5yDjZmZFc7BxszMCudgY2ZmhXOwMTOzwjVNdAcmi4ULF8bKlSsnuhtmZnXjnnvu2RURi2op62CTrFy5kg0bNkx0N8zM6oakrbWW9WU0MzMrnM9sxuGHj+3irdevp6WplH3KpSHXW0fIa2kq0dqvbPmEZQbVXy73rpdLmuh/FjOzQRxsxuFYV4Xj3dmHYxPdm0y5pEGBrXWYIDVSQGwuZ5+Wconmsmhu6ttuKqsqry+/3/YweQ6GZlOTg804/MbqRay/7hI27TxMZ3eFzu4KXd2R1qM3baj1ruHKVPqndY1Qx1D53ZXgaKWbo53dE/3PM6SS6AtUTSWaSilwNfUPUll+tt1UKtHSVB3ASrT0lO0JjKXqgNhXtqksmko9AVKUS1nZpt68LL+5nNJK6tsnrTen9HJJSA6WZmPhYDMOpZKY1lxmZuvk+GeMCLor0T8wVYYKbsMHwq6q/bq7g65KCqCVLK+ratmZ8rsrkcpkZfvK9aV1p35UIjsjPNY1ec4GR6M3AJVKlMt9gapc6gtK1UGrOQW66vyeM7zqeoYLgOUU9Eol9VuWldU56CNRLmfLofYrKWt3pP2bSiVKJfovhQOtjcvkOEpaLpQOJE1lmE55orszpO5KXwDrCVw9gay7KoANld+VAmBvUKsq1zkgEPakVdJ6d9Wyb70yYHvwek9/uytBJeitu4PKRP9TnnTlEYLccIFMEiVBSVnA61kvS6g3PS2HKVsqDcgbrpzoX3bI/YbJG1APEmmBUFr2bdO7XVWuqiwD8wbUQb/trFypNELdVeUZVF9VXpXBfxto2PyznzGHlqZix4vlEmwkza+hWCUi9uXRntWv7CBUZpKcDI5KJQYGob5gVKnQd5Y3bOCqDJvXVZVfGZDX024lGLAddFcYsJ2Vq0RWT3e/JQO203KIeqvLVqrKdBMwOa/Q2jjcfd0lLJ41rdA28vpf/pfpM9J5dhlYkVN7ZiddSaJUFs2T86SxUJXewDNUwEvpAwJZJbJLuz3BL9J2d9V6BahUsu0KWZAL+u/bfzlSvdkyeoLtoDb6t9fTTqUSBP33yfKByPKy9Wyfvry03bsOkTZ61vvVEdm/ZU8d9Jbpq7N3/559Sf3p10ZU5fW1Ua3/1qDsfqa1lGguFf8UTF7B5qGIOG+kApLuzaktMzvJegKtNZ4zl86iuVx8sMmrhRflVMbMzBpQLsEmIjryKGNmZo0p13MnSe2SvirpJ5Lul7RR0v05t3GZpEckbZJ07RD575b0YGr/Tkmn5dm+mZmNXt5jgr4A/BmwEfIfGyqpDHwceBmwDVgvaV1EPFhV7F6gPSKOSHo78PfA7+TdFzMzq13ewebpiFiXc53VLgA2RcRmAEk3A1cAvcEmIr5TVf4u4C0F9sfMzGqQd7B5r6RPA3dS9Xx4RHwlp/qXAU9UbW8DXjhC+bcB38ipbTMzG6O8g83VwJlAM32X0QLIK9jUTNJbgHbgN0YoswZYA7BihR8BMjMrSt7B5vyIOCPnOqttB06t2l6e0vqRdClwHfAbETHsDFwRsRZYC9De3j7CY09mZjYeeT/J80NJZ+VcZ7X1wGpJqyS1AG8A+t0jknQe8Cng1RGxs8C+mJlZjfI+s7kQ+KmkzWT3bARERJyTR+UR0SXpGuAOsulvro+IByS9H9iQBif8AzAT+Lc0S+0vIuLVebRvZmZjk3ewuWyItFwvT0XE7cDtA9L+smr90jzbMzOz8cs72Cwiu1dy2oC6czmzMTOz+lRXD3WamVl9qreHOs3MrA7V20OdZmZWhxr2oU4zM5s86u2hTjMzq0P19lCnmZnVoSIe6rxP0uMU8FCnmZnVp5PxUKeZmU1xuQabiNiaZ31mZtYYcrlnI+kneZQxM7PGlNeZzXMk3T9CvoA5ObVlZmZ1Jq9gc2YNZbpzasvMzOpMLsHG92rMzGwkeT9nY2ZmNoiDjZmZFa6QYCOpTVK5iLrNzKz+5DX0uSTpTZJuk7QTeBjYIelBSf8g6fQ82jEzs/qU15nNd4BnAe8BTomIUyNiMXARcBfwQUlvyaktMzOrM3kNfb40IjolvY7sLZ0ARMQe4FbgVknNObVlZmZ1Jpczm4joTKufB75Yfb9G0tUDypiZ2RST9wCBh4Hv0v9M5p05t2FmZnUm72ATEfHPZG/mXCdpOtlUNWZmNoXlHWz2AkTEjcBngNuAGXk2IOkySY9I2iTp2iHyWyV9KeX/WNLKPNs3M7PRyzXYRMQlVetfBj4MLMir/nQv6OPA5cBZwBuHeDPo24C9EXE68BHgg3m1b2ZmY5PXczZDXiqLiK9HxMKRyozSBcCmiNgcEceBm4ErBpS5ArghrX8ZuCSnts3MbIxye85G0jslrahOlNQi6WJJNwBX5dDOMuCJqu1tKW3IMhHRBewnx7MrMzMbvbyes7kM+D3gJkmrgH3ANKAM/Afw0Yi4N6e2ciNpDbAGYMWKFScobWZmY5XXKwY6gE8An0hDnhcCRyNiXx71V9kOnFq1vTylDVVmm6Qmspe27R6m32uBtQDt7e2Rc1/NzCzJdYCApGkR0RkROwoINADrgdWSVklqAd4ArBtQZh19l+xeB3w7IhxIzMwmUN5Dn++W9I9FTbyZ7sFcA9wBPATcEhEPSHq/pFenYp8BFkjaBLwbGDQ82szMTq687tn0OBf4TeAjkkrAJ4Hb8jyziIjbgdsHpP1l1XoH8Nt5tWdmZuOX95nNXOAB4H1kswj8PbA55zbMzKzO5H1mswv4EfAD4CDZzfcDObdhZmZ1Ju8zm3bg58DzgAeBf4qI63Nuw8zM6kze09X8JCKuBt4CnA58T9Jf5NmGmZnVn1wvo0n6LtBG3+SbFbLhx3+XZztmZlZf8r5ncyXZ7AEH/GyLmZn1yDXYRMTWPOszM7PGkPcAATMzs0EcbMzMrHCFBhtJSyW1FtmGmZlNfkWf2XweeFjShwpux8zMJrG8R6P1ExGXprdkDnx1s5mZTSF5P2fTCvwWsLK67oh4f57tmJlZfcn7zOZrZK9hvgc4lnPdZmZWp/IONssj4rKc6zQzszqX9wCBH0p6Xs51mplZncvlzEbSRiBSfVdL2kx2GU1ARMQ5ebRjZmb1Ka/LaK/KqR4zM2tAuVxGi4itaV60P+hZr07Low0zM6tfed+zedkQaZfn3IaZmdWZvO7ZvJ3sDOaZku6vyppF9opoMzObwvK6Z/NF4BvA/wKurUo/GBF7cmrDzMzqVC7BJiL2kz3M+cY86jMzs8aS12W0d4+UHxEfzqGN+cCXyKbC2QK8PiL2DihzLvBJYDbQDfxtRHxpvG2bmdn45DVAYFb6tANvB5alz/8AXpBTG9cCd0bEauBO+l+u63EEuDIingtcBnxU0tyc2jczszHK6zLa+wAkfQ94QUQcTNt/BdyWRxvAFcBL0voNwH8Cfz6gHz+vWv+lpJ3AImBfTn0wM7MxyHvo8xLgeNX28ZSWS90RsSOtP3mieiVdALQAj+XUvpmZjVHeE3HeCNwt6atkU9W8huwspCaSvgWcMkTWddUbERGSYoR6lpK9uO2qiKiMUG4NsAZgxYoVtXbTzMxGKddgExF/K+kbwEUp6aqIuG8U+186XJ6kpyQtjYgdKZjsHKbcbLJLd9dFxF0naG8tsBagvb192OBlZmbjk9dotP8XERdJOkg2Iaeq8iIiZufQzDrgKuADafm1IfrRAnwVuDEivpxDm2ZmloO85ka7KC1nRcTstOz55BFoIAsyL5P0KHBp2kZSu6RPpzKvB34deKuk+9Ln3JzaNzOzMcr7tdCfB74HfD8iHs6z7ojYDVwyRPoG4PfT+r8C/5pnu2ZmNn55j0a7HlgK/G9JmyXdKuldObdhZmZ1Ju8BAt9Jz9qcD7yU7KHO5wIfy7MdMzOrL3lfRrsTaAN+BHwfOD8ihhw1ZmaTUyWC7kpULcmWQ6UNVTb6ykZk9VUCIi2z9JQGfWXTegRUGHnfSO3EoLzB+2TLnnp72ggqpLzK4PZI/Yk0RjVSPimt33ZvudQfespldRA9ZaIqL5Wraieq9+/ZHrB/9b4D969FDFGwpanErW//VRbMbK2tkjHK+zmb+4FfAc4mm5hzn6QfRcTRnNsxK1xE0FXJDqI9y2y9QqUCXZXKEHl9692VyrB5XVX5lcpw7WR51QfvEx7oe7ajZ18G7JsdXPvyB+9vU093rdFqHPK+jPbHAJJmAW8FPkv2kGaxIdMmrYi+A2dXJejqrqRl0Fmp0NVdldYvP1t2dmcH965+y5HzO3sO5N0DD+KVEQ78/YNJz8F3qmoqifKAT1NJlJSWpf7LcqlEuUS2FL37lJR9lNKybXrTSyWQRLkqXRLl0uD13nqUtdu7nvbNyg5e761HpP0GttfXN6V9SwKRtgEJ6Let3nSlsgyVV7UPg+rsX3/ffoP3Z2B9w/RvIA1OGqIUzJ/RMr7/YGqQ92W0a4BfIzu72UI2YOD7ebZhg3VXgs7udKDtrqRP1Xqlf3rXUGV611NdlaA77ds1zAF92CCRDvI99dSzppJoKovmUolyWTSVSjSXs7SmUinll/qVaypnB7rmlN5cLmUH66r83v161nvrTXVV1TvsQb/cc/AvUSpBU78D/jD7DRM0esulA7lZ3vK+jDYN+DBwT0R05Vz3pHOsq5sn93ewY//RIQ/cXSc4oA958K8Mc/Aftq7KpP8LPDvwZgfdlnKJ5nKJ5qZsu7lUtd6bn7absu2mknrXe/PKJVqa+rabyiVaqvKay30H7+YTHdirgkh1kCinv4LNbPzyvoz2oTzrm+x++Nhurv7s+onuBhK0pINvazootzRVfXq3y7SUU5mRyvUeyLMDcN96/4N9bwCpLpvymqrKlf2XstmUl/eZzZQys7WJxbNa+x2sW0c4yPcGg4EH+qr11n5lyoPKDLV/k/8CN7NJzsFmHM5fOZ+7rxt27lAzM0vynkHAzMxsEAcbMzMrnIZ6onQqkvQ0sHWMuy8EduXYnXow1b7zVPu+4O88VYznO58WEYtqKehgkwNJGyKifaL7cTJNte881b4v+DtPFSfrO/sympmZFc7BxszMCudgk4+1E92BCTDVvvNU+77g7zxVnJTv7Hs2ZmZWOJ/ZmJlZ4RxszMyscA42ZmZWOAcbMzMrnIONmZkVzsHGzMwK52BjZmaFc7AxM7PCOdiYmVnhHGzMzKxwkyrYSLpe0k5JP6tKmy/pm5IeTct5KV2S/knSJkn3S3pB1T5XpfKPSrpqIr6LmZn1mVTBBvgccNmAtGuBOyNiNXBn2ga4HFidPmuAT0IWnID3Ai8ELgDe2xOgzMxsYkyqYBMR3wP2DEi+Arghrd8AvKYq/cbI3AXMlbQUeAXwzYjYExF7gW8yOICZmdlJNKmCzTCWRMSOtP4ksCStLwOeqCq3LaUNl25mZhOkaaI7MBoREZJyeyeCpDVkl+Boa2v7lTPPPDOvqs3MGt4999yzKyIW1VK2HoLNU5KWRsSOdJlsZ0rfDpxaVW55StsOvGRA+n8OVXFErCW9OKi9vT02bNiQb8/NzBqYpK21lq2Hy2jrgJ4RZVcBX6tKvzKNSrsQ2J8ut90BvFzSvDQw4OUpzczMJsikOrORdBPZWclCSdvIRpV9ALhF0tuArcDrU/HbgVcCm4AjwNUAEbFH0l8D61O590fEwEEHZmZ2Evm10Ikvo5mZjY6keyKivZay9XAZzczM6pyDjZmZFc7BxszMCudgY2ZmhXOwMTOzwjnYmJlZ4RxszMyscCd8qDNN2X8ilYjYl0N/zMysAdUyg8Av00cjlCkDK3LpkZmZNZxags1DEXHeSAUk3ZtTf8zMrAHVcs/mRTmVMTOzKeqEwSYiOvIoY2ZmU1fNsz5LageuA05L+4nsfWbnFNQ3MzNrEKN5xcAXgD8DNgKVYrpjZmaNaDTB5umIWFdYT8zMrGGNJti8V9KngTuBYz2JEfGV3HtlZmYNZTTB5mrgTKCZvstoARQebCRtAQ4C3UBXRLSnh02/BKwEtgCvj4i9kgR8jOwtnkeAt0bET4ruo5mZDW80web8iDijsJ6c2EsjYlfV9rXAnRHxAUnXpu0/By4HVqfPC4FPpqWZmU2Q0cyN9kNJZxXWk9G7Arghrd8AvKYq/cbI3AXMlbR0IjpoZmaZ0ZzZXAj8VNJmsns2J3PocwD/ISmAT0XEWmBJROxI+U8CS9L6MuCJqn23pbQdmJnZhBhNsLlsiLTIqyMncFFEbJe0GPimpIf7dSIiUiAaFUlrgDUAK1Z4ajczs6KMJtgsov9DnT0KP7OJiO1puVPSV4ELgKckLY2IHeky2c5UfDtwatXuy1PaUPWuBdYCtLe3n6zAaWY25Yzmns0XgM8CvwX8l6pPoSS1SZrVsw68HPgZsA64KhW7CvhaWl8HXKnMhcD+qsttZmY2Aerhoc4lwFezEc00AV+MiH+XtB64RdLbgK3A61P528mGPW8iG/p89cnvspmZVZv0D3VGxGbg+UOk7wYuGSI9gHcU2SczMxudunio08zM6ls9PdRpZmZ1qp4f6jQzszox2oc675P0OCf/oU4zM6tj432o08zM7IRqDjYRsbXIjpiZWeM64T0bSSecnr+WMmZmNnXVcmbzHEn3j5AvYE5O/TEzswZUS7A5s4Yy3ePtiJmZNa4TBhvfqzEzs/EazXM2ZmZmY+JgY2ZmhRt1sElT/peL6IyZmTWmWoY+lyS9SdJtknYCDwM7JD0o6R8knV58N83MrJ7VcmbzHeBZwHuAUyLi1IhYDFwE3AV8UNJbCuyjmZnVuVqGPl8aEZ2SXgds7EmMiD3ArcCtkpqL6qCZmdW/E57ZRERnWv088MXq+zWSrh5QZtKQdJmkRyRtknTtRPfHzGwqG80AgYeB79L/TOad+Xdp/FJA/DhwOXAW8Ea/HsHMbOKMJthERPwz2Zs510maTjZVzWR0AbApIjZHxHHgZuCKCe6TmdmUNZpXDOwFiIgbJR0BbgNmFNKr8VsGPFG1vQ144QT1xcxsyhvNKwYuqVr/sqQO4HNFdOpkkbQGWAOwYsWKCe6NmVnjquU5myEvlUXE1yNi4UhlJtB24NSq7eUprZ+IWBsR7RHRvmjRopPWOTOzqaam52wkvVNSvz/9JbVIuljSDcBVxXRvzNYDqyWtktQCvAFYN8F9MjObsmq5jHYZ8HvATZJWAfuAaUAZ+A/goxFxb3FdHL2I6JJ0DXAHWT+vj4gHJrhbZmZTVi2vGOgAPgF8Ig15XggcjYh9RXduPCLiduD2ie6HmZmNbjRaz8ObOwrqi5mZNahaBgj8rqSnJW2TdFVKu1DS30i6p/gumplZvatlgMBfAq8EzgVWSfom8G9AC/BHBfbNzMwaRC2X0Q5FxHoASe8DngKePdnv2ZiZ2eRRS7A5JT38+Ej6bHOgMTOz0agl2LwXeB7w5rScJelbwL3AvRHxxQL7Z2ZmDaCWoc9rq7clLScLOueQzarsYGNmZiMa1dBngIjYRjax5Tfy746ZmTWi0bxiwMzMbEwcbMzMrHAONmZmVrgxBxtJSyW15tkZMzNrTOM5s/k88LCkD+XVGTMza0yjHo3WIyIuTS9NOyvH/piZWQOqZSLO0yW9eIj0FwPP9HtizMzsRGq5jPZR4MAQ6QdSXmEk/ZWk7ZLuS59XVuW9R9ImSY9IekVV+mUpbZOka4vsn5mZ1aaWy2hLImLjwMSI2ChpZe49GuwjEdHvvpCks8he9fxc4BnAtyQ9O2V/HHgZ2YOn6yWti4gHT0I/zcxsGLUEm7kj5E3PqyOjdAVwc0QcAx6XtAm4IOVtiojNAJJuTmUdbMzMJlAtwWaDpP8WEf9SnSjp94GT8fK0ayRdCWwA/iQi9gLLgLuqymxLaQBPDEh/YVEde/CXB/ib2x5kZmsTM6c1Zcu0Pqu1ibZ+2820tZZ716c1l8jGV5iZNb5ags0fAV+V9Gb6gks72cvTXjveDqQZpE8ZIus64JPAXwORlv8I/N5426xqew2wBmDFihWj3v+pgx388LHdY2q7XBJtLWVmTWvuDUhtrVmQ6g1Q05qYM72ZOdObmT2tmTkz0jKlOWCZWb2oZdbnp4BflfRS4OyUfFtEfDuPDkTEpbWUk/QvwNfT5nbg1Krs5SmNEdKHanstsBagvb09auxyr3OXz+Vf3/ZCDh3r5NCxbg51dHLoWBcHj3Vx+FgXhzq6su2OLg4f7799rKvCgY4uDnR0jbbZXs1lZYFoel8Qmj29mTnTm/oCVG9aM/NmtDC/rYW5M5qZ1lwec7tmZqM1mudsfgQ8ndY3FdCXQSQtjYgdafO1wM/S+jrgi5I+TDZAYDVwNyBgtaRVZEHmDcCbiurfvLYWLlq9cEz7dnZXOJwCz6FjVZ+OvuWBjk72H+3kwNFsuf9oJwc6unrTjnVV2HXoOLsOHR91+20tZea1tTBvRgvz2lqYP6N5wHYL89qamZ/S5s5oprXJAcrMxuaEwUZSE/B3ZJevtpId0E+V9FnguojoLLB/fy/pXLLLaFuA/w4QEQ9IuoXsxn8X8I6I6E79vQa4AygD10/W54CayyXmzmhh7oyWMdfR0dnNgaOdvUEpC0JdVes9AaqTfUeyz54jx9l7+DiHj3dz+PhRtu09WnN7M1ubmN/WwsKZLSyc2crCWa0s6l1maYtmtbJwZittrWN+XtjMGpAiRr56JOkjwCzgjyPiYEqbDXwIOBoR7yq8lydBe3t7bNiwYaK7cVJEBAePdbHvcF/w2XP4OHuPZJ89hzuztCPH2dezfeQ43ZXarzROby6zcFYKQCkgLZzZyuJZrZwyexqnzJnGktnTWNDWQqnk+05m9UjSPRHRXlPZGoLNo8CzY0BBSWXg4YhYPeaeTiJTKdiMRURwoKOL3YeOpUt3x7LPwWM8fegYTx/M0p4+mKUf66rUVG9zWSyelQWfU2ZnAWjpnGksSdunzJ7G4tmtvsdkNgmNJtjUcq0jBgaalNgtadQ31a0+SeodBffMRSOXjQgOHeti16HjvcGnJxDtPHCMJw908NSBDnbs72D/0U627zvK9n0jX86b39bCsrnTWTZ3OsvnTWfZvOksnzcjS5s3nTnTm3P8tmaWt1qCzYOSroyIG6sTJb0FeLiYblk9k8Ssac3MmtbMqoVtI5Y9erybpw509AtAT+7v6Evb38FTB4+xJ13q27h9/5D1zJrWlALRDJbPSwFp7nRWLJjBaQvamOl7SGYTqpbLaMuArwBH6f+czXTgtREx7NDieuLLaJNXdyXYdegY2/ZmZ0Db9h5he+/6UbbvPcrRzu4R61g0q5WVKfCsXDCDlQvbWLmgjdMWzGDWNJ8VmY1Frvdsqiq9mGwuMshGgX0X+J2I+MKYejnJONjUr4hgz+Hj/YLP9n1HeWLPEbbuOcIvdh/hePfw95AWtLVwWgpAz1o0k9MXZ5/T5s+gqeyX2ZoNJ9d7Nmnk2TvIpoNZB3wTuAb4NHA/0BDBxuqXJBbMbGXBzFbOWT54Kr/uSvDkgQ627jrM47sPs3X3EbbsSsvdh9l9+Di7Dx/nJ7/Y12+/5rJYtbAtBZ9Z2XLRTJ65qM0DFsxGqZbLaF8D9pI91HkJsJjsWZt3RcR9hffwJPGZzdRUqQQ7Dx5jy+7DPL7rMI/tPMSjOw+xaeehYQctlASnzp/BGUtm8ZylsznrGbM5a+lsls+b7umDbErJe+jzxoh4XlovAzuAFRHRMe6eTiIONjbQ4WNdbH76MI/uPMimFIA27TzE1j1HhnzmaFZrE2cuzQJQz+eMJbOY3uKzIGtMeQ997p0hIA133tZogcZsKG2tTTxv+Ryet3xOv/RjXd08vuswD+84yEM7DvDgjgM8tOMAuw4dZ/2Wvazfsre3bEmwevEszlk+h3NOncvzl8/hzFNm09Lke0E2tdRyZtMNHO7ZJBuFdiStR0TMLrSHJ4nPbGy8dh7s4KEUgHo+jz19eNBZUEu5xHOeMZvnL5/DOcvncu6pc3jmwpmeScHqTiGj0Rqdg40VoaOzmwd+eYD7t+3j/m37+em2fWx++vCgcrOnNfErp82jfeV8zl85n3OWz/EgBJv0HGzGwMHGTpb9Rzv52fYs8Nz/RLbcsb//lemWconnLZ9D+8p5nH/afNpXzhvXpK1mRXCwGQMHG5tI2/cdZcOWPazfsocNW/byyFMHqf5fU4Kzls7mxacv5MWnL+T8lfOY0eJZEWxiOdiMgYONTSb7j3Tyk1/sZcPWPazfspf7frGv34OpzWVx3op5XHT6Ql58+gLOWT6XZj+AaieZg80YONjYZNbR2c2GLXv5wWO7+MGmXWzcvr/fmc+s1iZ+7dkLeekZi3npmYtZOLN14jprU4aDzRg42Fg92X+kkx9t3s0PNu3iB4/t6jfoQIJzls/l4jMWc/GZi3nuM2Z7pJsVou6CjaTfBv4KeA5wQURsqMp7D/A2oBv4w4i4I6VfBnyM7I2cn46ID6T0VcDNwAKyiUN/NyJO+N5kBxurZ0/sOcJ3HtnJnQ/t5Eebd3O86n1Ci2e18vLnLuGVZy/lglXzPd+b5aYeg81zgArwKeBPe4KNpLOAm4ALgGcA3wKenXb7OfAyYBuwHnhjRDyYXhf9lYi4WdI/Az+NiE+eqA8ONtYojhzv4oebdvPtR3by7Yd28uSBvpFu89taeMVzl3D52Ut50bMW+D6PjUveMwgULiIeAoaaV+oK4OaIOAY8LmkTWeAB2BSK0dE+AAAGdUlEQVQRm9N+NwNXSHoIuBh4UypzA9kZ0wmDjVmjmNHSxKVnLeHSs5YQrwke+OUBvvGzHdy+8Uke33WYm+5+gpvufoI505t5xXOX8JrzlnHhqgW+1GaFmhTBZgTLgLuqtrelNIAnBqS/kOzS2b6I6BqivNmUI4mzl83h7GVz+NOXn8EjTx3k9o1P8o2NO3h05yFu2bCNWzZs4xlzpvGa85bxX1+wjNMXz5roblsDOmnBRtK3gFOGyLouIr52svpRTdIaYE3aPCTpkTFWtRDYlU+v6sZU+84N/X23kk3r/uf9kxv6Ow/D33l0Tqu14EkLNhFx6Rh22w6cWrW9PKUxTPpuYK6kpnR2U11+qD6tBdaOoV/9SNpQ63XLRjHVvvNU+77g7zxVnKzvPNnvDq4D3iCpNY0yWw3cTTYgYLWkVZJagDcA6yIb7fAd4HVp/6uACTlrMjOzPpMi2Eh6raRtwIuA2yTdARARDwC3kL2G+t+Bd0REdzpruQa4A3gIuCWVhexKwLvTYIIFwGdO7rcxM7OBJsXQ53onaU26JDdlTLXvPNW+L/g7TxUn6zs72JiZWeEmxWU0MzNrbA424yDpMkmPSNok6dqJ7k8RJJ0q6TuSHpT0gKR3pfT5kr4p6dG0nDfRfc2bpLKkeyV9PW2vkvTj9Ht/KQ1OaRiS5kr6sqSHJT0k6UWN/jtL+uP03/XPJN0kaVqj/c6Srpe0U9LPqtKG/F2V+af03e+X9IK8+uFgM0aSysDHgcuBs4A3pul1Gk0X8CcRcRZwIfCO9D2vBe6MiNXAnWm70byLbABKjw8CH4mI04G9ZHP2NZKPAf8eEWcCzyf77g37O0taBvwh0B4RZ5PNs/gGGu93/hxw2YC04X7Xy8lG/a4mewYxt9lXHGzG7gLSlDlpos+byabXaSgRsSMifpLWD5IdgJaRfdcbUrEbgNdMTA+LIWk58JvAp9O2yKZC+nIq0lDfWdIc4NdJozcj4nhE7KPBf2eyZw2nS2oCZgA7aLDfOSK+B+wZkDzc73oFcGNk7iJ7bnFpHv1wsBm7ZQyeMqehp8aRtBI4D/gxsCQidqSsJ4ElE9StonwU+J9kE8RC40+FtAp4GvhsunT4aUltNPDvHBHbgQ8BvyALMvvJZopv5N+5x3C/a2HHNQcbq4mkmcCtwB9FxIHqvPQwbcMMa5T0KmBnRNwz0X05iZqAFwCfjIjzgMMMuGTWgL/zPLK/5FeRzSrfxuDLTQ3vZP2uDjZjN9JUOg1FUjNZoPlCRHwlJT/Vc3qdljsnqn8FeDHwaklbyC6PXkx2P2NuutwCjfd7bwO2RcSP0/aXyYJPI//OlwKPR8TTEdEJfIXst2/k37nHcL9rYcc1B5uxG3LKnAnuU+7SvYrPAA9FxIerstaRTQcEDTYtUES8JyKWR8RKst/12xHxZhp4KqSIeBJ4QtIZKekSspk7GvZ3Jrt8dqGkGem/857v3LC/c5Xhftd1wJVpVNqFwP6qy23j4oc6x0HSK8mu7ZeB6yPibye4S7mTdBHwfWAjffcv/oLsvs0twAqySYNfHxEDb0LWPUkvIXuh36skPZPsTGc+cC/wlvSupYYg6VyyAREtwGbgarI/SBv2d5b0PuB3yEZd3gv8Ptk9iob5nSXdBLyEbHbnp4D3Av+XIX7XFHT/D9nlxCPA1dVvTh5XPxxszMysaL6MZmZmhXOwMTOzwjnYmJlZ4RxszMyscA42ZmZWOAcbMzMrnIONmZkVzsHGbIJIWinpqKT7Cqh7uqT7JB2XtDDv+s1Gy8HGbGI9FhHn5l1pRBxN9f4y77rNxsLBxqxAkr6dzjDuk9Qh6fUjlF2Z3pL5OUk/l/QFSZdK+kF6o+IFoylnNpk42JgVKCIuTmcYnyKb5PDWE+xyOvCPwJnp8ybgIuBPyeakG205s0mh6cRFzGw8JF1J9rrd34qI7hMUfzwiNqb9HiB7dW9I2gisHEM5s0nBwcasQJJ+G3gzcEV6Z8qJVM8uXKnartD//9day5lNCv6P0qwg6Y2ffwC8KiI6Jro/ZhPJ92zMinMD2ZsOf5AGCLxtojtkNlH8PhuzCSJpJfD1iDi7wDa2AO0RsauoNsxq4TMbs4nTDcwp8qFOoJm+N6yaTRif2ZiZWeF8ZmNmZoVzsDEzs8I52JiZWeEcbMzMrHAONmZmVjgHGzMzK5yDjZmZFc7BxszMCvf/AecmHFPGyHD3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f529458d210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenmode properties\n",
      "--------------------\n",
      "At input mirror of cavity:\n",
      "\t q_in = (-100+200j)\n",
      "\t 1/q_in = (-0.002-0.004j)\n",
      "\n",
      "Waist location = left mirror + 100.0 mm\n",
      "Waist radius = 0.2523 mm\n",
      "\n",
      "ROC at left mirror = -500.0 mm\n",
      "Beam radius at left mirror = 0.282094791774 mm\n",
      "\n",
      "ROC at right mirror =  inf mm\n",
      "Beam radius at right mirror = 0.252313252202 mm\n",
      "\n",
      "FSR = 1.5 GHz\n",
      "Transverse mode spacing = 221.4 MHz\n"
     ]
    }
   ],
   "source": [
    "L = 100             # spacer length\n",
    "R1 = 500            # ROC for left mirror\n",
    "R2 = float(\"inf\")   # ROC for right mirror\n",
    "        \n",
    "system = Cascade( [FreeSpace(L),\n",
    "                    ThinLens(R2/2),\n",
    "                    FreeSpace(L),\n",
    "                    ThinLens(R1/2)] )\n",
    "\n",
    "# print(eigenmode(system.T))\n",
    "q_stable = eigenmode(system.T)       # = -z0 + 1j*z_R\n",
    "\n",
    "if len(q_stable)==0 or q_stable[0].imag==0: print(\"Unstable cavity\")\n",
    "else: \n",
    "    q_in = q_stable[0]\n",
    "    z = np.linspace(0,L,5000)\n",
    "    \n",
    "    q_array = np.array( [propagate( FreeSpace(zi).T, q_in )  for zi in z] )\n",
    "    w_array = w(q_array)\n",
    "    \n",
    "    gp = np.arctan2(q_array.real,q_array.imag)\n",
    "    diff_guoy_phase = (gp[z==L] - gp[z==0])[0]\n",
    "    free_spectral_range = speed_of_light/(2*L)\n",
    "    \n",
    "    fig,axes = plt.subplots(nrows=2,ncols=1,sharex=True)\n",
    "    axes[0].plot(z,w_array,z,-w_array,color='C0',lw=2)\n",
    "    axes[0].fill_between(z,w_array,-w_array,lw=2,alpha=0.2)\n",
    "    axes[1].plot(z,R(q_array),lw=2)    \n",
    "    axes[1].set_ylim(-1000,1000)\n",
    "    axes[0].set_ylabel(\"width, $w(z)$ [mm]\")\n",
    "    axes[1].set_ylabel(\"ROC, $R(z)$ [mm]\")\n",
    "    axes[1].set_xlabel(\"$z$ [mm]\")\n",
    "    plt.show() \n",
    "    \n",
    "    print( \"Eigenmode properties\" )\n",
    "    print( ''.join([\"-\"]*20) )\n",
    "    print(\"At input mirror of cavity:\")\n",
    "    print( \"\\t q_in = {}\".format(np.round(q_in,2)) )\n",
    "    print( \"\\t 1/q_in = {}\".format(np.round(1/q_in,4)) )\n",
    "    print()\n",
    "    print( \"Waist location = left mirror + {:.4} mm\".format( z[w_array == w_array.min()][0] ) )\n",
    "    print( \"Waist radius = {:.4} mm\".format( w_array.min() ) )\n",
    "    print()\n",
    "    print( \"ROC at left mirror = {:4} mm\".format( R(q_array[0]) ) )\n",
    "    print( \"Beam radius at left mirror = {:4} mm\".format( w(q_array[0]) ) )\n",
    "    print()\n",
    "    print( \"ROC at right mirror = {:4} mm\".format( R(q_array[z==L])[0] ) )\n",
    "    print( \"Beam radius at right mirror = {:4} mm\".format( w(q_array[z==L])[0] ) )\n",
    "    print()\n",
    "    print( \"FSR = {:.4} GHz\".format(free_spectral_range/1e9))\n",
    "    print( \"Transverse mode spacing = {:.4} MHz\".format(diff_guoy_phase*free_spectral_range/1e6/pi))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cavity eigenmode calculation: 10 cm confocal cavity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEOCAYAAABfM7oIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHhRJREFUeJzt3X2cHmV97/HPl4Tw1CBgIKRJ1gRJjSgIdgWsnFYharC20WoREIkI5tSCR+vRCvI6Um1t8aGIngI1B4OBgkh5qHkJijwpHhUhEcozugZSEgMh8ujBBEJ+54+5ltx7733vztw7s/fu7Pf9eu1r77nmmpnf7AR+9zVzXdcoIjAzM6vSdt0OwMzM6s/JxszMKudkY2ZmlXOyMTOzyjnZmJlZ5ZxszMysck42ZmZWOScbMzOrnJONmZlVzsnGzMwqN7nbAYwV06ZNizlz5nQ7DDOzcWPVqlUbI2LPPHWdbJI5c+awcuXKbodhZjZuSFqTt65vo5mZWeXcshmBH/dt5MTlt724LJT9Fg1lDZ/Tisay/oWW9QrsZ1tdtSjbVjqwrP1xpAFRDtpnJ/vJfd4tjtfuOM3bDjx2/v2oVcU8csyannde9bwTsEfOPebeX456+c+h3Fnk88VW7O8RLy7HgGWa1reqs20fMXC5KYTc2zWtbywdLt5txx76PNpt1xjPzX/7JvaauiNVcrIZga0RbHp+a7fDMDMb85xsRuCPXj6NWz91BGt+8+yLZa2+NTQauL75w+BvPk2rG/YTg9e32abVt82hvuW03U/OeLcdt8W+h4m33Tfj5m95jUut/lbtvhEP+hY4RLwiXyunaGNoyH3lrVfmQcl3Dvljy33U0moV/XMMakU3tbTV5g5B47KaNtq2bfM+1bTdMLG02mdTEM2t97axNMfeolX/ir2nMnXH6lOBk80ITNpO7DRlEjtuP6nboZiZdWSnKZNK//LSijsImJlZ5Qq1bCTtkaPa1oh4ssN4zMyshoreRvt1+hmqzTUJ6Ok4IjMzq52iyea+iDhoqAqSbh9BPGZmVkNFn9m8vqQ6ZmY2gRRKNhGxqYw6ZmY2sXTU9VlSL3A68LK0DwEREQeUGJuZmdVEp+NsLgY+AdwFeAi9mZkNqdNk81hErCg1EjMzq61Ok80Zks4HbgA29xdGxJWlRGVmZrXSabI5AZgPbM+222gBONmYmdkgnSab10XEK0qNpImkhcBXyAaJnh8RZzat/xhwErAFeAz4QESsSeteIHueBPBfEfHnVcZqZmZD63RutJ9I2q/USBpImgScAxwJ7Acc0+J4twO9qQfc5cAXGtb9LiIOTD9ONGZmXdZpy+ZQ4D8lrSZ7ZlN21+eDgb6IWA0g6VJgEXBvf4WIuKmh/i3AcSUd28zMStZpslnYoqzMV/TNBB5uWF4LHDJE/ROB7zYs7yhpJdkttjMj4j9KjM3MzArqNNnsycBBnf1GfVCnpOOAXuBPGopfFhHrJO0D3Cjproj4VYttlwBLAHp6PHeomVlVxuqgznXA7IblWalsAEkLyJLen0REYxfsden3akk/AA4CBiWbiFgKLAXo7e0t9+XpZmb2orE6qPM2YJ6kuWRJ5mjg2MYKkg4CvgYsjIgNDeW7A89GxGZJ04A3MLDzgJmZjbIxOagzIrZIOgW4lqzr87KIuEfSZ4GVKdF9Efg94N/TK037uzi/EviapK1kve3OjIh7Wx7IzMxGxZgd1BkR1wDXNJV9uuHzgjbb/QTYv6w4zMxs5MbsoE4zM6uPMTmo08zM6mUkgzrvkPQg1QzqNDOzGilzUKeZmVlLHSWb/gkvzczM8ij0zEbSz8uoY2ZmE0vRls0rJd05xHoBLxlBPGZmVkNFk838HHVe6CQQMzOrr0LJxs9qzMysE52OszEzM8vNycbMzCo3omQjaZf0CmczM7O2inZ93k7SsZKulrQBuB9YL+leSV+UtG81YZqZ2XhWtGVzE/By4DRg74iYHRF7AYcBtwCfT2/ONDMze1HRrs8LIuJ5Se8me0snABHxOHAFcIWk7csM0MzMxr9CLZuIeD59vAi4pPF5jaQTmuqYmZkBnXcQuB/4IQNbMh8uJyQzM6ubTpNNRMS/kr2Zc4WkncimqjEzMxuk02TzBEBEXAh8Hbga2LmsoAAkLZT0gKQ+Sae2WL+DpG+l9T+TNKdh3Wmp/AFJby0zLjMzK66jZBMRRzR8vhw4C3hpWUGlZ0HnAEcC+wHHtHgz6InAExGxL/Bl4PNp2/2Ao4FXkb1351yPBTIz666i42xa3iqLiO9ExLSh6hR0MNAXEasj4jngUmBRU51FwPL0+XLgiHTsRcClEbE5Ih4E+tL+zMysSwqPs5H0YUk9jYWSpkg6XNJyYHEJcc0EHm5YXpvKWtaJiC3AU2StqzzbmpnZKCo6zmYh8AHgm5LmAk8COwKTgO8DZ0fE7eWGWB1JS4AlAD09PcPUNjOzThV9xcAm4Fyy5yDbA9OA30XEkyXHtQ6Y3bA8K5W1qrNW0mSyl7b9Jue2AETEUmApQG9vb5QSuZmZDdJRBwFJO0bE8xGxvoJEA3AbME/SXElTyB74r2iqs4Jtt+zeDdwYEZHKj0691eYC84BbK4jRzMxyKnobrd+tkq4DzouIvjIDguwZjKRTgGvJbtEti4h7JH0WWBkRK8i6XF8kqQ94nCwhkepdBtwLbAFOjgi/PdTMrIuUNQYKbiRtB/wp2fOO7YDzgKujk52NEb29vbFy5crC2z2z6Xke2vhsBRGZmVVv/oypbD+psyGXklZFRG+eup0O6twNuAf4DNksAl8AVne4LzMzq7lOb6NtBH4K/Bh4huwh+9NlBWVmZvXSacumF/gFsD/Zs5GvRsSy0qIyM7Na6XS6mp9HxAnAccC+wM2SPlVqZGZmVhsd3UaT9ENgF7ZNvrmVrPvxP5YUl5mZ1Uinz2yOJ5s94Onx3APNzMxGR0fJJiLWlB2ImZnVV6cdBMzMzHJzsjEzs8qVkmwkzZC0Qxn7MjOz+imrZXMRcL+kL5W0PzMzq5FOe6MNEBEL0lsym1/dbGZm1vE4mx2AdwFzGvcREZ8tJywzM6uTTls23yZ7DfMqYHN54ZiZWR11mmxmRcTCUiMxM7Pa6rSDwE8k7V9qJGZmVluFWjaS7gIibXeCpNVkt9EEREQcUH6IZmY23hW9jfb2SqIwM7NaK3QbLSLWpHnR/rr/c2NZGQFJ2kPSdZJ+mX7v3qLOgZJ+KukeSXdKek/Dum9IelDSHennwDLiMjOzznX6zObNLcqOHEkgDU4FboiIecANabnZs8DxEfEqYCFwtqTdGtZ/IiIOTD93lBSXmZl1qOgzmw+RtWD2kXRnw6qpZK+ILsMi4I3p83LgB8AnGytExC8aPv9a0gZgT7LXHpiZ2RhT9JnNJcB3gX9iYIvjmYh4vKSYpkfE+vT5EWD6UJUlHQxMAX7VUPw5SZ8mtYwiwmOBzMy6qFCyiYinyAZzHjOSg0q6Hti7xarTm44Xktq+nE3SDLJ52RZHxNZUfBpZkpoCLCVrFbWc2UDSEmAJQE9PT8GzMDOzvIreRvvYUOsj4qw8+4mIBUMc41FJMyJifUomG9rU2xW4Gjg9Im5p2Hd/q2izpAuAjw8Rx1KyhERvb6/fOGpmVpGiHQSmpp9e4EPAzPTzV8BrS4ppBbA4fV5MNjXOAJKmAFcBF0bE5U3rZqTfAt4B3F1SXGZm1qGit9E+AyDpZuC1EfFMWv47slZGGc4ELpN0IrAGOCodoxf4q4g4KZX9MfBSSe9P270/9Ty7WNKeZANN7yBLhGZm1kWdzo02HXiuYfk5hnmQn1dE/AY4okX5SuCk9PnfgH9rs/3hZcRhZmbl6TTZXAjcKukqshbEO8i6KZuZmQ3SUbKJiM9J+i5wWCpa7MGTZmbWTtHeaP83Ig6T9AzZhJxqWBcRsWvZAZqZ2fhXtIPAYen31GrCMTOzOupobjRJF0n6oKT5ZQdkZmb10+lEnMuAGcD/lrRa0hWSPlJiXGZmViOddhC4KY21eR3wJrKxLK8CvlJibGZmVhMdJRtJNwC7AD8FfgS8LiJaTitjZmbW6W20O8kGcr4aOAB4taSdSovKzMxqpdPbaH8DIGkq8H7gArJZnHcoLTIzM6uNTm+jnQL8N+APgYfIOgz8qLywzMysTjqdrmZH4CxgVURsKTGecWXT8y/wyFOb2PjbzUR6QUHQ8KaCFh8HvMegxTbR4kUHMWA/McR+GrcZvKNoGdrgYw/cT6t4Wm09eJsif4vc8bY4dpG/RQz60Lifzt8yoW3jm4epV2inpVfNJkMv+fAVxFns+MXPqXmT/r+Lmur1f+i/voO2a9q+OaTBx1Ob8tZx0LyfdnE0xdl62+by7MNLdp7MzN12ZtJ2hf51FtbpbbQvlR3IeHTbQ4/zvq/f2u0wzMxG5NbTj2CvqTtWeoxOWzYGTJm0HXvvuu0CNX+bycoGf1toLGr1jWS4/TR/Q2n82LLeMMcZGFvTN6sO42W4/QwR78B67b+pDRtPm/3kPe+88raHWrXeypB3t0Vabrn3WeCUKjl+/sMPupPQ3JLvvz7blpuO0XZ90/7a7JfhtmsTR3NrvX2cLe5SNP2B2m27XYHWYaecbEbgkH1eyi2fGvQ2BDMza9Jp12czM7PcnGzMzKxyquo+8ngj6TGy11B3YhqwscRwxoOJds4T7XzB5zxRjOScXxYRe+ap6GRTAkkrI6K323GMpol2zhPtfMHnPFGM1jn7NpqZmVXOycbMzCrnZFOOpd0OoAsm2jlPtPMFn/NEMSrn7Gc2ZmZWObdszMysck42ZmZWOScbMzOrnJONmZlVzsnGzMwq52RjZmaVc7IxM7PKOdmYmVnlnGzMzKxyTjZmZla5MZVsJC2TtEHS3Q1le0i6TtIv0+/dU7kkfVVSn6Q7Jb22YZvFqf4vJS3uxrmYmdk2YyrZAN8AFjaVnQrcEBHzgBvSMsCRwLz0swQ4D7LkBJwBHAIcDJzRn6DMzKw7xlSyiYibgcebihcBy9Pn5cA7GsovjMwtwG6SZgBvBa6LiMcj4gngOgYnMDMzG0VjKtm0MT0i1qfPjwDT0+eZwMMN9damsnblZmbWJZO7HUARERGSSnsngqQlZLfg2GWXXf5w/vz5Ze3abMS2bA3uW/80k7cTr5yxa7fDMRtk1apVGyNizzx1x0OyeVTSjIhYn26TbUjl64DZDfVmpbJ1wBubyn/QascRsZT04qDe3t5YuXJluZGbjcDG326m9x+u56W7TGHl/3pzt8MxG0TSmrx1x8NttBVAf4+yxcC3G8qPT73SDgWeSrfbrgXeImn31DHgLanMzMy6ZEy1bCR9k6xVMk3SWrJeZWcCl0k6EVgDHJWqXwO8DegDngVOAIiIxyX9PXBbqvfZiGjudGBmZqNoTCWbiDimzaojWtQN4OQ2+1kGLCsxNDMzG4HxcBvNzMzGOScbMzOrnJONmZlVzsnGzMwq52RjZmaVc7IxM7PKOdmYmVnlhh1nk6bsH87WiHiyhHjMzKyG8gzq/HX60RB1JgE9pURkZma1kyfZ3BcRBw1VQdLtJcVjZmY1lOeZzetLqmNmZhPUsMkmIjaVUcfMzCau3BNxSuoFTgdelrYT2XyYB1QUm5mZ1USRWZ8vBj4B3AVsrSYcMzOroyLJ5rGIWFFZJGZmVltFks0Zks4HbgA29xdGxJWlR2VmZrVSJNmcAMwHtmfbbbQAKk82kh4CngFeALZERG8abPotYA7wEHBURDwhScBXyN7i+Szw/oj4edUxmplZe0WSzesi4hWVRTK8N0XExoblU4EbIuJMSaem5U8CRwLz0s8hwHnpt5mZdUmRudF+Imm/yiIpbhGwPH1eDryjofzCyNwC7CZpRjcCNDOzTJGWzaHAf0paTfbMZjS7PgfwfUkBfC0ilgLTI2J9Wv8IMD19ngk83LDt2lS2HjMz64oiyWZhi7IoK5BhHBYR6yTtBVwn6f4BQURESkSFSFoCLAHo6fHUbmZmVSmSbPZk4KDOfpW3bCJiXfq9QdJVwMHAo5JmRMT6dJtsQ6q+DpjdsPmsVNZqv0uBpQC9vb2jlTjNzCacIs9sLgYuAN4F/FnDT6Uk7SJpav9n4C3A3cAKYHGqthj4dvq8AjhemUOBpxput5mZWReMh0Gd04Grsh7NTAYuiYjvSboNuEzSicAa4KhU/xqybs99ZF2fTxj9kM3MrNGYH9QZEauB17Qo/w1wRIvyAE6uMiYzMytmXAzqNDOz8W08Deo0M7NxajwP6jQzs3Gi6KDOOyQ9yOgP6jQzs3FspIM6zczMhpU72UTEmioDMTOz+hr2mY2kYafnz1PHzMwmrjwtm1dKunOI9QJeUlI8ZmZWQ3mSzfwcdV4YaSBmZlZfwyYbP6sxM7ORKjLOxszMrCNONmZmVrnCySZN+T+pimDMzKye8nR93k7SsZKulrQBuB9YL+leSV+UtG/1YZqZ2XiWp2VzE/By4DRg74iYHRF7AYcBtwCfl3RchTGamdk4l6fr84KIeF7Su4G7+gsj4nHgCuAKSdtXFaCZmY1/w7ZsIuL59PEi4JLG5zWSTmiqM2ZIWijpAUl9kk7tdjxmZhNZkQ4C9wM/ZGBL5sPlhzRyKSGeAxwJ7Acc49cjmJl1T5FkExHxr2Rv5lwhaSeyqWrGooOBvohYHRHPAZcCi7ock5nZhFXkFQNPAETEhZKeBa4Gdq4kqpGbCTzcsLwWOKRLsZiZTXhFXjFwRMPnyyVtAr5RRVCjRdISYAlAT09Pl6MxM6uvPONsWt4qi4jvRMS0oep00TpgdsPyrFQ2QEQsjYjeiOjdc889Ry04M7OJJtc4G0kfljTgq7+kKZIOl7QcWFxNeB27DZgnaa6kKcDRwIoux2RmNmHluY22EPgA8E1Jc4EngR2BScD3gbMj4vbqQiwuIrZIOgW4lizOZRFxT5fDMjObsPK8YmATcC5wburyPA34XUQ8WXVwIxER1wDXdDsOMzMr1hutf/Dm+opiMTOzmsrTQeB9kh6TtFbS4lR2qKR/kLSq+hDNzGy8y9NB4NPA24ADgbmSrgP+HZgCfLTC2MzMrCby3Eb7bUTcBiDpM8CjwB+M9Wc2ZmY2duRJNnunwY8PpJ+1TjRmZlZEnmRzBrA/8N70e6qk64Hbgdsj4pIK4zMzsxrI0/V5aeOypFlkSecAslmVnWzMzGxIhbo+A0TEWrKJLb9bfjhmZlZHRV4xYGZm1hEnGzMzq5yTjZmZVa7jZCNphqQdygzGzMzqaSQtm4uA+yV9qaxgzMysngr3RusXEQvSS9P2KzEeMzOroTwTce4r6Q0tyt8A7OP3xJiZ2XDy3EY7G3i6RfnTaV1lJP2dpHWS7kg/b2tYd5qkPkkPSHprQ/nCVNYn6dQq4zMzs3zy3EabHhF3NRdGxF2S5pQe0WBfjogBz4Uk7Uf2qudXAb8PXC/pD9Lqc4A3kw08vU3Sioi4dxTiNDOzNvIkm92GWLdTWYEUtAi4NCI2Aw9K6gMOTuv6ImI1gKRLU10nGzOzLspzG22lpA82F0o6CRiNl6edIulOScsk7Z7KZgIPN9RZm8ralZuZWRfladl8FLhK0nvZllx6yV6e9s6RBpBmkN67xarTgfOAvwci/f5n4AMjPWbDsZcASwB6enrK2q2ZmTXJM+vzo8AfSXoT8OpUfHVE3FhGABGxIE89Sf8H+E5aXAfMblg9K5UxRHmrYy8FlgL09vZGzpDNzKygIuNsfgo8lj73VRDLIJJmRMT6tPhO4O70eQVwiaSzyDoIzANuBQTMkzSXLMkcDRw7GrGamVl7wyYbSZOBfyS7fbWG7H/osyVdAJweEc9XGN8XJB1IdhvtIeC/A0TEPZIuI3vwvwU4OSJeSPGeAlwLTAKWeRyQmVn35WnZfBGYCsyNiGcAJO0KfCn9fKSq4CLifUOs+xzwuRbl1wDXVBWTmZkVl6c32tuBD/YnGoCIeBr4EPC2tluZmZkleZJNRMSgh+fptpUfqpuZ2bDyJJt7JR3fXCjpOOD+8kMyM7O6yfPM5mTgSkkfYOA4m50oYZyNmZnVX55xNuuAQyQdTjYXGWQP4H8IvAe4uLrwzMysDvK8YmBXSacBfwE8APwL2TtsfgkcVW14ZmZWB3luo10EPEE2qPMk4FNkY23eGRF3VBibmZnVRJ5ks09E7A8g6XxgPdATEZsqjczMzGojT2+0F2cISN2d1zrRmJlZEXlaNq+R1P+mTgE7pWWRjcHZtbLozMysFvL0Rps0GoGYmVl95bmNZmZmNiJONmZmVjknGzMzq5yTjZmZVc7JxszMKjcmko2kv5R0j6Stknqb1p0mqU/SA5Le2lC+MJX1STq1oXyupJ+l8m9JmjKa52JmZoONiWQD3E0299rNjYWS9gOOJpsAdCFwrqRJkiYB5wBHks3TdkyqC/B54MsRsS/ZNDsnjs4pmJlZO2Mi2UTEfRHxQItVi4BLI2JzRDwI9AEHp5++iFgdEc8BlwKLJAk4HLg8bb8ceEf1Z2BmZkMZE8lmCDOBhxuW16ayduUvBZ6MiC1N5WZm1kV5pqsphaTrgb1brDo9Ir49WnE0krQEWJIWfyupVesqj2nAxnKiGjcm2jl37XzXAPp0N4484a4x+JyLelneiqOWbCJiQQebrQNmNyzPSmW0Kf8NsJukyal101i/VUxLgaUdxDWApJUR0Tt8zfqYaOc80c4XfM4TxWid81i/jbYCOFrSDpLmAvOAW4HbgHmp59kUsk4EKyIigJuAd6ftFwNdaTWZmdk2YyLZSHqnpLXA64GrJV0LEBH3AJcB9wLfA06OiBdSq+UU4FrgPuCyVBfgk8DHJPWRPcP5+uiejZmZNVPWGLCRkLQk3ZKbMCbaOU+08wWf80QxWufsZGNmZpUbE7fRzMys3pxsRqDdlDl1Imm2pJsk3ZumFPpIKt9D0nWSfpl+797tWMuWZqu4XdJ30nKtp0KStJukyyXdL+k+Sa+v+3WW9Dfp3/Xdkr4pace6XWdJyyRtkHR3Q1nL66rMV9O53ynptWXF4WTToWGmzKmTLcD/jIj9gEOBk9N5ngrcEBHzgBvSct18hKwDSr+6T4X0FeB7ETEfeA3Zudf2OkuaCfwPoDciXg1MIuvZWrfr/A2y6b4atbuuR5L1+p1HNgbxvLKCcLLpXMspc7ocU+kiYn1E/Dx9fobsf0Azyc51eapWu2mBJM0C/hQ4Py3XeiokSS8B/pjUezMinouIJ6n5dSYba7iTpMnAzsB6anadI+Jm4PGm4nbXdRFwYWRuIRu3OKOMOJxsOtduypzakjQHOAj4GTA9ItanVY8A07sUVlXOBv4W2JqW6z4V0lzgMeCCdOvwfEm7UOPrHBHrgC8B/0WWZJ4CVlHv69yv3XWt7P9rTjaWi6TfA64APhoRTzeuS4Npa9OtUdLbgQ0RsarbsYyiycBrgfMi4iDg/9F0y6yG13l3sm/yc4HfB3Zh8O2m2hut6+pk07mhptKpFUnbkyWaiyPiylT8aH/zOv3e0K34KvAG4M8lPUR2e/RwsucZu6XbLVC/670WWBsRP0vLl5Mlnzpf5wXAgxHxWEQ8D1xJdu3rfJ37tbuulf1/zcmmcy2nzOlyTKVLzyq+DtwXEWc1rFpBNh0Q1GxaoIg4LSJmRcQcsut6Y0S8lxpPhRQRjwAPS3pFKjqCbOaO2l5nsttnh0raOf077z/n2l7nBu2u6wrg+NQr7VDgqYbbbSPiQZ0jIOltZPf2JwHLIuJzXQ6pdJIOA34E3MW25xefIntucxnQQzYx8VER0fwQctyT9Ebg4xHxdkn7kLV09gBuB46LiM3djK9Mkg4k6xAxBVgNnED2hbS211nSZ4D3kPW6vB04iewZRW2us6RvAm8km935UeAM4D9ocV1T0v0XstuJzwInRMTKUuJwsjEzs6r5NpqZmVXOycbMzCrnZGNmZpVzsjEzs8o52ZiZWeWcbMzMrHJONmZmVjknG7MukTRH0u8k3VHBvneSdIek5yRNK3v/ZkU52Zh1168i4sCydxoRv0v7/XXZ+zbrhJONWYUk3ZhaGHdI2iTpqCHqzklvyfyGpF9IuljSAkk/Tm9UPLhIPbOxxMnGrEIRcXhqYXyNbJLDK4bZZF/gn4H56edY4DDg42Rz0hWtZzYmTB6+ipmNhKTjyV63+66IeGGY6g9GxF1pu3vIXt0bku4C5nRQz2xMcLIxq5CkvwTeCyxK70wZTuPswlsblrcy8L/XvPXMxgT/ozSrSHrj518Db4+ITd2Ox6yb/MzGrDrLyd50+OPUQeDEbgdk1i1+n41Zl0iaA3wnIl5d4TEeAnojYmNVxzDLwy0bs+55AXhJlYM6ge3Z9oZVs65xy8bMzCrnlo2ZmVXOycbMzCrnZGNmZpVzsjEzs8o52ZiZWeWcbMzMrHJONmZmVjknGzMzq9z/B4P9C+MJZdSYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f525a068090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenmode properties\n",
      "--------------------\n",
      "At input mirror of cavity:\n",
      "\t q_in = (-50+545.44j)\n",
      "\t 1/q_in = (-0.0002-0.0018j)\n",
      "\n",
      "Waist location = left mirror + 49.9 mm\n",
      "Waist radius = 0.4167 mm\n",
      "\n",
      "ROC at left mirror = -6000.0 mm\n",
      "Beam radius at left mirror = 0.418421458045 mm\n",
      "\n",
      "ROC at right mirror = 6000.0 mm\n",
      "Beam radius at right mirror = 0.418421458045 mm\n",
      "\n",
      "FSR = 1.5 GHz\n",
      "Transverse mode spacing = 87.29 MHz\n"
     ]
    }
   ],
   "source": [
    "L = 100      # spacer length\n",
    "R1 = 6000    # ROC for left mirror\n",
    "R2 = 6000    # ROC for right mirror\n",
    "        \n",
    "system = Cascade( [FreeSpace(L),\n",
    "                    ThinLens(R2/2),\n",
    "                    FreeSpace(L),\n",
    "                    ThinLens(R1/2)] )\n",
    "\n",
    "# print(eigenmode(system.T))\n",
    "q_stable = eigenmode(system.T)       # = -z0 + 1j*z_R\n",
    "\n",
    "if len(q_stable)==0 or q_stable[0].imag==0: print(\"Unstable cavity\")\n",
    "else: \n",
    "    q_in = q_stable[0]\n",
    "    z = np.linspace(0,L,500)\n",
    "    \n",
    "    q_array = np.array( [propagate( FreeSpace(zi).T, q_in )  for zi in z] )\n",
    "    w_array = w(q_array)\n",
    "    \n",
    "    gp = np.arctan2(q_array.real,q_array.imag)\n",
    "    diff_guoy_phase = (gp[z==L] - gp[z==0])[0]\n",
    "    free_spectral_range = speed_of_light/(2*L)\n",
    "    \n",
    "    fig,axes = plt.subplots(nrows=2,ncols=1,sharex=True)\n",
    "    axes[0].plot(z,w_array,z,-w_array,color='C0',lw=2)\n",
    "    axes[0].fill_between(z,w_array,-w_array,lw=2,alpha=0.2)\n",
    "    axes[1].plot(z,R(q_array),lw=2)    \n",
    "    axes[1].set_ylim(-1000,1000)\n",
    "    axes[0].set_ylabel(\"width, $w(z)$ [mm]\")\n",
    "    axes[1].set_ylabel(\"ROC, $R(z)$ [mm]\")\n",
    "    axes[1].set_xlabel(\"$z$ [mm]\")\n",
    "    plt.show() \n",
    "    \n",
    "    print( \"Eigenmode properties\" )\n",
    "    print( ''.join([\"-\"]*20) )\n",
    "    print(\"At input mirror of cavity:\")\n",
    "    print( \"\\t q_in = {}\".format(np.round(q_in,2)) )\n",
    "    print( \"\\t 1/q_in = {}\".format(np.round(1/q_in,4)) )\n",
    "    print()\n",
    "    print( \"Waist location = left mirror + {:.4} mm\".format( z[w_array == w_array.min()][0] ) )\n",
    "    print( \"Waist radius = {:.4} mm\".format( w_array.min() ) )\n",
    "    print()\n",
    "    print( \"ROC at left mirror = {:4} mm\".format( R(q_array[0]) ) )\n",
    "    print( \"Beam radius at left mirror = {:4} mm\".format( w(q_array[0]) ) )\n",
    "    print()\n",
    "    print( \"ROC at right mirror = {:4} mm\".format( R(q_array[z==L])[0] ) )\n",
    "    print( \"Beam radius at right mirror = {:4} mm\".format( w(q_array[z==L])[0] ) )\n",
    "    print()\n",
    "    print( \"FSR = {:.4} GHz\".format(free_spectral_range/1e9))\n",
    "    print( \"Transverse mode spacing = {:.4} MHz\".format(diff_guoy_phase*free_spectral_range/1e6/pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modematching to confocal cavity: Nagourney's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "laser beam radius = 0.52 mm\n",
      "100.0\n",
      "0.282094791774\n",
      "focal length ratio = 1.844\n",
      "optimized q @ cavity = (-54.665+10.008j)\n"
     ]
    }
   ],
   "source": [
    "lamda = 1e-3  # mm\n",
    "q_laser = 200 + 1j*50\n",
    "w_laser = w(q_laser)\n",
    "q_target = -50 + 50j\n",
    "print( \"laser beam radius = {} mm\".format(np.round(w_laser,3)) )\n",
    "d0,f_ratio = mode_matching_nagourney(w_laser,q_target)\n",
    "print(\"focal length ratio = {}\".format(np.round(f_ratio,3)))\n",
    "\n",
    "# test\n",
    "mm = make_mode_matching_system([223,100,50])\n",
    "print( \"optimized q @ cavity = {}\".format(np.round(propagate(mm.T,q_laser),3) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automated modematching"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 9.4185167729181974e-09\n",
      "     jac: array([ -2.32747211e-04,   2.22872237e-04,   2.12350915e-07,\n",
      "         0.00000000e+00])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 145\n",
      "     nit: 29\n",
      "    njev: 28\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([ 152.36149265,  150.80716567,   54.19587903])\n",
      "Target q @ cavity = (-50+50j)\n",
      "Optimized q @ cavity = (-50+50j)\n",
      "by hand: q @ cavity = (-45.013+51.123j)\n"
     ]
    }
   ],
   "source": [
    "lamda = 0.684e-3  # mm\n",
    "q_laser = 200 + 1j*50\n",
    "q_target = -50 + 50j\n",
    "init_params = (700,1000,50)        # f1,d1,f2,d2\n",
    "bnds = ( (50,1000),(50,1000),(40,60) )\n",
    "\n",
    "res = optimize.minimize( partial(mode_matching_function,q_laser,q_target),\n",
    "                    init_params,\n",
    "                    method='SLSQP',\n",
    "                    bounds=bnds ) \n",
    "\n",
    "print(res)\n",
    "\n",
    "q_opt = propagate( make_mode_matching_system(res.x).T, q_laser )\n",
    "print( \"Target q @ cavity = {}\".format(np.round(q_target,3)) )\n",
    "print( \"Optimized q @ cavity = {}\".format(np.round(q_opt,3)) )\n",
    "\n",
    "# testing tolerances by hand\n",
    "mm = make_mode_matching_system([150,150,55])\n",
    "print( \"by hand: q @ cavity = {}\".format(np.round(propagate(mm.T,q_laser),3) ) )\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}