• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Mar 28th, 2020 99 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. from pprint import pprint
2. import numpy as np
3. import scipy.constants as c
4. from numpy import exp
5. from math import sin, cos, atan, acos, sqrt, pi
7. i=complex(0,1)
8. def complex_round(x,n):
9.     ##print(type(x))
10.     reel=x.real
11.     imaaag=x.imag
12.     #print(x)
13.     #print(complex(round(reel,n), round(imaaag,n)))
14.     return complex(round(reel,n), round(imaaag,n))
15. def spherical_to_cartesian(r,phi,theta):
16.     x=r*sin(theta)*cos(phi)
17.     y=r*sin(phi)*sin(theta)
18.     z=r*cos(phi)
19.     if abs(x)< 10e-14: x=0
20.     if abs(y)< 10e-14: y=0
21.     if abs(z)< 10e-14: z=0
22.     #print(x,y,z)
23.     return x,y,z
24.     pass
25. def sci_round(x,n):
26.     return float(('{:0.'+str(n)+'e}').format(x))
27. def cartesian_to_spherical(x, y, z):
28.     r=round(sqrt(x**2+y**2+z**2),4)
29.     theta=round(atan(y/x),4) if x!=0  else (1.5708*(y/abs(y)) if y!=0 else  0)
30.     phi=round(acos(z/r),4) if r !=0 else 0
31.     #print(r,phi,theta)
32.     return r,phi,theta
33.     pass
34.
35.
38.     constant = 1
39.     if n == 1 and l == 0:
40.         return complex_round(2 * constant * exp(-(r / a)),5)
41.
42.     elif n == 2 and l == 0:
43.         return complex_round((1 / sqrt(2)) * constant * (1 - (r / (2 * a))) * exp(-(r / (2 * a))),5)
44.
45.     elif n == 2 and l == 1:
46.         return complex_round((1 / sqrt(24)) * constant * (r / a) * exp(-(r / (2 * a))),5)
47.
48.     elif n == 3 and l == 0:
49.         return complex_round((2 / (81 * sqrt(3))) * constant * (27 - 18 * (r / a) + 2 * (r / a) ** 2) * exp(-(r / (3 * a))),5)
50.
51.     elif n == 3 and l == 1:
52.         return complex_round((8 / (27 * sqrt(6))) * constant * (1 - (r / (6 * a))) * (r / a) * exp(-(r / (3 * a))),5)
53.
54.     elif n == 3 and l == 2:
55.         return complex_round((4 / (81 * sqrt(30))) * constant * ((r / a) ** 2) * exp(-(r / (3 * a))),5)
56.
57.     elif n == 4 and l == 0:
58.         return complex_round((1 / 4) * constant * (1 - (3 / 4) * (r / a) + (1 / 8) * (r / a) ** 2 - (1 / 192) * (r / a) ** 3) * exp(
59.             -(r / (4 * a))),5)
60.
61.     elif n == 4 and l == 1:
62.         return complex_round((sqrt(5) / (16 * sqrt(3))) * constant * (1 - (1 / 4) * (r / a) + (1 / 80) * (r / a) ** 2) * exp(
63.             -(r / (4 * a))),5)
64.
65.     elif n == 4 and l == 2:
66.         return complex_round((1 / (64 * sqrt(5))) * constant * (r / a) ** 2 * (1 - (1 / 12) * (r / a)) * exp(-(r / (4 * a))),5)
67.
68.     elif n == 4 and l == 3:
69.         return complex_round((1 / (768 * sqrt(35))) * constant * (r / a) ** 3 * exp(-(r / (4 * a))),5)
70.
71. def angular_wave_func(m, l, theta, phi):
72.     if l == 0 and m == 0:
73.         return complex_round(sqrt(1 / (4 * pi)),5)
74.     elif l == 1 and m == 1:
75.         #print((-sqrt(3 / (8 * pi)) * sin(theta) * exp(i * phi),5))
76.         return complex_round(-sqrt(3 / (8 * pi)) * sin(theta) * exp(i * phi),5)
77.     elif l == 1 and m == 0:
78.         return complex_round(sqrt(3 / (4 * pi)) * cos(theta),5)
79.         pass
80.     elif l == 1 and m == -1:
81.         return complex_round(sqrt(3 / 8 / pi) * sin(theta) * exp(-i * phi),5)
82.         pass
83.     elif l == 2 and m == 2:
84.         return complex_round(sqrt(15 / 32 / pi) * (sin(theta) ** 2) * exp(i * 2 * phi),5)
85.         pass
86.     elif l == 2 and m == 1:
87.         return complex_round((-sqrt(15 / 8 / pi)) * cos(theta) * sin(theta) * exp(i * phi),5)
88.         pass
89.     elif l == 2 and m == 0:
90.         return complex_round(sqrt(5 / 16 / pi) * (3 * (cos(theta) ** 2) - 1),5)
91.         pass
92.     elif l == 2 and m == -1:
93.         return complex_round(sqrt(15 / 8 / pi) * cos(theta) * sin(theta) * exp(-i * phi),5)
94.         pass
95.     elif l == 2 and m == -2:
96.         return complex_round(sqrt(15 / 32 / pi) * (sin(theta) ** 2) * exp(-i * 2 * phi),5)
97.         pass
98.     elif l == 3 and m == 3:
99.         return complex_round(-sqrt(35 / 64 / pi) * (sin(theta) ** 3) * exp(i * 3 * phi),5)
100.         pass
101.     elif l == 3 and m == 2:
102.         return complex_round(sqrt(105 / 32 / pi) * cos(theta) * (sin(theta) ** 2) * exp(i * phi),5)
103.         pass
104.     elif l == 3 and m == 1:
105.         return complex_round(-sqrt(21 / 64 / pi) * sin(theta) * (5 * (cos(theta) ** 2) - 1) * exp(i * phi),5)
106.         pass
107.     elif l == 3 and m == 0:
108.         return complex_round(sqrt(7 / 16 / pi) * (5 * (cos(theta) ** 3) - 3 * cos(theta)),5)
109.         pass
110.     elif l == 3 and m == -1:
111.         return complex_round(sqrt(21 / 64 / pi) * sin(theta) * (5 * (cos(theta) ** 2) - 1) * exp(-i * phi),5)
112.         pass
113.     elif l == 3 and m == -2:
114.         return complex_round(sqrt(105 / 32 / pi) * cos(theta) * (sin(theta) ** 2) * exp(-i * 2 * phi),5)
115.         pass
116.     elif l == 3 and m == -3:
117.         return complex_round(sqrt(35 / 64 / pi) * (sin(theta) ** 3) * exp(-i * 3 * phi),5)
118.         pass
119.     else:
120.         return
121.     pass
122. # this function accounts for linear combination of angular wave functions
123. # (week 2 Lesson 1 slide 27)
124. def lincomb_angular_wave_func(m, l, theta, phi):
125.     if m<0:
126.         return (1j / np.sqrt(2)) * (
127.                 angular_wave_func(m, l, theta, phi) - ((-1) ** float(m)) * angular_wave_func(-m, l, theta, phi))
128.         pass
129.     elif m>0:
130.         return (1 / np.sqrt(2)) * (
131.                     angular_wave_func(-m, l, theta, phi) + ((-1) **  float(m)) * angular_wave_func(m, l, theta, phi))
132.         pass
133.     elif m==0:
134.         return angular_wave_func(m, l, theta, phi)
135.
136. def mgrid3d(xstart, xend, xpoints,
137.             ystart, yend, ypoints,
138.             zstart, zend, zpoints):
139.     xr = []
140.     yr = []
141.     zr = []
142.     xval = xstart
143.     xcount = 0
144.
145.     # calculate the step size for each axis
146.     xstep = (abs(xend) + abs(xstart)) / (xpoints - 1)
147.     ystep = (abs(yend) + abs(ystart)) / (ypoints - 1)
148.     zstep = 0
149.     if zpoints != 0 or abs(abs(zend) - abs(zstart)) != 0:
150.         zstep = abs(abs(zend) - abs(zstart)) / (zpoints - 1)
151.     if zend < 0 or zstart < 0:
152.         zstep = abs(abs(zend) + abs(zstart)) / (zpoints - 1)
153.     #print(zend, zstart, zstep)
154.     while xcount < xpoints:
155.         # initialize y points
156.         yval = ystart
157.         ycount = 0
158.         # create temporary variable to store x, y and z list
159.         xrow = []
160.         yrow = []
161.         zrow = []
162.
163.         while ycount < ypoints:
164.             # initialize z points
165.             zval = zstart
166.             zcount = 0
167.             # create temporary variable to store the inner x, y, and z list
168.             inner_xrow = []
169.             inner_yrow = []
170.             inner_zrow = []
171.
172.             while zcount < zpoints:
173.                 # add the points x, y, and z to the inner most list
174.                 inner_xrow.append(xval)
175.                 inner_yrow.append(yval)
176.                 inner_zrow.append(zval)
177.
178.                 # increase z point
179.                 zval += zstep
180.                 zcount += 1
181.             # add the inner most lists to the second lists
182.             xrow.append(inner_xrow)
183.             yrow.append(inner_yrow)
184.             zrow.append(inner_zrow)
185.
186.             # increase y point
187.             yval += ystep
188.             ycount += 1
189.         # add the second lists to the returned lists
190.         xr.append(xrow)
191.         yr.append(yrow)
192.         zr.append(zrow)
193.
194.         # increase x point
195.         xval += xstep
196.         xcount += 1
197.     return (xr, yr, zr)
198.
199.
200. def mgrid2d(xstart, xend, xpoints, ystart, yend, ypoints):
201.     # initialize a list to store the grid points that will be returned
202.     xr = []
203.     yr = []
204.
205.     # initialize the first x value
206.     xval = xstart
207.
208.     # initialize a variable to count the number of x points
209.     xcount = 0
210.
211.     # Calculate the step size for x and y, replace None with the right expression
212.     xstep = (abs(xend) + abs(xstart)) / (xpoints - 1)
213.     ystep = (abs(yend) + abs(ystart)) / (ypoints - 1)
214.
215.     while xcount < xpoints:
216.         # initialize the first y value, replace None with the right value
217.         yval = ystart
218.
219.         # initialize the variable to count the number of y points, replace None with the right value
220.         ycount = 0
221.
222.         # initialize temporary lists
223.         xrow = []
224.         yrow = []
225.
226.         while ycount < ypoints:
227.             # write code to add the current x value to the temporary x list
228.             xrow.append(xval)
229.
230.             # write code to add the current y value to the temporary y list
231.             yrow.append(yval)
232.
233.             # increase the y value by the step size in y
234.             yval += ystep
235.
236.             # increase the counter for the number of y points
237.             ycount += 1
238.         # Add the temporary x list to the final x list
239.         xr.append(xrow)
240.
241.         # Add the temporary y list to the final y list
242.         yr.append(yrow)
243.
244.         # increase x value by the step size in x
245.         xval += xstep
246.
247.         # increase the counter for the number of x points
248.         xcount += 1
249.
250.     return xr, yr
251.
252. def hydrogen_wave_func(n,l,m,roa,Nx,Ny,Nz):
253.     #vectorizing functions to be used on arrays
254.     vec_cart_sphere = np.vectorize(cartesian_to_spherical)
256.     vec_linearcombi_angular_wave_func = np.vectorize(lincomb_angular_wave_func)
257.     array_complex_round = np.vectorize(complex_round)
258.     array_round = np.vectorize(round)
259.
260.     xstep = roa*2/(Nx-1)
261.     xrange = [round(-roa+xstep*i,5) for i in range(Nx)]
262.     ystep = roa*2/(Ny-1)
263.     yrange = [round(-roa+ystep*i ,5)for i in range(Ny)]
264.     zstep = roa*2/(Nz-1)
265.     zrange = [round(-roa+zstep*i,5) for i in range(Nz)]
266.
267.     xx = [[[x for z in zrange] for y in yrange] for x in xrange]
268.     yy = [[[y for z in zrange] for y in yrange] for x in xrange]
269.     zz = [[[z for z in zrange] for y in yrange] for x in xrange]
270.
271.     #convert array to spherical
272.     r,theta,phi = vec_cart_sphere(xx,yy,zz)
273.
275.