simple qx working code

.ipynb_checkpoints/Untitled4-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

.ipynb_checkpoints/my_tmm-checkpoint.ipynb

@@ -0,0 +1,274 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "code_folding": []
+   },
+   "outputs": [],
+   "source": [
+    "from __future__ import division, print_function, absolute_import\n",
+    "\n",
+    "from numpy import cos, inf, zeros, array, exp, conj, nan, isnan, pi, sin\n",
+    "\n",
+    "cimport numpy as np\n",
+    "cimport scipy as sp\n",
+    "\n",
+    "import sys\n",
+    "EPSILON = sys.float_info.epsilon # typical floating-point calculation error\n",
+    "\n",
+    "def make_2x2_array(float complex a, float complex b, float complex c, float complex d):\n",
+    "    \"\"\"\n",
+    "    Makes a 2x2 numpy array of [[a,b],[c,d]]\n",
+    "\n",
+    "    Same as \"numpy.array([[a,b],[c,d]], dtype=float)\", but ten times faster\n",
+    "    \"\"\"\n",
+    "    my_array = np.empty((2,2), dtype=complex)\n",
+    "    my_array[0,0] = a\n",
+    "    my_array[0,1] = b\n",
+    "    my_array[1,0] = c\n",
+    "    my_array[1,1] = d\n",
+    "    return my_array\n",
+    "\n",
+    "def is_forward_angle(n, theta):\n",
+    "    \"\"\"\n",
+    "    if a wave is traveling at angle theta from normal in a medium with index n,\n",
+    "    calculate whether or not this is the forward-traveling wave (i.e., the one\n",
+    "    going from front to back of the stack, like the incoming or outgoing waves,\n",
+    "    but unlike the reflected wave). For real n & theta, the criterion is simply\n",
+    "    -pi/2 < theta < pi/2, but for complex n & theta, it's more complicated.\n",
+    "    See https://arxiv.org/abs/1603.02720 appendix D. If theta is the forward\n",
+    "    angle, then (pi-theta) is the backward angle and vice-versa.\n",
+    "    \"\"\"\n",
+    "    assert n.real * n.imag >= 0, (\"For materials with gain, it's ambiguous which \"\n",
+    "                                  \"beam is incoming vs outgoing. See \"\n",
+    "                                  \"https://arxiv.org/abs/1603.02720 Appendix C.\\n\"\n",
+    "                                  \"n: \" + str(n) + \"   angle: \" + str(theta))\n",
+    "    ncostheta = n * np.cos(theta)\n",
+    "    if abs(ncostheta.imag) > 100 * EPSILON:\n",
+    "        # Either evanescent decay or lossy medium. Either way, the one that\n",
+    "        # decays is the forward-moving wave\n",
+    "        answer = (ncostheta.imag > 0)\n",
+    "    else:\n",
+    "        # Forward is the one with positive Poynting vector\n",
+    "        # Poynting vector is Re[n cos(theta)] for s-polarization or\n",
+    "        # Re[n cos(theta*)] for p-polarization, but it turns out they're consistent\n",
+    "        # so I'll just assume s then check both below\n",
+    "        answer = (ncostheta.real > 0)\n",
+    "    # convert from numpy boolean to the normal Python boolean\n",
+    "    answer = bool(answer)\n",
+    "    # double-check the answer ... can't be too careful!\n",
+    "#     error_string = (\"It's not clear which beam is incoming vs outgoing. Weird\"\n",
+    "#                     \" index maybe?\\n\"\n",
+    "#                     \"n: \" + str(n) + \"   angle: \" + str(theta))\n",
+    "#     if answer is True:\n",
+    "#         assert ncostheta.imag > -100 * EPSILON, error_string\n",
+    "#         assert ncostheta.real > -100 * EPSILON, error_string\n",
+    "#         assert (n * cos(theta.conjugate())).real > -100 * EPSILON, error_string\n",
+    "#     else:\n",
+    "#         assert ncostheta.imag < 100 * EPSILON, error_string\n",
+    "#         assert ncostheta.real < 100 * EPSILON, error_string\n",
+    "#         assert (n * cos(theta.conjugate())).real < 100 * EPSILON, error_string\n",
+    "    return answer\n",
+    "\n",
+    "def snell(n_1, n_2, th_1):\n",
+    "    \"\"\"\n",
+    "    return angle theta in layer 2 with refractive index n_2, assuming\n",
+    "    it has angle th_1 in layer with refractive index n_1. Use Snell's law. Note\n",
+    "    that \"angles\" may be complex!!\n",
+    "    \"\"\"\n",
+    "    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They\n",
+    "    # give different results e.g. for arcsin(2).)\n",
+    "    th_2_guess = sp.arcsin(n_1*np.sin(th_1) / n_2)\n",
+    "    if is_forward_angle(n_2, th_2_guess):\n",
+    "        return th_2_guess\n",
+    "    else:\n",
+    "        return pi - th_2_guess\n",
+    "\n",
+    "def list_snell(n_list, th_0):\n",
+    "    \"\"\"\n",
+    "    return list of angle theta in each layer based on angle th_0 in layer 0,\n",
+    "    using Snell's law. n_list is index of refraction of each layer. Note that\n",
+    "    \"angles\" may be complex!!\n",
+    "    \"\"\"\n",
+    "    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They\n",
+    "    # give different results e.g. for arcsin(2).)\n",
+    "    angles = sp.arcsin(n_list[0]*np.sin(th_0) / n_list)\n",
+    "    # The first and last entry need to be the forward angle (the intermediate\n",
+    "    # layers don't matter, see https://arxiv.org/abs/1603.02720 Section 5)\n",
+    "    if not is_forward_angle(n_list[0], angles[0]):\n",
+    "        angles[0] = pi - angles[0]\n",
+    "    if not is_forward_angle(n_list[-1], angles[-1]):\n",
+    "        angles[-1] = pi - angles[-1]\n",
+    "    return angles\n",
+    "\n",
+    "\n",
+    "def interface_r(polarization, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    reflection amplitude (from Fresnel equations)\n",
+    "\n",
+    "    polarization is either \"s\" or \"p\" for polarization\n",
+    "\n",
+    "    n_i, n_f are (complex) refractive index for incident and final\n",
+    "\n",
+    "    th_i, th_f are (complex) propegation angle for incident and final\n",
+    "    (in radians, where 0=normal). \"th\" stands for \"theta\".\n",
+    "    \"\"\"\n",
+    "    if polarization == 's':\n",
+    "        return ((n_i * np.cos(th_i) - n_f * np.cos(th_f)) /\n",
+    "                (n_i * np.cos(th_i) + n_f * np.cos(th_f)))\n",
+    "    else:\n",
+    "        return ((n_f * np.cos(th_i) - n_i * np.cos(th_f)) /\n",
+    "                (n_f * np.cos(th_i) + n_i * np.cos(th_f)))\n",
+    "\n",
+    "def interface_t(polarization, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    transmission amplitude (frem Fresnel equations)\n",
+    "\n",
+    "    polarization is either \"s\" or \"p\" for polarization\n",
+    "\n",
+    "    n_i, n_f are (complex) refractive index for incident and final\n",
+    "\n",
+    "    th_i, th_f are (complex) propegation angle for incident and final\n",
+    "    (in radians, where 0=normal). \"th\" stands for \"theta\".\n",
+    "    \"\"\"\n",
+    "    if polarization == 's':\n",
+    "        return 2 * n_i * np.cos(th_i) / (n_i * np.cos(th_i) + n_f * np.cos(th_f))\n",
+    "    else:\n",
+    "        return 2 * n_i * np.cos(th_i) / (n_f * np.cos(th_i) + n_i * np.cos(th_f))\n",
+    "\n",
+    "\n",
+    "\n",
+    "def coh_tmm(pol, n_list, d_list, float th_0, float lam_vac):\n",
+    "    \"\"\"\n",
+    "    Main \"coherent transfer matrix method\" calc. Given parameters of a stack,\n",
+    "    calculates everything you could ever want to know about how light\n",
+    "    propagates in it. (If performance is an issue, you can delete some of the\n",
+    "    calculations without affecting the rest.)\n",
+    "\n",
+    "    pol is light polarization, \"s\" or \"p\".\n",
+    "\n",
+    "    n_list is the list of refractive indices, in the order that the light would\n",
+    "    pass through them. The 0'th element of the list should be the semi-infinite\n",
+    "    medium from which the light enters, the last element should be the semi-\n",
+    "    infinite medium to which the light exits (if any exits).\n",
+    "\n",
+    "    th_0 is the angle of incidence: 0 for normal, pi/2 for glancing.\n",
+    "    Remember, for a dissipative incoming medium (n_list[0] is not real), th_0\n",
+    "    should be complex so that n0 sin(th0) is real (intensity is constant as\n",
+    "    a function of lateral position).\n",
+    "\n",
+    "    d_list is the list of layer thicknesses (front to back). Should correspond\n",
+    "    one-to-one with elements of n_list. First and last elements should be \"inf\".\n",
+    "\n",
+    "    lam_vac is vacuum wavelength of the light.\n",
+    "\n",
+    "    * R--reflected wave power (as fraction of incident)\n",
+    "    \"\"\"\n",
+    "    # Convert lists to numpy arrays if they're not already.\n",
+    "    n_list = array(n_list)\n",
+    "    d_list = array(d_list, dtype=float)\n",
+    "\n",
+    "    cdef int num_layers\n",
+    "    num_layers = n_list.size\n",
+    "    \n",
+    "    th_list = list_snell(n_list, th_0)\n",
+    "\n",
+    "    # kz is the z-component of (complex) angular wavevector for forward-moving\n",
+    "    # wave. Positive imaginary part means decaying.\n",
+    "    kz_list = 2 * np.pi * n_list * np.cos(th_list) / lam_vac\n",
+    "\n",
+    "    # delta is the total phase accrued by traveling through a given layer.\n",
+    "    # Ignore warning about inf multiplication\n",
+    "    delta = kz_list * d_list\n",
+    "\n",
+    "\n",
+    "    # t_list[i,j] and r_list[i,j] are transmission and reflection amplitudes,\n",
+    "    # respectively, coming from i, going to j. Only need to calculate this when\n",
+    "    # j=i+1. (2D array is overkill but helps avoid confusion.)\n",
+    "    t_list = zeros((num_layers, num_layers), dtype=complex)\n",
+    "    r_list = zeros((num_layers, num_layers), dtype=complex)\n",
+    "    for i in range(num_layers-1):\n",
+    "        t_list[i,i+1] = interface_t(pol, n_list[i], n_list[i+1],\n",
+    "                                    th_list[i], th_list[i+1])\n",
+    "        r_list[i,i+1] = interface_r(pol, n_list[i], n_list[i+1],\n",
+    "                                    th_list[i], th_list[i+1])\n",
+    "        \n",
+    "    # My M is a bit different than Sernelius's, but Mtilde is the same.\n",
+    "    M_list = zeros((num_layers, 2, 2), dtype=complex)\n",
+    "    for i in range(1, num_layers-1):\n",
+    "        M_list[i] = (1/t_list[i,i+1]) * np.dot(\n",
+    "            make_2x2_array(np.exp(-1j*delta[i]), 0, 0, np.exp(1j*delta[i]),dtype=complex),\n",
+    "            make_2x2_array(1, r_list[i,i+1], r_list[i,i+1], 1))\n",
+    "    Mtilde = make_2x2_array(1, 0, 0, 1)\n",
+    "    for i in range(1, num_layers-1):\n",
+    "        Mtilde = np.dot(Mtilde, M_list[i])\n",
+    "    Mtilde = np.dot(make_2x2_array(1, r_list[0,1], r_list[0,1], 1)/t_list[0,1], Mtilde)\n",
+    "\n",
+    "    # Net complex transmission and reflection amplitudes\n",
+    "    r = Mtilde[1,0]/Mtilde[0,0]\n",
+    "#     t = 1/Mtilde[0,0]\n",
+    "\n",
+    "#     # vw_list[n] = [v_n, w_n]. v_0 and w_0 are undefined because the 0th medium\n",
+    "#     # has no left interface.\n",
+    "#     vw_list = zeros((num_layers, 2), dtype=complex)\n",
+    "#     vw = array([[t],[0]])\n",
+    "#     vw_list[-1,:] = np.transpose(vw)\n",
+    "#     for i in range(num_layers-2, 0, -1):\n",
+    "#         vw = np.dot(M_list[i], vw)\n",
+    "#         vw_list[i,:] = np.transpose(vw)\n",
+    "\n",
+    "    return np.abs(r)**2\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.0"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

.ipynb_checkpoints/my_tmm-checkpoint.py

@@ -0,0 +1,396 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2018-11-01T07:40:08.911910Z",
+     "start_time": "2018-11-01T07:40:08.818797Z"
+    },
+    "code_folding": []
+   },
+   "outputs": [],
+   "source": [
+    "#from __future__ import division, print_function, absolute_import\n",
+    "\n",
+    "from numpy import cos, inf, zeros, array, exp, conj, nan, isnan, pi, sin\n",
+    "\n",
+    "import numpy as np\n",
+    "import scipy as sp\n",
+    "\n",
+    "import sys\n",
+    "EPSILON = sys.float_info.epsilon # typical floating-point calculation error\n",
+    "\n",
+    "def make_2x2_array(a, b, c, d, dtype=float):\n",
+    "    \"\"\"\n",
+    "    Makes a 2x2 numpy array of [[a,b],[c,d]]\n",
+    "\n",
+    "    Same as \"numpy.array([[a,b],[c,d]], dtype=float)\", but ten times faster\n",
+    "    \"\"\"\n",
+    "    my_array = np.empty((2,2), dtype=dtype)\n",
+    "    my_array[0,0] = a\n",
+    "    my_array[0,1] = bggvG\n",
+    "    my_array[1,0] = c\n",
+    "    my_array[1,1] = d\n",
+    "    return my_array\n",
+    "\n",
+    "def is_forward_angle(n, theta):\n",
+    "    \"\"\"\n",
+    "    if a wave is traveling at angle theta from normal in a medium with index n,\n",
+    "    calculate whether or not this is the forward-traveling wave (i.e., the one\n",
+    "    going from front to back of the stack, like the incoming or outgoing waves,\n",
+    "    but unlike the reflected wave). For real n & theta, the criterion is simply\n",
+    "    -pi/2 < theta < pi/2, but for complex n & theta, it's more complicated.\n",
+    "    See https://arxiv.org/abs/1603.02720 appendix D. If theta is the forward\n",
+    "    angle, then (pi-theta) is the backward angle and vice-versa.\n",
+    "    \"\"\"\n",
+    "    assert n.real * n.imag >= 0, (\"For materials with gain, it's ambiguous which \"\n",
+    "                                  \"beam is incoming vs outgoing. See \"\n",
+    "                                  \"https://arxiv.org/abs/1603.02720 Appendix C.\\n\"\n",
+    "                                  \"n: \" + str(n) + \"   angle: \" + str(theta))\n",
+    "    ncostheta = n * cos(theta)\n",
+    "    if abs(ncostheta.imag) > 100 * EPSILON:\n",
+    "        # Either evanescent decay or lossy medium. Either way, the one that\n",
+    "        # decays is the forward-moving wave\n",
+    "        answer = (ncostheta.imag > 0)\n",
+    "    else:\n",
+    "        # Forward is the one with positive Poynting vector\n",
+    "        # Poynting vector is Re[n cos(theta)] for s-polarization or\n",
+    "        # Re[n cos(theta*)] for p-polarization, but it turns out they're consistent\n",
+    "        # so I'll just assume s then check both below\n",
+    "        answer = (ncostheta.real > 0)\n",
+    "    # convert from numpy boolean to the normal Python boolean\n",
+    "    answer = bool(answer)\n",
+    "    # double-check the answer ... can't be too careful!\n",
+    "    error_string = (\"It's not clear which beam is incoming vs outgoing. Weird\"\n",
+    "                    \" index maybe?\\n\"\n",
+    "                    \"n: \" + str(n) + \"   angle: \" + str(theta))\n",
+    "    if answer is True:\n",
+    "        assert ncostheta.imag > -100 * EPSILON, error_string\n",
+    "        assert ncostheta.real > -100 * EPSILON, error_string\n",
+    "        assert (n * cos(theta.conjugate())).real > -100 * EPSILON, error_string\n",
+    "    else:\n",
+    "        assert ncostheta.imag < 100 * EPSILON, error_string\n",
+    "        assert ncostheta.real < 100 * EPSILON, error_string\n",
+    "        assert (n * cos(theta.conjugate())).real < 100 * EPSILON, error_string\n",
+    "    return answer\n",
+    "\n",
+    "def snell(n_1, n_2, th_1):\n",
+    "    \"\"\"\n",
+    "    return angle theta in layer 2 with refractive index n_2, assuming\n",
+    "    it has angle th_1 in layer with refractive index n_1. Use Snell's law. Note\n",
+    "    that \"angles\" may be complex!!\n",
+    "    \"\"\"\n",
+    "    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They\n",
+    "    # give different results e.g. for arcsin(2).)\n",
+    "    th_2_guess = sp.arcsin(n_1*np.sin(th_1) / n_2)\n",
+    "    if is_forward_angle(n_2, th_2_guess):\n",
+    "        return th_2_guess\n",
+    "    else:\n",
+    "        return pi - th_2_guess\n",
+    "\n",
+    "def list_snell(n_list, th_0):\n",
+    "    \"\"\"\n",
+    "    return list of angle theta in each layer based on angle th_0 in layer 0,\n",
+    "    using Snell's law. n_list is index of refraction of each layer. Note that\n",
+    "    \"angles\" may be complex!!\n",
+    "    \"\"\"\n",
+    "    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They\n",
+    "    # give different results e.g. for arcsin(2).)\n",
+    "    angles = sp.arcsin(n_list[0]*np.sin(th_0) / n_list)\n",
+    "    # The first and last entry need to be the forward angle (the intermediate\n",
+    "    # layers don't matter, see https://arxiv.org/abs/1603.02720 Section 5)\n",
+    "    if not is_forward_angle(n_list[0], angles[0]):\n",
+    "        angles[0] = pi - angles[0]\n",
+    "    if not is_forward_angle(n_list[-1], angles[-1]):\n",
+    "        angles[-1] = pi - angles[-1]\n",
+    "    return angles\n",
+    "\n",
+    "\n",
+    "def interface_r(polarization, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    reflection amplitude (from Fresnel equations)\n",
+    "\n",
+    "    polarization is either \"s\" or \"p\" for polarization\n",
+    "\n",
+    "    n_i, n_f are (complex) refractive index for incident and final\n",
+    "\n",
+    "    th_i, th_f are (complex) propegation angle for incident and final\n",
+    "    (in radians, where 0=normal). \"th\" stands for \"theta\".\n",
+    "    \"\"\"\n",
+    "    if polarization == 's':\n",
+    "        return ((n_i * cos(th_i) - n_f * cos(th_f)) /\n",
+    "                (n_i * cos(th_i) + n_f * cos(th_f)))\n",
+    "    elif polarization == 'p':\n",
+    "        return ((n_f * cos(th_i) - n_i * cos(th_f)) /\n",
+    "                (n_f * cos(th_i) + n_i * cos(th_f)))\n",
+    "    else:\n",
+    "        raise ValueError(\"Polarization must be 's' or 'p'\")\n",
+    "\n",
+    "def interface_t(polarization, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    transmission amplitude (frem Fresnel equations)\n",
+    "\n",
+    "    polarization is either \"s\" or \"p\" for polarization\n",
+    "\n",
+    "    n_i, n_f are (complex) refractive index for incident and final\n",
+    "\n",
+    "    th_i, th_f are (complex) propegation angle for incident and final\n",
+    "    (in radians, where 0=normal). \"th\" stands for \"theta\".\n",
+    "    \"\"\"\n",
+    "    if polarization == 's':\n",
+    "        return 2 * n_i * cos(th_i) / (n_i * cos(th_i) + n_f * cos(th_f))\n",
+    "    elif polarization == 'p':\n",
+    "        return 2 * n_i * cos(th_i) / (n_f * cos(th_i) + n_i * cos(th_f))\n",
+    "    else:\n",
+    "        raise ValueError(\"Polarization must be 's' or 'p'\")\n",
+    "\n",
+    "def R_from_r(r):\n",
+    "    \"\"\"\n",
+    "    Calculate reflected power R, starting with reflection amplitude r.\n",
+    "    \"\"\"\n",
+    "    return abs(r)**2\n",
+    "\n",
+    "def T_from_t(pol, t, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    Calculate transmitted power T, starting with transmission amplitude t.\n",
+    "\n",
+    "    n_i,n_f are refractive indices of incident and final medium.\n",
+    "\n",
+    "    th_i, th_f are (complex) propegation angles through incident & final medium\n",
+    "    (in radians, where 0=normal). \"th\" stands for \"theta\".\n",
+    "\n",
+    "    In the case that n_i,n_f,th_i,th_f are real, formulas simplify to\n",
+    "    T=|t|^2 * (n_f cos(th_f)) / (n_i cos(th_i)).\n",
+    "\n",
+    "    See manual for discussion of formulas\n",
+    "    \"\"\"\n",
+    "    if pol == 's':\n",
+    "        return abs(t**2) * (((n_f*cos(th_f)).real) / (n_i*cos(th_i)).real)\n",
+    "    elif pol == 'p':\n",
+    "        return abs(t**2) * (((n_f*conj(cos(th_f))).real) /\n",
+    "                                (n_i*conj(cos(th_i))).real)\n",
+    "    else:\n",
+    "        raise ValueError(\"Polarization must be 's' or 'p'\")\n",
+    "\n",
+    "def power_entering_from_r(pol, r, n_i, th_i):\n",
+    "    \"\"\"\n",
+    "    Calculate the power entering the first interface of the stack, starting with\n",
+    "    reflection amplitude r. Normally this equals 1-R, but in the unusual case\n",
+    "    that n_i is not real, it can be a bit different than 1-R. See manual.\n",
+    "\n",
+    "    n_i is refractive index of incident medium.\n",
+    "\n",
+    "    th_i is (complex) propegation angle through incident medium\n",
+    "    (in radians, where 0=normal). \"th\" stands for \"theta\".\n",
+    "    \"\"\"\n",
+    "    if pol == 's':\n",
+    "        return ((n_i*cos(th_i)*(1+conj(r))*(1-r)).real\n",
+    "                     / (n_i*cos(th_i)).real)\n",
+    "    elif pol == 'p':\n",
+    "        return ((n_i*conj(cos(th_i))*(1+r)*(1-conj(r))).real\n",
+    "                      / (n_i*conj(cos(th_i))).real)\n",
+    "    else:\n",
+    "        raise ValueError(\"Polarization must be 's' or 'p'\")\n",
+    "\n",
+    "def interface_R(polarization, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    Fraction of light intensity reflected at an interface.\n",
+    "    \"\"\"\n",
+    "    r = interface_r(polarization, n_i, n_f, th_i, th_f)\n",
+    "    return R_from_r(r)\n",
+    "\n",
+    "def interface_T(polarization, n_i, n_f, th_i, th_f):\n",
+    "    \"\"\"\n",
+    "    Fraction of light intensity transmitted at an interface.\n",
+    "    \"\"\"\n",
+    "    t = interface_t(polarization, n_i, n_f, th_i, th_f)\n",
+    "    return T_from_t(polarization, t, n_i, n_f, th_i, th_f)\n",
+    "\n",
+    "def coh_tmm(pol, n_list, d_list, th_0, lam_vac):\n",
+    "    \"\"\"\n",
+    "    Main \"coherent transfer matrix method\" calc. Given parameters of a stack,\n",
+    "    calculates everything you could ever want to know about how light\n",
+    "    propagates in it. (If performance is an issue, you can delete some of the\n",
+    "    calculations without affecting the rest.)\n",
+    "\n",
+    "    pol is light polarization, \"s\" or \"p\".\n",
+    "\n",
+    "    n_list is the list of refractive indices, in the order that the light would\n",
+    "    pass through them. The 0'th element of the list should be the semi-infinite\n",
+    "    medium from which the light enters, the last element should be the semi-\n",
+    "    infinite medium to which the light exits (if any exits).\n",
+    "\n",
+    "    th_0 is the angle of incidence: 0 for normal, pi/2 for glancing.\n",
+    "    Remember, for a dissipative incoming medium (n_list[0] is not real), th_0\n",
+    "    should be complex so that n0 sin(th0) is real (intensity is constant as\n",
+    "    a function of lateral position).\n",
+    "\n",
+    "    d_list is the list of layer thicknesses (front to back). Should correspond\n",
+    "    one-to-one with elements of n_list. First and last elements should be \"inf\".\n",
+    "\n",
+    "    lam_vac is vacuum wavelength of the light.\n",
+    "\n",
+    "    Outputs the following as a dictionary (see manual for details)\n",
+    "\n",
+    "    * r--reflection amplitude\n",
+    "    * t--transmission amplitude\n",
+    "    * R--reflected wave power (as fraction of incident)\n",
+    "    * T--transmitted wave power (as fraction of incident)\n",
+    "    * power_entering--Power entering the first layer, usually (but not always)\n",
+    "      equal to 1-R (see manual).\n",
+    "    * vw_list-- n'th element is [v_n,w_n], the forward- and backward-traveling\n",
+    "      amplitudes, respectively, in the n'th medium just after interface with\n",
+    "      (n-1)st medium.\n",
+    "    * kz_list--normal component of complex angular wavenumber for\n",
+    "      forward-traveling wave in each layer.\n",
+    "    * th_list--(complex) propagation angle (in radians) in each layer\n",
+    "    * pol, n_list, d_list, th_0, lam_vac--same as input\n",
+    "\n",
+    "    \"\"\"\n",
+    "    # Convert lists to numpy arrays if they're not already.\n",
+    "    n_list = array(n_list)\n",
+    "    d_list = array(d_list, dtype=float)\n",
+    "\n",
+    "    # Input tests\n",
+    "    # if ((hasattr(lam_vac, 'size') and lam_vac.size > 1)\n",
+    "    #       or (hasattr(th_0, 'size') and th_0.size > 1)):\n",
+    "    #     raise ValueError('This function is not vectorized; you need to run one '\n",
+    "    #                      'calculation at a time (1 wavelength, 1 angle, etc.)')\n",
+    "    # if (n_list.ndim != 1) or (d_list.ndim != 1) or (n_list.size != d_list.size):\n",
+    "    #     raise ValueError(\"Problem with n_list or d_list!\")\n",
+    "    # assert d_list[0] == d_list[-1] == inf, 'd_list must start and end with inf!'\n",
+    "    # assert abs((n_list[0]*np.sin(th_0)).imag) < 100*EPSILON, 'Error in n0 or th0!'\n",
+    "    # assert is_forward_angle(n_list[0], th_0), 'Error in n0 or th0!'\n",
+    "    num_layers = n_list.size\n",
+    "\n",
+    "    # th_list is a list with, for each layer, the angle that the light travels\n",
+    "    # through the layer. Computed with Snell's law. Note that the \"angles\" may be\n",
+    "    # complex!\n",
+    "    th_list = list_snell(n_list, th_0)\n",
+    "\n",
+    "    # kz is the z-component of (complex) angular wavevector for forward-moving\n",
+    "    # wave. Positive imaginary part means decaying.\n",
+    "    kz_list = 2 * np.pi * n_list * cos(th_list) / lam_vac\n",
+    "\n",
+    "    # delta is the total phase accrued by traveling through a given layer.\n",
+    "    # Ignore warning about inf multiplication\n",
+    "    olderr = sp.seterr(invalid='ignore')\n",
+    "    delta = kz_list * d_list\n",
+    "    sp.seterr(**olderr)\n",
+    "\n",
+    "    # For a very opaque layer, reset delta to avoid divide-by-0 and similar\n",
+    "    # errors. The criterion imag(delta) > 35 corresponds to single-pass\n",
+    "    # transmission < 1e-30 --- small enough that the exact value doesn't\n",
+    "    # matter.\n",
+    "    # for i in range(1, num_layers-1):\n",
+    "    #     if delta[i].imag > 35:\n",
+    "    #         delta[i] = delta[i].real + 35j\n",
+    "    #         if 'opacity_warning' not in globals():\n",
+    "    #             global opacity_warning\n",
+    "    #             opacity_warning = True\n",
+    "    #             print(\"Warning: Layers that are almost perfectly opaque \"\n",
+    "    #                   \"are modified to be slightly transmissive, \"\n",
+    "    #                   \"allowing 1 photon in 10^30 to pass through. It's \"\n",
+    "    #                   \"for numerical stability. This warning will not \"\n",
+    "    #                   \"be shown again.\")\n",
+    "\n",
+    "    # t_list[i,j] and r_list[i,j] are transmission and reflection amplitudes,\n",
+    "    # respectively, coming from i, going to j. Only need to calculate this when\n",
+    "    # j=i+1. (2D array is overkill but helps avoid confusion.)\n",
+    "    t_list = zeros((num_layers, num_layers), dtype=complex)\n",
+    "    r_list = zeros((num_layers, num_layers), dtype=complex)\n",
+    "    for i in range(num_layers-1):\n",
+    "        t_list[i,i+1] = interface_t(pol, n_list[i], n_list[i+1],\n",
+    "                                    th_list[i], th_list[i+1])\n",
+    "        r_list[i,i+1] = interface_r(pol, n_list[i], n_list[i+1],\n",
+    "                                    th_list[i], th_list[i+1])\n",
+    "    # At the interface between the (n-1)st and nth material, let v_n be the\n",
+    "    # amplitude of the wave on the nth side heading forwards (away from the\n",
+    "    # boundary), and let w_n be the amplitude on the nth side heading backwards\n",
+    "    # (towards the boundary). Then (v_n,w_n) = M_n (v_{n+1},w_{n+1}). M_n is\n",
+    "    # M_list[n]. M_0 and M_{num_layers-1} are not defined.\n",
+    "    # My M is a bit different than Sernelius's, but Mtilde is the same.\n",
+    "    M_list = zeros((num_layers, 2, 2), dtype=complex)\n",
+    "    for i in range(1, num_layers-1):\n",
+    "        M_list[i] = (1/t_list[i,i+1]) * np.dot(\n",
+    "            make_2x2_array(exp(-1j*delta[i]), 0, 0, exp(1j*delta[i]),\n",
+    "                           dtype=complex),\n",
+    "            make_2x2_array(1, r_list[i,i+1], r_list[i,i+1], 1, dtype=complex))\n",
+    "    Mtilde = make_2x2_array(1, 0, 0, 1, dtype=complex)\n",
+    "    for i in range(1, num_layers-1):\n",
+    "        Mtilde = np.dot(Mtilde, M_list[i])\n",
+    "    Mtilde = np.dot(make_2x2_array(1, r_list[0,1], r_list[0,1], 1,\n",
+    "                                   dtype=complex)/t_list[0,1], Mtilde)\n",
+    "\n",
+    "    # Net complex transmission and reflection amplitudes\n",
+    "    r = Mtilde[1,0]/Mtilde[0,0]\n",
+    "    t = 1/Mtilde[0,0]\n",
+    "\n",
+    "    # vw_list[n] = [v_n, w_n]. v_0 and w_0 are undefined because the 0th medium\n",
+    "    # has no left interface.\n",
+    "    vw_list = zeros((num_layers, 2), dtype=complex)\n",
+    "    vw = array([[t],[0]])\n",
+    "    vw_list[-1,:] = np.transpose(vw)\n",
+    "    for i in range(num_layers-2, 0, -1):\n",
+    "        vw = np.dot(M_list[i], vw)\n",
+    "        vw_list[i,:] = np.transpose(vw)\n",
+    "\n",
+    "    # Net transmitted and reflected power, as a proportion of the incoming light\n",
+    "    # power.\n",
+    "    #R = R_from_r(r)\n",
+    "\n",
+    "    return np.abs(r)**2\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.0"
+  },
+  "latex_envs": {
+   "LaTeX_envs_menu_present": true,
+   "autoclose": false,
+   "autocomplete": true,
+   "bibliofile": "biblio.bib",
+   "cite_by": "apalike",
+   "current_citInitial": 1,
+   "eqLabelWithNumbers": true,
+   "eqNumInitial": 1,
+   "hotkeys": {
+    "equation": "Ctrl-E",
+    "itemize": "Ctrl-I"
+   },
+   "labels_anchors": false,
+   "latex_user_defs": false,
+   "report_style_numbering": false,
+   "user_envs_cfg": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

de2.py

@@ -1,7 +1,8 @@
 import numpy as np
+import random
 
 
-def de_c(
+def de_cp(
     fobj,  # objective function, serial eval
     bounds,  # bounds array for the population
     pop=[],  # already initialized population
@@ -14,6 +15,61 @@ def de_c(
     **kwargs
 ):
     """
+        This function performs diff evolution using fobj evaluated in parallel  
+        returns the  last pop, fitness of pop, best individual and opt history 
+    """
+    dimensions = len(bounds)
+    history=[]
+    min_b, max_b = np.asarray(bounds).T
+    if pop == []: 
+        pop = np.random.rand(popsize, dimensions)
+    diff = np.fabs(min_b - max_b)
+    pop_denorm = min_b + pop * diff
+    fitness = np.asarray(fobj(pop_denorm, **kwargs))
+    best_idx = np.argmin(fitness)
+    best = pop_denorm[best_idx]
+    for i in range(its):
+        trialarr = np.zeros((popsize, dimensions))
+        for j in range(popsize):
+            idxs = [idx for idx in range(popsize) if idx != j]
+            a, b, c = pop[np.random.choice(idxs, 3, replace = False)]
+            #mutant = np.clip(a + mut * (b - c), 0, 1)
+            mutant = np.clip(pop[best_idx] + mut * (b - c), 0, 1)
+            cross_points = np.random.rand(dimensions) < crossp
+            if not np.any(cross_points):
+                cross_points[np.random.randint(0, dimensions)] = True
+            trialarr[j] = np.where(cross_points, mutant, pop[j])
+        trial_denorm = min_b + trialarr * diff
+        f = fobj(trial_denorm, **kwargs)
+        for j in range(popsize):
+            if f[j] < fitness[j]:
+                fitness[j] = f[j]
+                pop[j] = trialarr[j]
+                if f[j] < fitness[best_idx]:
+                    best_idx = j
+                    best = trial_denorm[j]
+        if i%5 == 0:
+            print(i, fitness[best_idx])
+        history.append([i, fitness[best_idx]])
+    return best, fitness[best_idx], history
+
+    
+    
+    
+
+def de_c(
+    fobj,  # objective function, serial eval
+    bounds,  # bounds array for the population
+    pop=[],  # already initialized population
+    history=[],  # previous history that will be appended to
+    it_start=0,  # the start index of the iteration
+    mut=0.8,  # mutation rate
+    crossp=0.8,  # crossover
+    popsize=20,  # the population size
+    its=1000,  # the number of iterations needed to run for
+    **kwargs
+):
+    """
         This function performs diff evolution using fobj evaluated serially.
         returns the  best, fitness of the best, and opt history
     """
@@ -30,8 +86,13 @@ def de_c(
         for j in range(popsize):
             idxs = [idx for idx in range(popsize) if idx != j]
             a, b, c = pop[np.random.choice(idxs, 3, replace=False)]
-            #TODO: Here we can chage the strategy for instance
-            mutant = np.clip(pop[best_idx] + mut * (b - c), 0, 1)
+            mutant = np.clip(a + mut * (b - c), 0, 1)
+#             #TODO: Here we can chage the strategy for instance
+#             if i >=150:
+#                 if random.random() >=0.82: 
+#                     mutant = np.clip(a + mut * (b - c), 0, 1)
+#             else:
+#                 mutant = np.clip(pop[best_idx] + mut * (b - c), 0, 1)
             cross_points = np.random.rand(dimensions) < crossp
             if not np.any(cross_points):
                 cross_points[np.random.randint(0, dimensions)] = True
@@ -44,8 +105,8 @@ def de_c(
                 if f < fitness[best_idx]:
                     best_idx = j
                     best = trial_denorm
-#         if i % 20 == 0:
-#             print(i, fitness[best_idx])
+        if i % 5 == 0:
+            print(i, fitness[best_idx])
         history.append([i + it_start, fitness[best_idx]])
     return best, fitness[best_idx], history
 

makeqx.py

@@ -6,6 +6,11 @@ import random
 from tmm import coh_tmm
 import multiprocessing as mp
 
+#import pyximport; pyximport.install(pyimport=True)
+#import mtmm as mtmm
+
+
+
 #make a materials dictionary
 matsdict = {
   1: './materials/gold.dat',
@@ -86,7 +91,7 @@ def make_qx(num_layers):
     rwalk = bet01(rwalk)
 
     # some random number seeds
-    samps = int(np.random.uniform(num_layers/10, 9*num_layers/10, 1))
+    samps = int(np.random.uniform(int(num_layers/10)+2, int(9*num_layers/10), 1))
     scales = np.random.uniform(0.05, 0.75, 1)
 
     s = np.linspace(-1, num_layers + 1, samps)
@@ -97,6 +102,9 @@ def make_qx(num_layers):
         pers = np.random.uniform(0.5, 2.2, 1)
         fr_p[ctr] = np.clip((0.5 + pnoise1(x, octaves=octa, persistence=pers)), 0, 1)
 
+    
+    
+
     if np.random.randint(2):
         f = interpolate.interp1d(s, fr, kind="quadratic", assume_sorted=True)
     else:
@@ -108,18 +116,21 @@ def make_qx(num_layers):
     return rwalk 
 
 
-def digitize_qx(q_x):
+def digitize_qx(q_x, dlevels=2):
     """TODO: Docstring for digitize_qx.
     :q_x: TODO
     :returns: TODO
     """
-    bins = np.array([-0.02, 0.2, 0.4, 0.6, 0.8, 1.01])
+    #bins = np.array([-0.02, 0.2, 0.4, 0.6, 0.8, 1.01])
+    bins = np.linspace(0,1, num=dlevels+1, endpoint=True)
+    bins[0] = -0.02
+    bins[-1] = 1.01
     dig = (np.digitize(q_x, bins, right=True)) - 1
-    act_r = np.linspace(0, 1, 5, endpoint=True)
+    act_r = np.linspace(0, 1, num=dlevels, endpoint=True)
     return act_r[dig]
 
 
-def make_nxdx(qx, wavelen, n_substrate=1.52, layer_thickness=5.0):
+def make_nxdx(qx, cthick, wavelen, n_substrate=1.52):
     """TODO: Docstring for make_nxdx.
     :qx: TODO
     :n_substrate: TODO
@@ -127,17 +138,56 @@ def make_nxdx(qx, wavelen, n_substrate=1.52, layer_thickness=5.0):
     :returns: TODO
     """
     num_layers = qx.size
-    d_x = num_layers*[layer_thickness]
+    d_x = num_layers*[cthick/num_layers]
     d_x = [np.inf] + d_x + [np.inf]
     sde = LL_mixing(fH = qx, 
-                    n_H = get_nk(datafile=matsdict[3], wavelengths=wavelen), 
-                    n_L = get_nk(datafile=matsdict[4], wavelengths=wavelen))
+                    n_H = 2.58, #get_nk(datafile=matsdict[3], wavelengths=wavelen), 
+                    n_L = 1.45) #get_nk(datafile=matsdict[4], wavelengths=wavelen
+               #   ))
     n_x = [1.0] + sde.tolist() + [n_substrate]
     return d_x, n_x
 
+def tmm_eval2(
+        qx, 
+        cthick,
+        lam_low=400, #in nm 
+        lam_high=1200,
+        lam_pts=100,
+        ang_low=0,   #degrees
+        ang_high=90, 
+        ang_pts=25,
+        n_subs=1.52,
+        ):
+    """TODO: Docstring for tmm_eval.
+    :qx: TODO
+    :lam_low: TODO
+    :#in nm 
+        lam_high: TODO
+    :lam_pts: TODO
+    :ang_low: TODO
+    :#degrees
+        ang_high: TODO
+    :ang_pts: TODO
+    :returns: TODO
+
+    """
+    degree = np.pi/180
+    lams = np.linspace(lam_low, lam_high, endpoint=True, num=lam_pts)
+    thetas = np.linspace(ang_high, ang_low, endpoint=True, num=ang_pts)
+    Rs = np.zeros((thetas.size, lams.size))
+    Rp = np.zeros((thetas.size, lams.size))
+
+    for tidx, theta in enumerate(thetas):
+        for lidx, lam in enumerate(lams):
+            d_x, n_x = make_nxdx(qx=qx, cthick=cthick, wavelen=lam, n_substrate=n_subs)
+            Rs[tidx, lidx] = 100*coh_tmm('s',n_x,d_x, th_0=theta*degree,lam_vac=lam)
+            Rp[tidx, lidx] = 100*coh_tmm('p',n_x,d_x, th_0=theta*degree,lam_vac=lam)
+    return Rs, Rp 
+
 
 def tmm_eval(
         qx, 
+        cthick,
         lam_low=400, #in nm 
         lam_high=1200,
         lam_pts=100,
@@ -145,7 +195,6 @@ def tmm_eval(
         ang_high=90, 
         ang_pts=25,
         n_subs=1.52,
-        l_thick=5
         ):
     """TODO: Docstring for tmm_eval.
     :qx: TODO
@@ -168,38 +217,36 @@ def tmm_eval(
 
     for tidx, theta in enumerate(thetas):
         for lidx, lam in enumerate(lams):
-            d_x, n_x = make_nxdx(qx=qx, wavelen=lam, n_substrate=n_subs, 
-                    layer_thickness=l_thick)
-            Rs[tidx, lidx] = coh_tmm('s',n_x,d_x, th_0=theta*degree,lam_vac=lam)
-            Rp[tidx, lidx] = coh_tmm('p',n_x,d_x, th_0=theta*degree,lam_vac=lam)
+            d_x, n_x = make_nxdx(qx=qx, cthick=cthick, wavelen=lam, n_substrate=n_subs)
+            Rs[tidx, lidx] = 100*coh_tmm('s',n_x,d_x, th_0=theta*degree,lam_vac=lam)
+            Rp[tidx, lidx] = 100*coh_tmm('p',n_x,d_x, th_0=theta*degree,lam_vac=lam)
     return Rs, Rp 
 
 def tmm_eval_wbk(
         qx, 
+        cthick,
         inc_ang,
         lam,
-        n_subs=1.52,
-        l_thick=5
+        n_subs=1.52
         ):
     """TODO: Docstring for tmm_eval.
     :qx: TODO
     :returns: TODO
     """
     degree = np.pi/180
-    d_x, n_x = make_nxdx(qx=qx, wavelen=lam, n_substrate=n_subs, 
-            layer_thickness=l_thick)
-    Rs = coh_tmm('s',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
-    Rp = coh_tmm('p',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
-    return Rs, Rp  
+    d_x, n_x = make_nxdx(qx=qx, cthick=cthick, wavelen=lam, n_substrate=n_subs)
+    Rs = 100*mtmm.coh_tmm('s',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
+#     Rp = 100*coh_tmm('p',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
+    return Rs
     
 def tmm_eval_wsweep(
         qx, 
+        cthick,
         inc_ang,
         lam_low=400, #in nm 
         lam_high=1200,
         lam_pts=100,
         n_subs=1.52,
-        l_thick=5
         ):
     """TODO: Docstring for tmm_eval.
     :qx: TODO
@@ -211,14 +258,13 @@ def tmm_eval_wsweep(
     degree = np.pi/180
     lams = np.linspace(lam_low, lam_high, endpoint=True, num=lam_pts)
     Rs = np.zeros(lams.size)
-    Rp = np.zeros(lams.size)
+    #Rp = np.zeros(lams.size)
 
     for lidx, lam in enumerate(lams):
-            d_x, n_x = make_nxdx(qx=qx, wavelen=lam, n_substrate=n_subs, 
-                    layer_thickness=l_thick)
-            Rs[lidx] = coh_tmm('s',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
-            Rp[lidx] = coh_tmm('p',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
-    return Rs, Rp  
+            d_x, n_x = make_nxdx(qx=qx, cthick=cthick, wavelen=lam, n_substrate=n_subs)
+            Rs[lidx] = 100*coh_tmm('s',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
+            #Rp[lidx] = 100*coh_tmm('p',n_x,d_x, th_0=inc_ang*degree,lam_vac=lam)
+    return Rs  
     
 def tmm_wrapper2(arg):
     args, kwargs = arg
@@ -226,6 +272,7 @@ def tmm_wrapper2(arg):
 
 def tmm_lam_parallel(
     q_x, 
+    cthick,
     inc_ang, 
     n_par=12,
     lam_low=400, 
@@ -237,9 +284,10 @@ def tmm_lam_parallel(
     pool=mp.Pool(n_par)
     lams = np.linspace(lam_low, lam_high, endpoint=True, num=lam_pts)
     for lam in lams:
-        jobs.append((q_x, inc_ang, lam))
+        jobs.append((q_x, cthick, inc_ang, lam))
     arg = [(j, kwargs) for j in jobs]
     answ = np.array(pool.map(tmm_wrapper2, arg))
+    pool.close()
     return answ[:,0], answ[:,1] 
     
     
@@ -248,14 +296,15 @@ def tmm_wrapper(arg):
     args, kwargs = arg
     return tmm_eval_wsweep(*args, **kwargs)
 
-def tmm_evol_parallel(q_x, n_ang= 25, n_par=10, **kwargs):
+def tmm_eval_parallel(q_x, cthick, n_ang= 25, n_par=10, **kwargs):
     jobs = []
     pool=mp.Pool(n_par)
     angs = np.linspace(90, 0, endpoint=True, num=n_ang)
     for ang in angs:
-        jobs.append((q_x, ang))
+        jobs.append((q_x, cthick, ang))
     arg = [(j, kwargs) for j in jobs]
     answ = np.array(pool.map(tmm_wrapper, arg))
+    pool.close()
     return answ[:,0,:], answ[:,1,:] 
 
 

mtmm.pyx

@@ -0,0 +1,215 @@
+
+from __future__ import division, print_function, absolute_import
+
+from numpy import cos, inf, zeros, array, exp, conj, nan, isnan, pi, sin
+
+import numpy as np
+import scipy as sp
+
+import sys
+EPSILON = sys.float_info.epsilon # typical floating-point calculation error
+
+
+def make_2x2_array(a, b, c, d, dtype=float):
+    """
+    Makes a 2x2 numpy array of [[a,b],[c,d]]
+    Same as "numpy.array([[a,b],[c,d]], dtype=float)", but ten times faster
+    """
+    my_array = np.empty((2,2), dtype=dtype)
+    my_array[0,0] = a
+    my_array[0,1] = b
+    my_array[1,0] = c
+    my_array[1,1] = d
+    return my_array
+
+
+
+def is_forward_angle(n, theta):
+    """
+    if a wave is traveling at angle theta from normal in a medium with index n,
+    calculate whether or not this is the forward-traveling wave (i.e., the one
+    going from front to back of the stack, like the incoming or outgoing waves,
+    but unlike the reflected wave). For real n & theta, the criterion is simply
+    -pi/2 < theta < pi/2, but for complex n & theta, it's more complicated.
+    See https://arxiv.org/abs/1603.02720 appendix D. If theta is the forward
+    angle, then (pi-theta) is the backward angle and vice-versa.
+    """
+    assert n.real * n.imag >= 0, ("For materials with gain, it's ambiguous which "
+                                  "beam is incoming vs outgoing. See "
+                                  "https://arxiv.org/abs/1603.02720 Appendix C.\n"
+                                  "n: " + str(n) + "   angle: " + str(theta))
+    ncostheta = n * np.cos(theta)
+    if abs(ncostheta.imag) > 100 * EPSILON:
+        # Either evanescent decay or lossy medium. Either way, the one that
+        # decays is the forward-moving wave
+        answer = (ncostheta.imag > 0)
+    else:
+        # Forward is the one with positive Poynting vector
+        # Poynting vector is Re[n cos(theta)] for s-polarization or
+        # Re[n cos(theta*)] for p-polarization, but it turns out they're consistent
+        # so I'll just assume s then check both below
+        answer = (ncostheta.real > 0)
+    # convert from numpy boolean to the normal Python boolean
+    answer = bool(answer)
+    # double-check the answer ... can't be too careful!
+#     error_string = ("It's not clear which beam is incoming vs outgoing. Weird"
+#                     " index maybe?\n"
+#                     "n: " + str(n) + "   angle: " + str(theta))
+#     if answer is True:
+#         assert ncostheta.imag > -100 * EPSILON, error_string
+#         assert ncostheta.real > -100 * EPSILON, error_string
+#         assert (n * cos(theta.conjugate())).real > -100 * EPSILON, error_string
+#     else:
+#         assert ncostheta.imag < 100 * EPSILON, error_string
+#         assert ncostheta.real < 100 * EPSILON, error_string
+#         assert (n * cos(theta.conjugate())).real < 100 * EPSILON, error_string
+    return answer
+
+def snell(n_1, n_2, th_1):
+    """
+    return angle theta in layer 2 with refractive index n_2, assuming
+    it has angle th_1 in layer with refractive index n_1. Use Snell's law. Note
+    that "angles" may be complex!!
+    """
+    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They
+    # give different results e.g. for arcsin(2).)
+    th_2_guess = sp.arcsin(n_1*np.sin(th_1) / n_2)
+    if is_forward_angle(n_2, th_2_guess):
+        return th_2_guess
+    else:
+        return pi - th_2_guess
+
+def list_snell(n_list, th_0):
+    """
+    return list of angle theta in each layer based on angle th_0 in layer 0,
+    using Snell's law. n_list is index of refraction of each layer. Note that
+    "angles" may be complex!!
+    """
+    # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They
+    # give different results e.g. for arcsin(2).)
+    angles = sp.arcsin(n_list[0]*np.sin(th_0) / n_list)
+    # The first and last entry need to be the forward angle (the intermediate
+    # layers don't matter, see https://arxiv.org/abs/1603.02720 Section 5)
+    if not is_forward_angle(n_list[0], angles[0]):
+        angles[0] = pi - angles[0]
+    if not is_forward_angle(n_list[-1], angles[-1]):
+        angles[-1] = pi - angles[-1]
+    return angles
+
+
+def interface_r(polarization, n_i, n_f, th_i, th_f):
+    """
+    reflection amplitude (from Fresnel equations)
+
+    polarization is either "s" or "p" for polarization
+
+    n_i, n_f are (complex) refractive index for incident and final
+
+    th_i, th_f are (complex) propegation angle for incident and final
+    (in radians, where 0=normal). "th" stands for "theta".
+    """
+    if polarization == 's':
+        return ((n_i * np.cos(th_i) - n_f * np.cos(th_f)) /
+                (n_i * np.cos(th_i) + n_f * np.cos(th_f)))
+    else:
+        return ((n_f * np.cos(th_i) - n_i * np.cos(th_f)) /
+                (n_f * np.cos(th_i) + n_i * np.cos(th_f)))
+
+def interface_t(polarization, n_i, n_f, th_i, th_f):
+    """
+    transmission amplitude (frem Fresnel equations)
+
+    polarization is either "s" or "p" for polarization
+
+    n_i, n_f are (complex) refractive index for incident and final
+
+    th_i, th_f are (complex) propegation angle for incident and final
+    (in radians, where 0=normal). "th" stands for "theta".
+    """
+    if polarization == 's':
+        return 2 * n_i * np.cos(th_i) / (n_i * np.cos(th_i) + n_f * np.cos(th_f))
+    else:
+        return 2 * n_i * np.cos(th_i) / (n_f * np.cos(th_i) + n_i * np.cos(th_f))
+
+
+
+def coh_tmm(pol, n_list, d_list, float th_0, float lam_vac):
+    """
+    Main "coherent transfer matrix method" calc. Given parameters of a stack,
+    calculates everything you could ever want to know about how light
+    propagates in it. (If performance is an issue, you can delete some of the
+    calculations without affecting the rest.)
+
+    pol is light polarization, "s" or "p".
+
+    n_list is the list of refractive indices, in the order that the light would
+    pass through them. The 0'th element of the list should be the semi-infinite
+    medium from which the light enters, the last element should be the semi-
+    infinite medium to which the light exits (if any exits).
+
+    th_0 is the angle of incidence: 0 for normal, pi/2 for glancing.
+    Remember, for a dissipative incoming medium (n_list[0] is not real), th_0
+    should be complex so that n0 sin(th0) is real (intensity is constant as
+    a function of lateral position).
+
+    d_list is the list of layer thicknesses (front to back). Should correspond
+    one-to-one with elements of n_list. First and last elements should be "inf".
+
+    lam_vac is vacuum wavelength of the light.
+
+    * R--reflected wave power (as fraction of incident)
+    """
+    # Convert lists to numpy arrays if they're not already.
+    n_list = array(n_list)
+    d_list = array(d_list, dtype=float)
+
+    cdef int num_layers
+    num_layers = n_list.size
+    
+    th_list = list_snell(n_list, th_0)
+
+    # kz is the z-component of (complex) angular wavevector for forward-moving
+    # wave. Positive imaginary part means decaying.
+    kz_list = 2 * np.pi * n_list * np.cos(th_list) / lam_vac
+
+    # delta is the total phase accrued by traveling through a given layer.
+    # Ignore warning about inf multiplication
+    delta = kz_list * d_list
+
+
+    # t_list[i,j] and r_list[i,j] are transmission and reflection amplitudes,
+    # respectively, coming from i, going to j. Only need to calculate this when
+    # j=i+1. (2D array is overkill but helps avoid confusion.)
+    t_list = zeros((num_layers, num_layers), dtype=complex)
+    r_list = zeros((num_layers, num_layers), dtype=complex)
+    for i in range(num_layers-1):
+        t_list[i,i+1] = interface_t(pol, n_list[i], n_list[i+1],
+                                    th_list[i], th_list[i+1])
+        r_list[i,i+1] = interface_r(pol, n_list[i], n_list[i+1],
+                                    th_list[i], th_list[i+1])
+        
+    # My M is a bit different than Sernelius's, but Mtilde is the same.
+    M_list = zeros((num_layers, 2, 2), dtype=complex)
+    for i in range(1, num_layers-1):
+        M_list[i] = (1/t_list[i,i+1]) * np.dot(
+            make_2x2_array(np.exp(-1j*delta[i]), 0, 0, np.exp(1j*delta[i]), dtype=complex),
+            make_2x2_array(1, r_list[i,i+1], r_list[i,i+1], 1, dtype=complex))
+    Mtilde = make_2x2_array(1, 0, 0, 1)
+    for i in range(1, num_layers-1):
+        Mtilde = np.dot(Mtilde, M_list[i])
+    Mtilde = np.dot(make_2x2_array(1, r_list[0,1], r_list[0,1], 1, dtype=complex)/t_list[0,1], Mtilde)
+
+    # Net complex transmission and reflection amplitudes
+    r = Mtilde[1,0]/Mtilde[0,0]
+#     t = 1/Mtilde[0,0]
+
+#     # vw_list[n] = [v_n, w_n]. v_0 and w_0 are undefined because the 0th medium
+#     # has no left interface.
+#     vw_list = zeros((num_layers, 2), dtype=complex)
+#     vw = array([[t],[0]])
+#     vw_list[-1,:] = np.transpose(vw)
+#     for i in range(num_layers-2, 0, -1):
+#         vw = np.dot(M_list[i], vw)
+#         vw_list[i,:] = np.transpose(vw)
+
+    return np.abs(r)**2

qxplots.py

@@ -9,7 +9,10 @@ def plot_qx(q_x, with_dig=True, with_swatch=True):
     uni = np.linspace(0, q_x.size, q_x.size, endpoint=False)
     q_xd = mkq.digitize_qx(q_x)
     ax = plt.subplot(211)
-    gratim = np.tile(q_x, (int(q_x.size/40),1))
+    if q_x.size >=100:
+        gratim = np.tile(q_x, (int(q_x.size/40),1))
+    else:
+        gratim = np.tile(q_x, (int(q_x.size/4),1))
     #ax.axis('tight')
     ax.axis('off')
     ax.set_axis_off()

rectangletry/array-python-cpp

@@ -0,0 +1 @@
+Subproject commit 0923d8a97a9c7d2909cc7efe60f2e05a28f8f014

setup.py

@@ -0,0 +1,7 @@
+from distutils.core import setup
+from Cython.Build import cythonize
+
+
+setup(
+        ext_modules = cythonize("mtmm.pyx")
+)