Beispiel Brechungsindex
Ausarbeitung Helmut Hörner 8850092

Polarisierung eines angetriebenen Oszillators

α=(( ((ωp)^2))/((((ω0)^2)-(ω^2))+I ω/τ)) ; (* Amplitude *)

Et = E0 Cos[ω t]; (* Treibendes E-Feld *)

P[t]=FullSimplify[ϵ0 α Et,ComplexityFunction->LeafCount] (* Polarierung *)

χ=ϵ-1; (* definiere Suszeptibilität *)

Ermittle n² aus der Beziehung P[t] = ϵ0*χ*E[t]

nSquared=Apply[List,Flatten[FullSimplify[Solve[P[t]==ϵ0 χ Et ,ϵ],ComplexityFunction->LeafCount]][[1]]][[2]] (* Berechne n² *)

Daraus folgt bei Substitution ωp->(Sqrt[((N (q^2))/(ϵ0 m))]) und τ->(1/γ) 
die allgemeine Formel für n² gem. Demtröder 2 (8.32):

nSquared2=FullSimplify[nSquared /. { ωp->(Sqrt[((N (q^2))/(ϵ0 m))]),τ->(1/γ)},ComplexityFunction->LeafCount]

Für (n-1)<<1 gilt (n²-1) -> 2(n-1). 
Daraus folgt die Näherungsformel für dünnes Medium (zwei Varianten)
siehe Demtröder 2, (8.12a)

(* Abhängig von ωp, τ *)
nDuenn=Apply[List,Flatten[FullSimplify[Solve[2(n-1)==nSquared-1,n],ComplexityFunction->LeafCount]][[1]]][[2]] 
(* Abhängig von N, q, m, γ *)
nDuenn2=Apply[List,Flatten[FullSimplify[Solve[2(n-1)==nSquared2-1,n],ComplexityFunction->LeafCount]][[1]]][[2]]

Dispersionsrelationen: n = Re[nDuenn] - Im[nDuenn] (zwei Varianten)
siehe Demtröder 2, (8.21a.b)

NurRealeParameter={Element[{N,q,ω,m,ϵ0,γ,ω0, τ, ωp}, Reals],N>0,q>0,ω>0,m>0,ϵ0>0,γ>0,ω0>0, τ>0, ωp>0};

(* Abhängig von ωp, τ *)
nRe = FullSimplify[Re[ComplexExpand[nDuenn]],Assumptions->NurRealeParameter]
nIm=FullSimplify[Im[ComplexExpand[nDuenn]],Assumptions->NurRealeParameter]
(* Abhängig von N, q, m, γ *)
nRe2 = FullSimplify[Re[ComplexExpand[nDuenn2]],Assumptions->NurRealeParameter]
nIm2=FullSimplify[Im[ComplexExpand[nDuenn2]],Assumptions->NurRealeParameter]

Bei Metallen: Elektronen frei beweglich; daher k=0 und daher ω0=k/m=0. 
Daraus folgt die allgemein gültige n²-Formel für Metalle (drei Varianten):

nMetallSquared = Limit[nSquared,ω0->0] (* Abhängig von ωp, τ *)
nMetallSquared2 = Limit[nSquared2,ω0->0] (* Abhängig von N, q, m, γ *)
nMetallSquared3 = FullSimplify[Limit[nSquared2,ω0->0]/.{ ((N (q^2))/m)->  σ γ} /.{γ-> (1/τ)},ComplexityFunction->LeafCount](* Abhängig von σ, τ *)

Plot Re[n] und Im[n] für Silber

SilberDaten = {τ->1.13997 (10^(-14)), ωp-> 1.23026 (10^16), ω0->0};

Plot[{Re[(Sqrt[nMetallSquared]) /. SilberDaten ],Im[(Sqrt[nMetallSquared]) /. SilberDaten ]},{ω,0.9 (10^16), 1.5 (10^16)}, PlotLegends->{"Re(ω)","Im(ω)"}]

Manipulate[Plot[{nRe /. {τ->Tau,ωp-> Plasmafreq, ω0->Resonanzfreq} ,nIm /. {τ->Tau,ωp-> Plasmafreq, ω0->Resonanzfreq}},{ω,0.8 (10^16), 1.5 (10^16)}],{Tau,0,5 (10^(-14))},{Plasmafreq,0,(10^17)}, {Resonanzfreq,0,(10^17)}]