G2A Many GEOs
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
  6. a=c.physical_constants['Bohr radius'][0]
  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.  
  36. def radial_wave_func(n, l, r):
  37.     a = c.physical_constants['Bohr radius'][0]
  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)
  255.     vec_radial_wave_func = np.vectorize(radial_wave_func)
  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.  
  274.     a = c.physical_constants['Bohr radius'][0]
  275.  
  276.     #get radial function output
  277.     radial = vec_radial_wave_func(n,l,r*a)
  278.  
  279.     #get angular function output, linearly combined
  280.     angular = vec_linearcombi_angular_wave_func(m,l,theta,phi)
  281.  
  282.     #get the square of the magnitude of the wavefunction, rounded to 5
  283.     mag = array_round(abs(np.multiply(radial, angular))**2,5)
  284.  
  285.     #convert nparray back to python list
  286.     return np.round((xx,yy,zz,mag.tolist()),5)
RAW Paste Data
Ledger Nano X - The secure hardware wallet
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top