Mathematica Intro : 3

- Reminders of some basics
- Wavepacket evolution
- Black body laws

(* Choose this depending whether presenting/printing *)
plotsize=750

Reminders

Tips on using Mathematica

Run pieces of code by pressing Shift+Enter with the cursor on the code.

Use curly brackets to define lists of anything:    {1,2,3,4}      {"One","Two","Three"}      {(x^2),12,"Quantum"}

Use square brackets when calling and defining functions

Use the esc key for special characters eg esc-p-esc gives π

Immediate/delayed assignment

(* Immediate assignment *)
variable=53;

(* Delayed assignment *)
func[x_]:=(x^2)

variable

Lists

primeList=Table[Prime[n],{n,1,200}]

Plots

plotParams={Frame->True,FrameStyle->Directive[16,Black],FrameLabel->{"Which prime","The prime","The first 200 prime numbers",""},ImageSize->450};

primePlot=ListPlot[primeList,plotParams]

fullPath="/Users/peterleek/Desktop/";
Export[fullPath<>"primePlot.pdf",primePlot]

Replacement rules

x

x/.x->y

x

<<PhysicalConstants`

SpeedOfLight

kb T

kb T/ElectronCharge/.{kb->BoltzmannConstant,T->300 Kelvin}

kb T

Help system

Use the Help system!!
Just press F1 on a function name, and the correct Help page will appear.

E.g.

Integrate
Sum
Fourier
Expand

Wavepacket Motion

Clear[x,t,k,ω]

g[k_]:=(((((a^2)/(2π))))^(1/4))Exp[-(a^2)(((k-k0))^2)/4]
Ψ[x_,t_,ω_]:=(1/(Sqrt[2π]))(∫_(-∞)^∞)g[k]Exp[I(k x -ω t)]\[DifferentialD]k

Ψ[x,0,(ℏ/2m)(k^2)]

Ψt0[a_,k0_,x_,t_]:=(1/(((2 π))^(3/4)))((((a^2)))^(1/4))((2 (E^(x (I k0-(x/(a^2))))) (Sqrt[π]))/(Sqrt[(a^2)]))

Ψ[x,t,(ℏ/2m)(k^2)]

Ψt[a_,k0_,x_,t_,m_]:=(1/(((2 π))^(3/4)))((((a^2)))^(1/4))((2 (E^(-((I (a^2) k0 (k0 m t-2 x)+2 (x^2))/(2 ((a^2)+2 I m t))))) (Sqrt[π]))/(Sqrt[(a^2)+2 I m t]))

Manipulate[
Plot[Re[Ψt[a,k0,x,t,m]]/.{a->1,k0->5,m->1},{x,-5,100},PlotRange->{All,{-1,1}},PlotPoints->200,ImageSize->plotsize,PlotStyle->Directive[Blue,Thickness[0.001]],Frame->True,FrameStyle->Directive[Black,20],FrameLabel->{"x","Ψ[x,t]","",""},Axes->False],
{{t,1,"Time"},0,10},LabelStyle->{Black,Large,Italic}]

Ψ[x,t,c k]

Ψtmassless[a_,k0_,x_,t_]:=(1/(((2 π))^(3/4)))((((a^2)))^(1/4))((2 (E^(((5 t-x) (-I (a^2) k0-5 t+x))/(a^2))) (Sqrt[π]))/(Sqrt[(a^2)]))

Manipulate[
Plot[Re[Ψtmassless[a,k0,x,t]]/.{a->1,k0->5,m->1},{x,-5,100},PlotRange->{All,{-1,1}},PlotPoints->200,ImageSize->plotsize,PlotStyle->Directive[Blue,Thickness[0.001]],Frame->True,FrameStyle->Directive[Black,20],FrameLabel->{"x","Ψ[x,t]","",""},Axes->False],
{{t,1,"Time"},0,10},LabelStyle->{Black,Large,Italic}]

Manipulate[
Plot[{Re[Ψt[a,k0,x,t,m]]/.{a->1,k0->5,m->1},Re[Ψtmassless[a,k0,x,t]]/.{a->1,k0->5,m->1}},{x,-5,100},PlotRange->{All,{-1,1}},PlotPoints->200,ImageSize->plotsize,PlotStyle->Directive[Blue,Thickness[0.001]],Frame->True,FrameStyle->Directive[Black,20],FrameLabel->{"x","Ψ[x,t]","",""},Axes->False],
{{t,1,"Time"},0,10},LabelStyle->{Black,Large,Italic}]

Old code from last year

(* General form of wavepacket as Fourier transform *)
(∫_(-∞)^∞)g[k]Exp[I (k x -ω[k] t)]\[DifferentialD]k

(* Expand ω[k] around k=k0, set g[k]=1 *)
(∫_(-∞)^∞)Exp[-(((δk)^2)/(2(Δ^2)))+I(x-vg t)δk-I ((σ t)/2)((δk)^2)]\[DifferentialD]δk

<<PhysicalConstants`

Clear[k]
ℏ=1;PlanckConstant/(2 π Joule Second);
c=5;SpeedOfLight Second/Meter;
ωmassive[k_,m_]=(ℏ/2m)(k^2)
ωmassless[k_]=c k
vgmassive[k_,m_]=(∂_k) ωmassive[k,m]
σmassive[k_,m_]=(∂_k) vgmassive[k,m]
vgmassless[k_]=(∂_k) ωmassless[k]
σmassless[k_]=(∂_k) vgmassless[k]

(∫_(-∞)^∞)Exp[-(((δk)^2)/(2(Δ^2)))+I(x-vg t)δk-I ((σ t)/2)((δk)^2)]\[DifferentialD]δk

IntResult=Assuming[Re[(1/(Δ^2))]>Im[t σ],(∫_(-∞)^∞)Exp[-(((δk)^2)/(2(Δ^2)))+I(x-vg t)δk-I ((σ t)/2)((δk)^2)]\[DifferentialD]δk]

Envelopemassive[x_,t_]=IntResult/.{vg->vgmassive[k0,m],σ->σmassive[k,m],ω0->ωmassive[k0,m]}
WPmassive[x_,t_]=Exp[I(k0 x-ω0 t)]×Envelopemassive[x,t]/.ω0->ωmassive[k0,m]

Envelopemassless[x,t]=IntResult/.{vg->vgmassless[k0],σ->σmassless[k],ω0->ωmassless[k0]}
WPmassless[x_,t_]=Exp[I(k0 x-ω0 t)]×Envelopemassless[x,t]/.ω0->ωmassless[k0]

Manipulate[
Plot[Re[WPmassive[x,t]]/.{Δ->1,k0->5,m->1},{x,-5,100},PlotRange->{All,{-3,3}},PlotPoints->200,ImageSize->plotsize,PlotStyle->Directive[Blue,Thickness[0.001]],Frame->True,FrameStyle->Directive[Black,20],FrameLabel->{"x","Ψ[x,t]","",""},Axes->False],
{{t,1,"Time"},0,10},LabelStyle->{Black,Large,Italic}]

Manipulate[
Plot[Re[WPmassless[x,t]]/.{Δ->1,k0->5},{x,-5,100},PlotRange->{All,{-3,3}},PlotPoints->200,ImageSize->plotsize,PlotStyle->Directive[Red,Thickness[0.001]],Frame->True,FrameStyle->Directive[Black,20],FrameLabel->{"x","Ψ[x,t]","",""},Axes->False],
{{t,10,"Time"},0,10},LabelStyle->{Black,Large,Italic}]

Manipulate[
Plot[{Re[WPmassive[x,t]]/.{Δ->1,k0->5,m->1},Re[WPmassless[x,t]]/.{Δ->1,k0->5}},{x,-5,100},PlotRange->{All,{-3,3}},PlotPoints->200,ImageSize->plotsize,PlotStyle->{Directive[Blue,Thickness[0.001]],Directive[Red,Thickness[0.001]]},Frame->True,FrameStyle->Directive[Black,20],FrameLabel->{"x","Ψ[x,t]","",""},Axes->False],
{{t,1,"Time"},0,10},LabelStyle->{Black,Large,Italic}]

Black Body Radiation

(* 1D Planck distribution *)
Clear[h,c,kB]
u1D[ν_]:=((2 h)/c) (ν/(Exp[h ν/(kB T)]-1))

(* Integrate to get Stefan-Boltzmann law *)
(∫_0^∞)u1D[ν]\[DifferentialD]ν

(* With assumptions *)
SBLaw1D[T_]=Assuming[{h>0,kB>0,T>0},(∫_0^∞)u1D[ν]\[DifferentialD]ν]

(* 2D *)
u2D[ν_]:=((8 π h)/(c^2)) ((ν^2)/(Exp[h ν/(kB T)]-1))
SBLaw2D[T_]=Assuming[{h>0,kB>0,T>0},(∫_0^∞)u2D[ν]\[DifferentialD]ν]

(* 3D *)
u3D[ν_]:=((8 π h)/(c^3)) ((ν^3)/(Exp[h ν/(kB T)]-1))
SBLaw3D[T_]=Assuming[{h ϵ Reals,h>0,kB ϵ Reals,kB>0,T ϵ Reals,T>0},(∫_0^∞)u3D[ν]\[DifferentialD]ν]

(* But NB this is the total energy density inside a BB *)
(* To get the power per unit area emitted, extra geometrical considerations are necessary (c/4 factor for 3D) *)

Clear[h,c,kB]
pars={h->PlanckConstant/(Joule Second),c->SpeedOfLight/(Meter/Second),kB->BoltzmannConstant/(Joule/Kelvin)}

u1D[f]/.pars

Plot Planck distribution

Plot[u1D[f]/.T->300/.pars,{f,0,(10^14)},PlotRange->All,Frame->True,FrameStyle->Directive[16,Black],FrameLabel->{"Frequency [Hz]","Energy density [(Hz^(-1))(m^(-n))]","Planck distribution in nD",""},ImageSize->plotsize]

Plot[u2D[f]/.T->300/.pars,{f,0,(10^14)},PlotRange->All,Frame->True,FrameStyle->Directive[16,Black],FrameLabel->{"Frequency [Hz]","Energy density [(Hz^(-1))(m^(-n))]","Planck distribution in nD",""},ImageSize->plotsize]

Plot[u3D[f]/.T->300/.pars,{f,0,(10^14)},PlotRange->All,Frame->True,FrameStyle->Directive[16,Black],FrameLabel->{"Frequency [Hz]","Energy density [(Hz^(-1))(m^(-n))]","Planck distribution in nD",""},ImageSize->plotsize]

Plot[{u3D[f]/(T^4)/.T->500/.pars,
u3D[f]/(T^4)/.T->300/.pars,
u3D[f]/(T^4)/.T->77/.pars},
{f,0,(10^14)},PlotRange->All,Frame->True,FrameStyle->Directive[16,Black],PlotPoints->200,FrameLabel->{"Frequency [Hz]","Energy density [(Hz^(-1))(m^(-n))]","Planck distribution in nD",""},ImageSize->plotsize]

Wien Displacement Law

(* Find the maximum of the Planck distribution *)

(∂_ν)u1D[ν]

(∂_ν)u2D[ν]

(∂_ν)u3D[ν]

Solve[(∂_ν)u1D[ν]==0,ν]
Solve[(∂_ν)u2D[ν]==0,ν]
Solve[(∂_ν)u3D[ν]==0,ν]

NSolve[(∂_ν)u1D[ν]==0,ν]
NSolve[(∂_ν)u2D[ν]==0,ν]
NSolve[(∂_ν)u3D[ν]==0,ν]

u3D[c/λ](∂_λ)(c/λ)

pars

NSolve[(∂_λ)(u3D[c/λ](∂_λ)(c/λ))==0.0,λ]