Skip to content

Commit e2f0f9d

Browse files
com-pycom-py
authored
Jupyter notebook version added
1 parent 750f583 commit e2f0f9d

File tree

2 files changed

+190
-0
lines changed

2 files changed

+190
-0
lines changed

‎ch12/Program_12.1_deflect.ipynb

+93
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"#\n",
10+
"# Program 12.1: Deflection function (deflect.ipynb)\n",
11+
"# J Wang, Computational modeling and visualization with Python\n",
12+
"#\n",
13+
"\n",
14+
"import matplotlib.pyplot as plt\n",
15+
"import numpy as np, rootfinder as rtf, integral as itg\n",
16+
"%matplotlib notebook\n",
17+
"\n",
18+
"def V(r): # plum potential\n",
19+
" return Z*(3-r*r/(R*R))/(R+R) if (r<R) else Z/r\n",
20+
"\n",
21+
"def fu(u): # f(u), turning point eqn\n",
22+
" return 1 - V(1./u)/E - b*b*u*u\n",
23+
" \n",
24+
"def fx(x): # integrand, called by gauss\n",
25+
" u = umin - x*x\n",
26+
" return 2*x*b/np.sqrt(fu(u))\n",
27+
"\n",
28+
"def xection(theta, ba, da): # cross section\n",
29+
" cs = 0.0 # ba=impact para, da=deflection angle\n",
30+
" for i in range(len(ba)-1):\n",
31+
" if ((theta-da[i])*(theta-da[i+1]) < 0.): # theta bracketed \n",
32+
" db = ba[i+1] - ba[i]\n",
33+
" cs += (ba[i] + db/2.)*abs(db/(da[i+1]-da[i]))\n",
34+
" return cs/np.sin(theta)\n",
35+
" \n",
36+
"Z, R = 1.0, 1.0 # nuclear charge Z, radius of plum potential\n",
37+
"E, b, bmax = 1.2, 0.01, 20. # energy, initial b, bmax\n",
38+
"eps, tiny =1.E-14, 1.E-5 # rel error, u-limit\n",
39+
"ba, theta = [], [] # impact para, deflection\n",
40+
"while (b <= bmax):\n",
41+
" umin = rtf.bisect(fu, tiny, 1./b, eps) # find turning pt\n",
42+
" alpha = itg.gauss(fx, 0., np.sqrt(umin))\n",
43+
" ba.append(b), theta.append(np.pi - 2*alpha)\n",
44+
" b *= 1.02\n",
45+
" \n",
46+
"plt.figure(), plt.plot(ba, theta) # plot deflection function\n",
47+
"plt.xlabel('$b$ (a.u.)'), plt.ylabel('$\\Theta$', rotation=0)\n",
48+
"plt.yticks([0, np.pi/2, np.pi], ['$0$','$\\pi/2$', '$\\pi$'])\n",
49+
"plt.xlim(0,3)\n",
50+
"\n",
51+
"cs, sa = [], np.linspace(0.01, np.pi-.01, 500) # scattering angle\n",
52+
"for x in sa: # calc, plot cross section\n",
53+
" cs.append(xection(x, ba, theta))\n",
54+
"plt.figure(), plt.plot(sa, cs)\n",
55+
"plt.xlabel(r'$\\theta$'), plt.ylabel(r'$\\sigma(\\theta)$')\n",
56+
"plt.xticks([0, np.pi/2, np.pi], ['$0$','$\\pi/2$', '$\\pi$'])\n",
57+
"plt.xlim(0, np.pi), plt.semilogy()\n",
58+
"\n",
59+
"plt.show()"
60+
]
61+
},
62+
{
63+
"cell_type": "code",
64+
"execution_count": null,
65+
"metadata": {
66+
"collapsed": true
67+
},
68+
"outputs": [],
69+
"source": []
70+
}
71+
],
72+
"metadata": {
73+
"kernelspec": {
74+
"display_name": "Python 3",
75+
"language": "python",
76+
"name": "python3"
77+
},
78+
"language_info": {
79+
"codemirror_mode": {
80+
"name": "ipython",
81+
"version": 3
82+
},
83+
"file_extension": ".py",
84+
"mimetype": "text/x-python",
85+
"name": "python",
86+
"nbconvert_exporter": "python",
87+
"pygments_lexer": "ipython3",
88+
"version": "3.6.1"
89+
}
90+
},
91+
"nbformat": 4,
92+
"nbformat_minor": 2
93+
}

‎ch12/Program_12.2_qmscatt.ipynb

+97
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"#\n",
10+
"# Program 12.2: Partial wave method (qmscatt.ipynb)\n",
11+
"# J Wang, Computational modeling and visualization with Python\n",
12+
"#\n",
13+
"\n",
14+
"from scipy.special import sph_jn, sph_yn, lpn\n",
15+
"import matplotlib.pyplot as plt, numpy as np, ode\n",
16+
"%matplotlib notebook\n",
17+
"\n",
18+
"def V(r): # Yukawa potential\n",
19+
" Z, sa = 2., 1.0 # nuclear charge, screening length\n",
20+
" return -Z*np.exp(-r/sa)/r\n",
21+
" \n",
22+
"def f(r): # Sch eqn in Numerov form\n",
23+
" return 2*(E - V(r)) - L*(L+1)/(r*r)\n",
24+
"\n",
25+
"def wf(M, xm): # find w.f. and deriv at xm\n",
26+
" c = (h*h)/6.\n",
27+
" wfup, nup = ode.numerov(f, [0,.1], M, xL, h) # 1 step past xm\n",
28+
" dup = ((1+c*f(xm+h))*wfup[-1] - (1+c*f(xm-h))*wfup[-3])/(h+h)\n",
29+
" return wfup, dup/wfup[-2]\n",
30+
"\n",
31+
"xL, a, M = 0., 10., 200 # limits, matching point\n",
32+
"h, Lmax, E =(a-xL)/M, 15, 2. # step size, max L, energy\n",
33+
"\n",
34+
"k, ps = np.sqrt(2*E), np.zeros(Lmax+1) # wave vector, phase shift\n",
35+
"jl, dj = sph_jn(Lmax, k*a) # (j_l, j_l') tuple \n",
36+
"nl, dn = sph_yn(Lmax, k*a) # (n_l, n_l') \n",
37+
"\n",
38+
"for L in range(Lmax+1):\n",
39+
" u, g = wf(M, a) # g= u'/u\n",
40+
" x = np.arctan(((g*a-1)*jl[L] - k*a*dj[L])/ # phase shift \n",
41+
" ((g*a-1)*nl[L] - k*a*dn[L]))\n",
42+
" while (x < 0.0): x += np.pi # handle jumps by pi \n",
43+
" ps[L] = x\n",
44+
"\n",
45+
"theta, sigma = np.linspace(0., np.pi, 100), []\n",
46+
"cos, La = np.cos(theta), np.arange(1,2*Lmax+2,2)\n",
47+
"for x in cos: # calc cross section\n",
48+
" pl = lpn(Lmax, x)[0] # Legendre polynomial \n",
49+
" fl = La*np.exp(1j*ps)*np.sin(ps)*pl # amplitude \n",
50+
" sigma.append(np.abs(np.sum(fl))**2/(k*k))\n",
51+
" \n",
52+
"plt.figure() # plot phase shift vs L\n",
53+
"plt.plot(range(Lmax+1), ps, '-o')\n",
54+
"plt.xlabel('$l$'), plt.ylabel(r'$\\delta_l$', fontsize=16)\n",
55+
"\n",
56+
"plt.figure() \n",
57+
"plt.plot(theta*180/np.pi, sigma) # plot cross sections \n",
58+
"xtck = [0, 30, 60, 90, 120, 150, 180]\n",
59+
"plt.xticks(xtck, [repr(i) for i in xtck]) # custom ticks \n",
60+
"plt.xlabel(r'$\\theta$ (deg)')\n",
61+
"plt.ylabel(r'$\\sigma(\\theta)$ (a.u.)'), plt.semilogy()\n",
62+
"\n",
63+
"plt.show()"
64+
]
65+
},
66+
{
67+
"cell_type": "code",
68+
"execution_count": null,
69+
"metadata": {
70+
"collapsed": true
71+
},
72+
"outputs": [],
73+
"source": []
74+
}
75+
],
76+
"metadata": {
77+
"kernelspec": {
78+
"display_name": "Python 3",
79+
"language": "python",
80+
"name": "python3"
81+
},
82+
"language_info": {
83+
"codemirror_mode": {
84+
"name": "ipython",
85+
"version": 3
86+
},
87+
"file_extension": ".py",
88+
"mimetype": "text/x-python",
89+
"name": "python",
90+
"nbconvert_exporter": "python",
91+
"pygments_lexer": "ipython3",
92+
"version": "3.6.1"
93+
}
94+
},
95+
"nbformat": 4,
96+
"nbformat_minor": 2
97+
}

0 commit comments

Comments
 (0)