Interferenz an Einzel- und Doppelspalt

<<PhysicalConstants`

(* Festlegen von globalen Optionen für Plot-Funktion *)

SetOptions[Plot,BaseStyle->{FontSize->14,FontFamily->"Arial"}];

Huygens :

E=(Subscript[E, 0])((Subscript[r, 1]))(e^(i (k (Subscript[r, 1]) - ω t)))+(Subscript[E, 0])((Subscript[r, 2]))(e^(i (k (Subscript[r, 2]) - ω t)))+...+(Subscript[E, 0])((Subscript[r, N]))(e^(i (k (Subscript[r, 2N]) - ω t)))

mit (Subscript[E, 0])(r)=const. = (Subscript[E, 0]):

E=(Subscript[E, 0])(e^(-i ω t))(e^(-i k (Subscript[r, 1])))(1+(e^(i k ((Subscript[r, 2])-(Subscript[r, 1]) )))+...+(e^(i k ((Subscript[r, N])-(Subscript[r, 1]) ))))

(* Hecht - Optics *)

Phasendifferenz der einzelnen Quellen :

δ = kΔ mit Δ = d sinθ

Daher

δ=k ((Subscript[r, 2])-(Subscript[r, 1])), 2δ=k ((Subscript[r, 3])-(Subscript[r, 1])),3δ=k ((Subscript[r, 4])-(Subscript[r, 1])),...

und das Feld kann wie folgt beschrieben werden:

E=(Subscript[E, 0])(e^(-i ω t))(e^(-i k (Subscript[r, 1])))(1+(e^(i δ))+((((e^(i δ))))^2)...+((((e^(i δ))))^(N-1)))

Sum[(E^(j I δ)),{j,0,N-1}]

E=(Subscript[E, 0])(e^(-i ω t))(e^(-i k (Subscript[r, 1])))(((e^(i δ N))-1)/((e^(i δ))-1))

Distanz vom Zentrum des Spaltes zum Beobachtungspunkt:

R=(1/2)(N-1)d sin θ +(Subscript[r, 1])

k[λ_]:=2π/λ
r1[R_,N_,d_,θ_]:=R-(1/2)(N-1)d Sin[θ]
δ[λ_,d_,θ_]:=k[λ]d Sin[θ]
ω[λ_]:=SpeedOfLight k[λ]Second /Meter
Efield[λ_,R_,N_,a_,θ_,t_,E0_:1]:=E0 Exp[-I ω[λ]t]Exp[I k[λ] r1[R,N,a,θ]](((E^(I δ[λ,a,θ] N))-1)/((E^(I δ[λ,a,θ]))-1))+Conjugate[E0 Exp[-I ω[λ]t]Exp[I k[λ] r1[R,N,a,θ]](((E^(I δ[λ,a,θ] N))-1)/((E^(I δ[λ,a,θ]))-1))]

nm=1*^-9;μm=1*^-6;

Module[{Num,R,D,d,λ},
D=10μm;
Num=10;
(*Spalt mit fixer Breite D=N*d *)
d=D/Num;
R=1;
λ=500nm;
I0=((Abs[Limit[Efield[λ,R,Num,d,θ,1],θ->0]])^2);
Plot[((Abs[Efield[λ,R,Num,d,ArcSin[(y/(Sqrt[(R^2)+(y^2)]))],1]])^2)/I0,{y,-0.1,0.1},Frame->True,FrameLabel->{"y [m]","I/I0"},PlotRange->All,PlotPoints->30]]

Beugung am Einzelspalt

(* Funktion für Intensität *)
k[λ_]:=2π/λ
sinθ[y_,R_]:=(y/(Sqrt[(R^2)+(y^2)]))
β[λ_,d_,y_,R_]:=((k[λ] d)/2)sinθ[y,R]
Intensity[λ_,d_,y_,R_]= (((((Sin[β[λ,d,y,R]])/(β[λ,d,y,R]))))^2);

Plot[Intensity[λ nm,10 μm,y,R]/.{λ->500,R->1,d->10},{y,-0.2,0.2},Frame->True,FrameLabel->{"y [m]","I/(Subscript[I, 0])"}]

(* Intensität als Funktion der Distanz, der Spaltbreite und der Wellenlänge *)
Manipulate[
Plot[Intensity[λ nm,d μm,y,R],{y,-0.2,0.2},
PlotStyle->{Thickness->0.002,Blue},
Filling->Axis,
FillingStyle->Directive[Blue,Opacity[0.2]],
PlotRange->All,
Axes->None,
Frame->True,FrameStyle->Directive[Black,20],
FrameLabel->{"y [m]","I/I0","",""}],
{{R,1,"Distance [m]"},0.1,10},
{{d,10,"Slit width [μm]"},1,100},
{{λ,600,"Wavelength [nm]"},400,800},
LabelStyle->{Black,Italic,Large}]

Beugung am Doppelspalt

α[λ_,a_,y_,R_]:=((k[λ] a)/2)sinθ[y,R]
IntensityDS[λ_,d_,y_,R_,a_]= (((((Sin[β[λ,d,y,R]])/(β[λ,d,y,R]))))^2)(((Cos[α[λ,a,y,R]]))^2);

Manipulate[

Plot[IntensityDS[λ nm,d μm,y,R,a μm],{y,-0.2,0.2},PlotStyle->{Thickness->0.002,Red},
PlotRange->All,
PlotPoints->100,
Filling->Axis,
FillingStyle->Directive[Red,Opacity[0.2]],
Axes->None,
Frame->True,FrameStyle->Directive[Black,20],
FrameLabel->{"y [m]","I/I0","",""}],
{{R,1,"Distance [m]"},0.1,10},
{{d,5,"Slit width [μm]"},1,10},
{{λ,600,"Wavelength [nm]"},400,1000},
{{a,55,"Slit separation [μm]"},10,100},
LabelStyle->{Black,Italic,Large}]

Doppelspaltinterferenzmuster mit Einhüllender vom Einzelspalt

(* Display with dotted single slit envelope *)
Manipulate[Plot[{IntensityDS[λ nm,d μm,y,R,a μm],Intensity[λ nm,d μm,y,R]},{y,-0.2,0.2},PlotStyle->{{Thickness->0.002,Red},{Thickness->0.002,Dashed,Black}},
PlotRange->All,
PlotPoints->100,
Filling->{2->1,1->Axis},
FillingStyle->{None,Directive[Red,Opacity[0.2]]},
Axes->None,
Frame->True,FrameStyle->Directive[Black,20],
FrameLabel->{"y [m]","I/I0","",""}],
{{R,1,"Distance [m]"},0.1,10},
{{d,5,"Slit width [μm]"},1,10},
{{λ,600,"Wavelength [nm]"},400,1000},
{{a,55,"Slit separation [μm]"},10,100},
LabelStyle->{Black,Italic,Large}]

Beugungsgitter

IntensityDG[λ_,d_,y_,R_,a_,N_]:= (((((Sin[β[λ,d,y,R]])/(β[λ,d,y,R]))))^2)(((((Sin[N α[λ,a,y,R]])/(Sin[α[λ,a,y,R]]))))^2)

Manipulate[Plot[Evaluate[IntensityDG[λ nm,d μm,y,R,a μm,Num]],{y,-0.2,0.2},PlotStyle->{Thickness->0.002,Red},
PlotRange->All,
PlotPoints->1000,
Filling->Axis,
FillingStyle->Directive[Red,Opacity[0.2]],
Axes->None,
Frame->True,FrameStyle->Directive[Black,20],
FrameLabel->{"y [m]","I/I0","",""},
AspectRatio->0.5],
{{R,1,"Distance [m]"},0.1,10},
{{d,5,"Slit width [μm]"},1,10},
{{λ,600,"Wavelength [nm]"},400,1000},
{{a,55,"Slit separation [μm]"},10,100},
{{Num,2,"Number of slits"},2,6,1},
LabelStyle->{Black,Italic,Large},ContinuousAction->False]

Pseudo-Experiment zum Doppelspaltversuch

Generierung von Zufallszahlen mit einer bestimmten Verteilung :

z.B. Gaussverteilung :

gauss[x_,σ_]:=(1/(Sqrt[2π (σ^2)]))Exp[-(x^2)/(2(σ^2))];
gaussplot=Plot[gauss[x,1],{x,-5,5}]

(* Generiere Paare von Zufallszahlen: (x,y) *)
randompairs=Table[{RandomReal[{-4,4}],RandomReal[]},{i,1,1000}];
ListPlot[randompairs]

(* Wähle diejenigen x-Werte aus, für die y innerhalb der Gausskurve liegt : *)

SelectX[{x_,y_},σ_]:={x,y}/;(gauss[x,σ]>y)
SelectX[{x_,y_},σ_]:=NaN

(* Map-Funktion oder "/@" wendet eine Funktion auf alle Listeneinträge an:)

f/@{1,2,3,32}

(* Argumente können mit # als Platzhalter und anschliessendem "&" übergeben werden *)

f[a,#]&/@{1,12,23}

f[a,#[[2]],#[[1]]]&/@{{1,3},{2,12},{3,23}}

(* zurück zum Beispiel: *)

xvals=SelectX[#,1]&/@randompairs;

(* Eliminiere die ungültigen x-Werte *)

xvalsselected=Cases[xvals,_?ListQ]

ListPlot[xvalsselected];
Show[%,gaussplot]

Die x - Werte haben nun näherungsweise eine Gaussverteilung :

Histogram[xvalsselected[[All,1]],{-4,4,1}]

Pseudo Doppelspalt - Experiment :

(* Generiere Paare von Zufallszahlen: (zufälliger Ort y, zufällige normierte Intensität) *)
random2D[ym_,Nr_]:=Table[{RandomReal[{-ym,ym}],RandomReal[]},{i,1,Nr}]

random2D[0.1,20]

(* wenn Zufallszahl für Intensität kleiner als Intensität am Ort y, dann registriere dieses Event *)
(* Otherwise record a non-event as some number way outside the range (-123) *)
ClickFunction[λ_,d_,y_,R_,a_,randI_]:={y,randI}/;(IntensityDS[λ,d,y,R,a]>randI)
ClickFunction[λ_,d_,y_,R_,a_,randI_]:=NaN
counts[λ_,d_,ym_,R_,a_,Nr_]:= (ClickFunction[λ,d,#[[1]],R,a,#[[2]]])&/@random2D[ym,Nr]

counts[500nm,5μm,0.1,1,20μm,20]

countsselected[λ_,d_,ym_,R_,a_,Nr_]:=Cases[counts[λ,d,ym,R,a,Nr],_?ListQ]

countsselected[500nm,5μm,0.1,1,20μm,120]

Histogram[countsselected[500nm,5μm,0.1,1,20μm,10000][[All,1]],{-0.1,0.1,0.001}]

Fitten der Pseudodaten:

dataDS=BinCounts[countsselected[500nm,5μm,0.1,1,20μm,40000][[All,1]],{-0.1,0.1,0.001}]

Generiere eine Liste von (x, y) - Werten :

xaxis=Table[i,{i,-0.1+0.001/2,0.1,0.001}];

dataDSxy=Transpose[{xaxis,dataDS}];

Fit :

fitparams=FindFit[dataDSxy,Nmax IntensityDS[λfit,5μm,y,1,20μm],{{Nmax,30},{λfit,500 nm}},y]

(* Note the use of replacement rules here *)
fitplot=Plot[ Nmax IntensityDS[λfit,5μm,y,1,20μm]/.fitparams,{y,-0.1,0.1},PlotStyle->Directive[Red,Thickness[0.004]]];
dataplot=ListPlot[dataDSxy];

Show[fitplot,dataplot,Frame->True,FrameLabel->{"Position am Detektionsschirm [m]","Photonenzahl"},ImageSize->500,BaseStyle->{FontSize->14,FontFamily->"Arial"}]

Numerische Lösung einfacher Differentialgleichungen

Zerfallsgesetz :

DSolve[Nr'[t]==-λ Nr[t],Nr[t],t]

mit Randbedingung :

zerfall=DSolve[{Nr'[t]==-λ Nr[t],Nr[0]==N0},Nr[t],t]

Nr[t]/.First@zerfall/.N0->1

(E^(-1))//N

Plot[Nr[t]/.First@zerfall/.{N0->1,λ->0.1},{t,0,10},Frame->True,FrameLabel->{"t [s]","N"}]

Harmonischer Oscillator

solution0=DSolve[{y''[t]+y[t]==0,y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.First@solution0,{t,0,10},PlotRange->All,Frame->True,Axes->True,FrameLabel->{"t [s]","y[t]"}]

Gedämpfter harmonischer Oszillator

solution1=DSolve[{y''[t]+y'[t]/10+y[t]==0,y[0]==1,y'[0]==0},y[t],t]
solution1=DSolve[{y''[t]+y'[t]/10.+y[t]==0,y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.First@solution1,{t,0,10},PlotRange->All,Frame->True,Axes->True,FrameLabel->{"t [s]","y[t]"}]

Getriebener gedämpfter harmonischer Oszillator

solution2=DSolve[{y''[t]+0.1y'[t]+y[t]==a Cos[0.3 t],y[0]==1,y'[0]==0},y[t],t]

Plot[y[t]/.solution2/.a->0.1,{t,0,130},PlotRange->All,Frame->True,FrameLabel->{"t [s]","y[t]"}]

DensityPlot[y[t]/.solution2,{t,0,50},{a,0,1},PlotPoints->40,FrameLabel->{"t","a"},BaseStyle->{FontSize->14,FontFamily->"Arial"}]

Zusätzliche "Komplikation" durch zeitabhängige Dämpfung :

solution3=DSolve[{y''[t]+0.1Cos[0.5 t]y'[t]+y[t]==0.1 Cos[0.3 t],y[0]==1,y'[0]==0},y[t],t]

DSolve kann diese Art von Differentialgleichung nicht mehr lösen, numerisch funktioniert das jedoch mit NDSolve.

solution3=NDSolve[{y''[t]+0.1Cos[0.5 t]y'[t]+y[t]==0.1Cos[0.3 t],y[0]==1,y'[0]==0},y[t],{t,0,50}]

Plot[y[t]/.solution3,{t,0,50},PlotRange->All,Frame->True]