Mathematica Introduction : 2

- Reminders from last time
- Some general tips
- Nice graphics demonstrations for wave interference patterns
- Pseudo data generation
- Data fitting
- Wavepacket evolution

(* 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]

(* Fill in your local file path here *)
fullPath="/Users/peterleek/Desktop/Mathematica Teaching/";
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

Interference

Single slit diffraction

(* Define function for intensity *)
k:=2π/λ
sinθ:=(y/(Sqrt[(R^2)+(y^2)]))
β[λ_,d_,y_,R_]:=((k d)/2)sinθ
Intensity[λ_,d_,y_,R_]= (((((Sin[β[λ,d,y,R]])/(β[λ,d,y,R]))))^2);

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

Plot[Intensity[λ nm,d μm,y,R]/.{λ->500,R->1,d->10},{y,-0.2,0.2}]

(* Plot intensity as function of screen position *)
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,
ImageSize->plotsize,
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}]

Double slit diffraction

α:=((k a)/2)sinθ
IntensityDS[λ_,d_,y_,R_,a_]= (((((Sin[β[λ,d,y,R]])/(β[λ,d,y,R]))))^2)(((Cos[α]))^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,
ImageSize->plotsize,
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}]

(* 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,
ImageSize->plotsize,
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}]

Diffraction gratings

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

Manipulate[Plot[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,
ImageSize->plotsize,
Frame->True,FrameStyle->Directive[Black,20],
FrameLabel->{"y [m]","I/I0","",""},
AspectRatio->0.2],
{{R,0.5,"Distance [m]"},0.1,5},
{{d,0.5,"Slit width [μm]"},0.1,1},
{{λ,450,"Wavelength [nm]"},400,1000},
{{a,20,"Slit separation [μm]"},5,25},
{{Num,2,"Number of slits"},2,10},
LabelStyle->{Black,Italic,Large}]

A pseudo-experiment

(* Generate random complex numbers (i.e. pairs of random numbers) *)
(* The range for the Im part is 0 to 1, so that it scales to a probability *)
random2D=Table[RandomComplex[{-0.1,0.1+I}],{i,20}]

Clear[R,d,λ,a]
params={R->1,d->1 μm,λ->100 nm,a->2μm};

(* Register a photon 'click' at position Re[#] if IntensityDS[Re[#]]>Im[#] *)
(* Otherwise record a non-event as some number way outside the range (-123) *)
counts=If[(IntensityDS[λ,d,Re[#],R,a]/.params)>Im[#],Re[#],-123]&/@random2D

(* Count up the photon clicks in binned regions of screen position *)
(* Using function 'BinCounts' *)
binwidth=0.0005;
binleft=-0.1;
binright=0.1;
binnedcounts={Table[x,{x,binleft+binwidth/2,binright,binwidth}],BinCounts[counts,{binleft,binright,binwidth}]}//Transpose

(* Now put this all together in a 'pseudo' experiment *)
Clear[R,d,λ,a]
params={R->1,d->1 μm,λ->100 nm,a->2μm};

Manipulate[
random2D=Table[RandomComplex[{-0.1,0.1+I}],{i,2M}];

counts=If[(IntensityDS[λ,d,Re[#],R,a]/.params)>Im[#],Re[#],-123]&/@random2D;
binnedcounts={Table[x,{x,binleft+binwidth/2,binright,binwidth}],BinCounts[counts,{binleft,binright,binwidth}]}//Transpose;
dataplot=ListPlot[binnedcounts,
PlotStyle->{Black,PointSize[0.001]},
PlotRange->All,
Frame->True,
Filling->Axis,
FillingStyle->Directive[Thickness[0.002],Darker[Blue]],
FrameLabel->{"Position on screen [m]","Photon count","",""},
FrameStyle->Directive[20,Black],
Axes->False,
ImageSize->plotsize],{{M,10,"Trials"},10,5000},LabelStyle->{Black,Italic,Large}]

Fit the pseudo data

Clear[x,d,a,R,λ]
fitparams=FindFit[binnedcounts,Nmax IntensityDS[λfit,d,x,R,a]/.params,{{Nmax,30},{λfit,100 nm}},x]

(* Note the use of replacement rules here *)
fitplot=Plot[ Nmax IntensityDS[λfit,d,x,R,a]/.params/.fitparams,{x,-0.1,0.1},PlotStyle->Directive[Red,Thickness[0.005]]];

Show[dataplot,fitplot,Frame->True,FrameLabel->{"Position on screen [m]","Photon count","Photon count data fitted to wave interference",""},FrameStyle->Directive[20,Black]]

For full analysis of data fitting, use package NonLinearRegress.  Here you can for example find confidence intervals for the fit parameters.

Wavepacket Motion (problem set 4)

<<PhysicalConstants`

(∂_x) (x^4)

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]

(* Mathematica can solve this integral, but only if it makes certain assumptions *)
(* Hence it gives the answer inside an 'If' statement *)
(* In general, Mathematica does not assume that variables are purely real, unless you tell it explicitly that this is the case *)
(∫_(-∞)^∞)Exp[-(((δk)^2)/(2(Δ^2)))+I(x-vg t)δk-I ((σ t)/2)((δk)^2)]\[DifferentialD]δk

(* We can also use the 'Assuming' function to do the integral with certain assumptions *)
(* For this particular case, we don't need to worry because all our parameters are real *)
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}]