SHARE
TWEET

Untitled

a guest Oct 21st, 2019 81 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import scipy as sp
  3. import matplotlib.pyplot as plt
  4. from matplotlib import gridspec
  5. from sympy import *
  6. from functions import *
  7.  
  8. def Vparcut(x):
  9.     return max(0., 0.4 - 0.001 * x**2);
  10.  
  11. hsq = 3.81 # h^2/2m_e
  12.  
  13. a0 = 70.    #width of potential well
  14. d = 200.    #total width
  15.  
  16. V0 = 0.3    #well's depth
  17. m = 0.3 #effective mass
  18.  
  19. a = -115.   #starting point
  20. b = 115.        #ending point
  21. N = 300     #number of blocks
  22.  
  23.  
  24. def x(i, start, end, blocks):
  25.     return float(start + i * (end-start) / blocks);
  26.  
  27. def blocked(potential, start, end, blocks):
  28.     blocked = np.empty([blocks]) # nazywanie zmiennej tak samo jak funkcje to zla praktyka
  29.     for i in range(blocks):
  30.         blocked[i] = float((potential(x(i+1, start, end, blocks)) + potential(x(i, start, end, blocks))) / 2.)# fix potencjalnego buga - moze rzutowac na inta
  31.     return blocked;
  32.    
  33. def transmission(potential, start, end, blocks):
  34.     E = 0.35
  35.     m2 = 1
  36.     m1 = 1
  37.     V = blocked(potential, start, end, blocks)
  38.    
  39.     M = np.empty([blocks, 2, 2], dtype=np.complex)
  40.     M[0,] = np.identity(2)
  41.  
  42.     k0 = sp.sqrt(m / hsq * (E-V[0]))
  43.     trans = np.empty([blocks])
  44.     refl = np.empty([blocks])
  45.    
  46.     for i in range(blocks - 1): # na pewno -1? to nie sugestia zmiany tylko upewnienia sie, i te wszystkie nizej roznania polecam sprawdzic z naciskiem na +1, m2 m1 itd, tak w ogole nie polecam nazywac tak zmiennych, bardzo latwo sie jebnac
  47.         k1 = sp.sqrt(m / hsq * (E - V[i]))
  48.         k2 = sp.sqrt(m / hsq * (E - V[i + 1]))
  49.      # (k1 * m2 + k2 * m1) dalabym w inna zmienna bo sie powtarza, a moze to zle i blad?
  50.         C = (k1 * m2 + k2 * m1) / (2 * k1 * m2) + 0j
  51.         D = (k1 * m2 - k2 * m1) / (2 * k1 * m2) + 0j
  52.  
  53.         x0 = x(i + 1, start, end, blocks)
  54.     # potwora nizej to wgl rozbic na zmienne
  55.         M[i + 1,] = np.matmul(M[i], ([C * sp.exp(1j * (k2 - k1) * x0), D * sp.exp(-1j * (k2 + k1) * x0)], [D * sp.exp(1j*(k2 + k1) * x0), C * sp.exp(-1j * (k2 - k1) * x0)]))
  56.  
  57.         if i > 135 and i < 155 :
  58.             print('\n i=' + str(i) + ', x=' + str(x0) + ', k1=' + str(k1)+', k2=' + str(k2))
  59.             print(1 / Abs(M[i,0,0])**2 + Abs(M[i,1,0] / M[i,0,0])**2)
  60.       # to nizej tez rozbic na zmienne
  61.             print([C * sp.exp(1j * (k2 - k1) * x0), D * sp.exp(-1j * (k2 + k1) * x0)], [D * sp.exp(1j * (k2 + k1) * x0) , C * sp.exp(-1j * (k2 - k1) * x0)])
  62.             print(M[i,])
  63.        
  64.     print((1 / Abs(M[blocks-1,0,0])**2 + Abs(M[blocks-1,1,0])**2 / Abs(M[blocks-1,0,0])**2))
  65.    
  66.     for i in range(blocks):
  67.         trans[i] = Abs(sp.sqrt(m / hsq * (E - V[i])) / (k0 * M[i,0,0]**2))
  68.         refl[i] = Abs(M[i,1,0] / M[i,0,0])**2
  69.        
  70.  
  71.     print('transmission= ' + str(trans[blocks - 1]))
  72.     print('reflection= ' + str(refl[blocks - 1]))
  73.    
  74.     grid = np.linspace(start, end, blocks)
  75.     gs = gridspec.GridSpec(2, 1, height_ratios=[2,1])
  76.  
  77.     a0 = plt.subplot(gs[0])
  78.     a0.plot(grid, trans, label='transmission')
  79.     a0.plot(grid, refl, label='reflection')
  80.     a0.legend(loc='upper left')
  81.     #plt.title('Propability density')
  82.     plt.xlabel('x [' + u"\u212B" +"]")
  83.     #plt.ylabel(u"\u03A8"+' ^2')
  84.     a0.xaxis.set_label_coords(1.05, -0.025)
  85.  
  86.     a1 = plt.subplot(gs[1])
  87.     plt.plot(grid, V)
  88.     plt.title('Potential')
  89.     plt.ylabel('V(x)')
  90.     plt.xlabel('x [' + u"\u212B" +"]")
  91.     a1.xaxis.set_label_coords(1.05, -0.025)
  92.    
  93.     plt.subplots_adjust(hspace=0.32)
  94.     plt.savefig(potential.__name__+'.png', dpi=300)
  95.     plt.close()
  96.     return;
  97.    
  98. transmission(Vparcut, a, b, N)
  99. #transmission(Vsq,a,b,N)
RAW Paste Data
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