SHARE
TWEET

Untitled

a guest Oct 21st, 2019 73 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])
  29.     for i in range(blocks):
  30.         blocked[i]= float((potential(x(i+1,start,end,blocks))+potential(x(i,start,end,blocks)))/2)
  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):
  47.         k1=sp.sqrt(m/hsq*(E-V[i]))
  48.         k2=sp.sqrt(m/hsq*(E-V[i+1]))
  49.         C=(k1*m2+k2*m1)/(2*k1*m2) +0j
  50.         D=(k1*m2-k2*m1)/(2*k1*m2) +0j
  51.  
  52.         x0=x(i+1,start,end,blocks)
  53.         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)]))
  54.  
  55.         if i>135 and i <155 :
  56.             print('\n i='+str(i) +', x=' +str(x0) + ', k1=' + str(k1)+', k2=' + str(k2) )
  57.             print(1/Abs(M[i,0,0])**2+ Abs(M[i,1,0]/M[i,0,0])**2)
  58.             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)])
  59.             print(M[i,])
  60.        
  61.     print((1/Abs(M[blocks-1,0,0])**2+ Abs(M[blocks-1,1,0])**2/Abs(M[blocks-1,0,0])**2))
  62.    
  63.     for i in range(blocks):
  64.         trans[i]=Abs(sp.sqrt(m/hsq*(E-V[i]))/(k0*M[i,0,0]**2))
  65.         refl[i]=Abs(M[i,1,0]/M[i,0,0])**2
  66.        
  67.  
  68.     print('transmission= ' +str(trans[blocks-1]))
  69.     print('reflection= ' +str(refl[blocks-1]))
  70.    
  71.     grid = np.linspace(start,end,blocks)
  72.     gs= gridspec.GridSpec(2,1,height_ratios=[2,1])
  73.  
  74.     a0=plt.subplot(gs[0])
  75.     a0.plot(grid, trans, label='transmission')
  76.     a0.plot(grid, refl, label='reflection')
  77.     a0.legend(loc='upper left')
  78.     #plt.title('Propability density')
  79.     plt.xlabel('x [' + u"\u212B" +"]")
  80.     #plt.ylabel(u"\u03A8"+' ^2')
  81.     a0.xaxis.set_label_coords(1.05, -0.025)
  82.  
  83.     a1=plt.subplot(gs[1])
  84.     plt.plot(grid, V)
  85.     plt.title('Potential')
  86.     plt.ylabel('V(x)')
  87.     plt.xlabel('x [' + u"\u212B" +"]")
  88.     a1.xaxis.set_label_coords(1.05, -0.025)
  89.    
  90.     plt.subplots_adjust(hspace=0.32)
  91.     plt.savefig(potential.__name__+'.png', dpi=300)
  92.     plt.close()
  93.     return;
  94.    
  95. transmission(Vparcut,a,b,N)
  96. #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