h_orbit_analytisch.py (11126B)
1 from scipy.constants import physical_constants 2 import numpy as np 3 from scipy.special import assoc_laguerre as lag #zugehörige Laguerre Polynome 4 from scipy.special import sph_harm as kugelfl #Kugelflächenfunktionen 5 import matplotlib.pyplot as plt 6 7 (a0,blub,bla)=physical_constants["Bohr radius"] 8 a0=0.529177210903 #in Angsröm 9 10 def fak(n): 11 #Beschreibung: Diese Funktion ermittelt zur natürliche Zahl n die zugehörige Fakultät (n!) 12 # Variablenname: Datentyp: Beschreibung: 13 #Input: n integer Zu Berechnende Fakultät. ACHTUNG: es gilt n>=0 und n element N 14 #Output: - integer n! 15 #Fehleroutput integer=-1 Falsche Eingabe für n 16 17 if (n<0) or (n-int(n)>0): 18 print("Upps, das hätte leider nicht passieren dürfen.\nn muss eine natürliche Zahl sein.") 19 return -1 20 elif (n==0) or (n==1): 21 return 1 22 else: 23 return n*fak(n-1) 24 25 def transform(x,y,z): 26 # Beschreibung: Diese Funktion Transfromiert Kartesische Koordinaten in Kugelkoordinaten 27 # Variablenname: Datentyp: Beschreibung: 28 # Input: x real kartesische x Koordinate 29 # y real kratesische y Koordinate 30 # z real kartesische z koordinate 31 # Output: r real Radius 32 # theta real Polarwinkel mit theta\in[0,pi] 33 # phi real Azimutwinkel mit phi\in[0,2*pi] 34 35 # Berechnung von r 36 r=np.sqrt(x**2+y**2+z**2) 37 #Problembehandlung: Koordinatenursprung 38 if(r==0): 39 phi=0. 40 theta=0. 41 else: 42 #Berechnung von theta 43 theta=np.arccos(z/r) 44 45 # Berechnung von phi 46 if (x>0)and(y>=0): 47 phi=np.arctan(y/x) 48 elif(x==0)and(y>0): 49 phi=np.pi/2 50 elif(x<0): 51 phi=np.arctan(y/x)+np.pi 52 elif(x==0)and(y<0): 53 phi=(3/2)*np.pi 54 elif(x>0)and(y<0): 55 phi=2*np.pi+np.arctan(y/x) 56 else:#if(x==0)and(y==0) 57 phi=0. 58 return r,theta,phi 59 60 def norm(Z,n,l): 61 #Beschreibung: Diese Funktion ermittelt den Normierungsfaktor der Wellenfunktion für ein Elektron um einen Atomkeern der Ladung Z*e 62 # Variablenname: Datentyp: Beschreibung: 63 #Input: Z integer Anzahl der Protonen im Atomkern 64 # n integer Hauptquantenzahl 65 # l integer Nebenquantenzahl 66 #Output: - real Normierungsfaktor der Wellenfunktion 67 68 return np.sqrt((((2*Z)/(n*a0))**3)*(fak(n-l-1)/(2*n*fak(n+l)))) 69 70 def radial(Z,n,l,r,normierung): 71 #Beschreibung: Diese Funktion ermittelt den Radialen Anteil R_nl(r) der Wellenfunktion für das Wasserstoffatom 72 # Variablenname: Datentyp: Beschreibung: 73 #Input Z integer Anzahl der Protonen im Atomkern 74 # n integer Hauptquantenzahl 75 # l integer Nebenquantenzahl 76 # r real Abstand des Elektrons zum Atomkern 77 # normierung real Normierungsfaktor der Wellenfunktion 78 #Output: - real Radialer Anteil R_nl(r) der Wellenfunktion für das Wasserstoffatom 79 80 rho=(2*Z*r)/(n*a0) 81 return normierung*np.exp(-rho/2)*(rho**l)*lag(rho,n-l-1,2*l+1) 82 83 def aufenthalt_radial(Z,n,l,r,normierung): 84 #Beschreibung: Diese Funktion ermittelt die Radiale Aufenthaltswahrscheinlichkeitsdichte des Wasserstoffatoms 85 # Variablenname: Datentyp: Beschreibung: 86 #Input Z integer Anzahl der Protonen im Atomkern 87 # n integer Hauptquantenzahl 88 # l integer Nebenquantenzahl 89 # r real Abstand des Elektrons zum Atomkern 90 # normierung real Normierungsfaktor der Wellenfunktion 91 #Output: - real Radiale Aufenthaltswahrscheinlichkeitsdichte des Wasserstoffatoms 92 93 return np.sqrt(radial(Z,n,l,r,normierung)**2)**2 94 95 def aufenthalt_welle(Z,n,l,m,r,normierung,phi,theta): 96 #Beschreibung: Diese Funktion ermittelt die Radiale Aufenthaltswahrscheinlichkeitsdichte des Wasserstoffatoms 97 # Variablenname: Datentyp: Beschreibung: 98 #Input Z integer Anzahl der Protonen im Atomkern 99 # n integer Hauptquantenzahl 100 # l integer Nebenquantenzahl 101 # m integer magnetische Quantenzahl des Drehimpulses 102 # r real Abstand des Elektrons zum Atomkern 103 # normierung real Normierungsfaktor der Wellenfunktion 104 # phi real Azimutwinkel mit phi\in[0,2*pi] 105 # theta real Polarwinkel mit theta\in[0,pi] 106 #Output: result real Aufenthaltswahrscheinlichkeitsdichte des Wasserstoffatoms 107 108 result=radial(Z,n,l,r,normierung)*kugelfl(m,l,phi,theta) 109 return np.sqrt(result.conjugate()*result)**2 110 111 112 def plotaufenthalt_radial(Z,n,l,aufl,name_nl,xmax): 113 #Beschreibung: Diese Funktion plottet die Radialenaufenthaltswahrscheinlichkeitsdichte*r**2, wobei nach Funktionsaufruf plt.show() ausgeführt werden muss 114 # Variablenname: Datentyp: Beschreibung: 115 #Input Z integer Anzahl der Protonen im Atomkern 116 # n integer,array Hauptquantenzahl, das Array muss dieselbe Größe wie l und name_nl haben 117 # l integer,array Nebenquantenzahl, das Array muss dieselbe Größe wie n und name_nl haben 118 # aufl integer Auflösung des Graphen 119 # name_nl characther,array Legende des Plots, das Array muss dieseleb Größe wie n und l haben 120 # xmax real maximaler Abstand des Elektrons zum Atomkern, wobei innerhald der Funktion xmax in xmax*a0 umgerechnet wird 121 #Output: - - Plotvorbereitung für die Radialenaufenthaltswahrscheinlichkeitsdichte*r**2 122 123 x=np.linspace(0,xmax,aufl) 124 x=x*a0 125 anz=np.shape(n)[0] #Anzahl der zu plottenden Orbitale 126 y=np.zeros((anz,aufl)) 127 r=x 128 for i in range(0,anz): 129 normi=norm(Z,n[i],l[i]) 130 y[i,:]=aufenthalt_radial(Z,n[i],l[i],r[:],normi)*r**2 131 132 # Plot 133 plt.figure() 134 for i in range(0,anz): 135 plt.plot(x,y[i,:],label=name_nl[i]) 136 plt.title("Radiale Aufenthaltswahrscheinlichkeitsdichte *r²") 137 plt.ylabel("|R_nl(r)|*r²") 138 plt.xlabel("r*a0 [Angström]") 139 plt.legend() 140 141 def plot2d_aufenthalt(Z,n,l,m,aufl,xmin,xmax,ymin,ymax,name): 142 #Beschreibung: Diese Funktion plottet die x-z Ebene der Aufenthaltswahrscheinlichkeitsdichte und die Aufenthaltswahrscheinlichkeitsdichte*r**2, 143 # wobei nach Funktionsaufruf plt.show() ausgeführt werden muss 144 # Variablenname: Datentyp: Beschreibung: 145 #Input Z integer Anzahl der Protonen im Atomkern 146 # n integer,array Hauptquantenzahl, das Array muss dieselbe Größe wie l und name haben 147 # l integer,array Nebenquantenzahl, das Array muss dieselbe Größe wie n und namehaben 148 # aufl integer Auflösung 149 # xmin,xmax real Gittergrenzen in x Richtung 150 # ymin,ymax real Gittergrenzen in y Richtung 151 # name character,array enthält die Spezifikation von n,l und m die geplottet werden und so in der Überschrift angezeigt werden, 152 # das Array muss dieselbe Größe wie n und l haben 153 #Output: - - Plotvorbereitung für die Aufenthaltswahrscheinlichkeitsdichte und die Aufenthaltswahrscheinlichkeitsdichte*r**2 154 155 #Definiere Gitter, Achtung y entspricht hier den Werten auf der z Achse 156 x=np.linspace(xmin,xmax,aufl)*a0 157 y=np.linspace(ymin,ymax,aufl)*a0 158 X,Y=np.meshgrid(x,y) 159 160 blub=np.shape(X)#Form von np.shape=(Anzahl y Werte, Anzahl x Werte, Anzahl z Werte) 161 162 #Bereite Transformation in Kugelkoordinaten vor 163 r=np.zeros(blub);theta=np.zeros(blub);phi=np.zeros(blub) 164 #Koordinatentransformation in Kugelkoordinaten 165 for i in range(0,blub[0]): #Schleife für z Werte 166 for j in range(0,blub[1]): #Schleife für x Werte 167 (r[i,j],theta[i,j],phi[i,j])=transform(X[i,j],0,Y[i,j]) 168 169 #Plottvorbereitungen 170 anz=np.shape(n)[0] #Anazhl der zu Plottenden Orbitale 171 172 for number in range(0,anz): 173 #Berechne Normaisierungsfaktor 174 normi=norm(Z,n[number],l[number]) 175 #Definiere Plot Array 176 orbit=np.zeros(np.shape(X)) 177 #Berechne Raumpunkte 178 orbit=np.real(aufenthalt_welle(Z,n[number],l[number],m[number],r,normi,phi,theta)) 179 180 #Plot 181 #Plot Aufenthaltswahrscheinlichkeitsdichte 182 plt.figure() 183 plt.imshow(orbit,cmap="gnuplot",extent=[x.min(),x.max(),y.min(),y.max()]) 184 plt.title("Aufenthaltswahrscheinlichkeitsdichte\nfür "+name[number]) 185 plt.xlabel("[Angström]") 186 plt.ylabel("[Angström]") 187 plt.colorbar() 188 #Plot Aufenthaltswahrscheinlichkeitsdichte*r² 189 plt.figure() 190 plt.imshow(orbit*r**2,cmap="gnuplot",extent=[x.min(),x.max(),y.min(),y.max()]) 191 plt.title("Aufenthaltswahrscheinlichkeitsdichte*r²\n für "+name[number]) 192 plt.xlabel("[Angström]") 193 plt.ylabel("[Angström]") 194 plt.colorbar() 195 196 197 #Plot der Radialenaufenthaltswahrscheinlichkeitsdichte*r**2 198 Z=1 #Anzahl der Protonen im Kern 199 n=np.array([1,2,2,]) 200 l=np.array([0,0,1]) 201 aufl=200 #Auflösung für den Plot der Radialenaufenthaltswahrscheinlichkeitsdichte*r**2 202 name_nl=np.array(["n=1,l=0","n=2,l=0","n=2,l=1"]) #Legende des Plots 203 xmax=20 204 205 plotaufenthalt_radial(Z,n,l,aufl,name_nl,xmax) 206 207 208 #Plot x-z Ebene 209 Z=1 210 n=np.array([1,2,2]) 211 l=np.array([0,0,1]) 212 m=np.array([0,0,0]) 213 #Definiere Gitterauflösung 214 aufl=403 215 #Definiere Gittergrenzen 216 xmin=-20 217 xmax=20 218 ymin=-20 219 ymax=20 220 name=np.array(["n=1, l=0, m=0","n=2, l=0, m=0","n=2, l=1, m=0"]) #Quantenzahlen, die in der Überschrift angezeigt werden 221 222 plot2d_aufenthalt(Z,n,l,m,aufl,xmin,xmax,ymin,ymax,name) 223 224 plt.show()