import numpy  # Enthält Datentypen für numerisch aufwendige Berechnungen, ist ein Teil von
import scipy  # welches fast alles enthält, was man als (numerischer) Physiker braucht, z.B.
import scipy.constants  # Zahlenwerte für Naturkonstanten und
import scipy.integrate  # Methoden zu numerischen Integration

h = scipy.constants.physical_constants['Planck constant'][0] # So erhält man Werte...
type(scipy.constants.physical_constants) # Ist ein dict, dessen Schlüssel man mit .keys() erfrägt, und
scipy.constants.physical_constants['Planck constant'] # dessen Werte als (Zahl, Einheit, Unsicherheit)-Tupel enthalten sind

%latex
\noindent Wir wollen als erstes $\vert \Phi(p_x) \vert^2$ normieren, also müssen wir das Integral
\[ I_\mathrm{Normierung} = \int_{-\pi/2}^{\pi/2} \left( \frac{\sin(\frac{\pi d}{\lambda} \sin(\phi))}{\frac{\pi d}{\lambda} \sin(\phi)} \right)^2 d\phi
\]
berechnen.

pi = scipy.constants.pi  # Benutze diese Definition, da potentiell schneller
Lambda = 5.6e-7  # Meter, wir bleiben in SI-Einheiten, damit wir dadurch kein Kopfzerbrechen haben. Außerdem ist lambda ein Python-Keyword, deswegen besser Lambda verwenden
d = 4e-5 # Meter
a = pi * d / Lambda
# Definiere den ersten Integranden:
Phipx = lambda x: (numpy.sin(a*numpy.sin(x)) / (a*numpy.sin(x)))**2 # Numpy enthält auch spezielle Funktionen, die schneller rechnen

# Nun kann das Integral berechnet werden, und da dieses symmetrisch ist, sparen wir uns die halbe Rechenzeit:
zwischenwert = scipy.integrate.quad(Phipx, 0, pi/2)
zwischenwert # (Ergebnis, Genauigkeit)
normierung = zwischenwert[0]*2
normierung

%latex
\noindent Nun können wir den Erwartungswert $<p_x^2>$ berechnen:
\[ \frac{1}{I_\mathrm{Normierung}}  \int_{-\pi/2}^{\pi/2} \left( \frac{h}{\lambda} \frac{\sin(\frac{\pi d}{\lambda} \sin(\phi))}{\frac{\pi d}{\lambda} } \right)^2 d\phi
\]
berechnen.

# Definiere den neuen Integranden:
px2int = lambda x : ( h/Lambda *  numpy.sin(a*numpy.sin(x)) / a)**2
# Auch dieser ist symmetrisch
zwischenwert = scipy.integrate.quad(px2int, 0, pi/2)
px2 = zwischenwert[0]*2 / normierung
px2

# Nun kommen wir zum Ergebnis:
deltapx = numpy.sqrt(px2)
deltapx
# Und wenn man h als Faktor herauszieht:
deltapx/h

# Und schließlich:
produkt = d*deltapx
produkt
#bzw.
produkt / h

# Schlussbemerkung: h im Integral zu belassen, war zwar anschaulich, aber für die Genauigkeit schlecht. Besser ist, alles ohne h zu berechnen, da uns sowieso nur Vielfache von h interessieren.
︠9768af58-3914-4131-959b-f8a74f51489d︠