fix for spherical Bessel functions (sph_jn/yn) compatibility

fix for spherical Bessel functions (sph_jn/yn) compatibility between scipy ver. 0.x.x and 1.x.x

Author: com-py

ch12/Program_12.2_qmscatt.ipynb

@@ -0,0 +1,108 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "#\n",
+    "# Program 12.2: Partial wave method (qmscatt.ipynb)\n",
+    "# J Wang, Computational modeling and visualization with Python\n",
+    "#\n",
+    "\n",
+    "import scipy\n",
+    "if scipy.__version__[0] < '1':  # sph. Bessel func changed after ver 1.0.0\n",
+    "    from scipy.special import sph_jn, sph_yn, lpn\n",
+    "else:\n",
+    "    from scipy.special import spherical_jn as sph_jn, spherical_yn as sph_yn, lpn\n",
+    "import matplotlib.pyplot as plt, numpy as np, ode\n",
+    "%matplotlib notebook\n",
+    "\n",
+    "def V(r):                       # Yukawa potential\n",
+    "    Z, sa = 2., 1.0             # nuclear charge, screening length\n",
+    "    return -Z*np.exp(-r/sa)/r\n",
+    "    \n",
+    "def f(r):                       # Sch eqn in Numerov form\n",
+    "    return 2*(E - V(r)) - L*(L+1)/(r*r)\n",
+    "\n",
+    "def wf(M, xm):                  # find w.f. and deriv at xm\n",
+    "    c = (h*h)/6.\n",
+    "    wfup, nup = ode.numerov(f, [0,.1], M, xL, h)    # 1 step past xm\n",
+    "    dup = ((1+c*f(xm+h))*wfup[-1] - (1+c*f(xm-h))*wfup[-3])/(h+h)\n",
+    "    return wfup, dup/wfup[-2]\n",
+    "\n",
+    "xL, a, M = 0., 10., 200                 # limits, matching point\n",
+    "h, Lmax, E =(a-xL)/M, 15, 2.            # step size, max L, energy\n",
+    "\n",
+    "k, ps = np.sqrt(2*E), np.zeros(Lmax+1)  # wave vector, phase shift\n",
+    "if scipy.__version__[0] < '1':\n",
+    "    jl, dj = sph_jn(Lmax, k*a)          # (j_l, j_l') tuple     \n",
+    "    nl, dn = sph_yn(Lmax, k*a)          # (n_l, n_l')\n",
+    "else:\n",
+    "    Lrange = range(Lmax + 1)\n",
+    "    jl, dj = sph_jn(Lrange, k*a, False), sph_jn(Lrange, k*a, True) # (j_l, j_l')\n",
+    "    nl, dn = sph_yn(Lrange, k*a, False), sph_yn(Lrange, k*a, True) # (n_l, n_l')\n",
+    "\n",
+    "for L in range(Lmax+1):\n",
+    "    u, g = wf(M, a)                     # g= u'/u\n",
+    "    x = np.arctan(((g*a-1)*jl[L] - k*a*dj[L])/    # phase shift \n",
+    "                  ((g*a-1)*nl[L] - k*a*dn[L]))\n",
+    "    while (x < 0.0): x += np.pi         # handle jumps by pi \n",
+    "    ps[L] = x\n",
+    "\n",
+    "theta, sigma = np.linspace(0., np.pi, 100), []\n",
+    "cos, La = np.cos(theta), np.arange(1,2*Lmax+2,2)\n",
+    "for x in cos:                               # calc cross section\n",
+    "    pl = lpn(Lmax, x)[0]                    # Legendre polynomial \n",
+    "    fl = La*np.exp(1j*ps)*np.sin(ps)*pl     # amplitude \n",
+    "    sigma.append(np.abs(np.sum(fl))**2/(k*k))\n",
+    "        \n",
+    "plt.figure()                                # plot phase shift vs L\n",
+    "plt.plot(range(Lmax+1), ps, '-o')\n",
+    "plt.xlabel('$l$'), plt.ylabel(r'$\\delta_l$', fontsize=16)\n",
+    "\n",
+    "plt.figure()        \n",
+    "plt.plot(theta*180/np.pi, sigma)            # plot cross sections   \n",
+    "xtck = [0, 30, 60, 90, 120, 150, 180]\n",
+    "plt.xticks(xtck, [repr(i) for i in xtck])   # custom ticks \n",
+    "plt.xlabel(r'$\\theta$ (deg)')\n",
+    "plt.ylabel(r'$\\sigma(\\theta)$ (a.u.)'), plt.semilogy()\n",
+    "\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}