Using Mathematica in Advanced Undergraduate Physics and Astrophysics Instruction and Research

David Meisel*, Charles Freeman*, and Kirk Anne**

SUNY Geneseo
Geneseo, NY 14454


*Department of Physics and Astronomy
**Computing and Information Technology

Background

According to the American Institute of Physics figures for 2013, the physics department at SUNY-Geneseo is ranked 1st in the nation in terms of the number of Bachelor’s physics degrees produced by an undergraduate department. Because a large fraction of our graduating students are bound for graduate and professional schools, the upper division physics instruction is particularly rigorous and includes a large component of graduate level research activity cleverly disguised as “undergraduate”. Within this context, it has been difficult to support these activities with the usual programming courses and computational packages that are typical of liberal arts colleges. About 5 years ago, we started using Mathematica as the sole computational support in several of our advanced undergraduate courses (statistical thermodynamics, advanced mechanics, and advanced electricity and magnetism) with great success. Since then we have managed to convince the entire department to switch to Mathematica instead of MathCad. This was brought about largely through student pressure on the faculty as many students were already making extensive use of Wolfram/Alpha on their own. This activity has already led to one textbook “Astrophysics Through Computation with Mathematica Suppport”, Koberlein and Meisel (Cambridge, 2013). Mathematica use now explicitly includes four upper division physics courses, advanced laboratories (above freshman level) and research.

Anatomy of our Mathematica Adoption

1999 - The first shoot out. Meisel needed solution to old problem (Watson, Theoretical Astronomy, 1869, Dover reprint 1964) involving hyperbolic orbit determination that was not covered in modern books. Tried Maple®, MathCad®, MatLab®, IDL® and Mathematica. Only Mathematica produced a useful solution in record time and without extraordinary programming effort.

2000 - Meisel began using Mathematica exclusively for personal and professional computations in physics, electrical engineering, astronomy and astrophysics and encouraged others to do likewise. 

2001 - Meisel began to use Mathematica in senior level astrophysics and statistical thermodynamics by providing physics students with unrefined Mathematica notebooks to be used as templates for the solutions of their own problems.

2003-2008 - Started discussions with Anne about extending use of Mathematica within Physics department. 2005 - Meisel retired (he thought...), 2007 - Meisel resumed physics research / teaching. Continued using Mathematica in classes ~20% of physics majors having some exposure to Mathematica. New astronomy faculty person Aaron Steinhauer successfully used some of the astrophysics notebooks in class and student research. Meisel using and developing more notebooks in statistical thermodynamics.

2009 - Savi Iyer and Geneseo student Edward Hansen study analytical solutions for light deflection around spinning balck holes using Mathematica (Iyer and Hansen, 2009). Meisel started working with Brian Koberlein on astrophysics book, “Astrophysics through Computation” for Cambridge. The success in these other courses persuaded Freeman to use Mathematica in the advanced mechanics course for both demos and student problems.  Visit to campus by Wolfram Academic Representative who gave public talk for students and faculty. Wolfram/Alpha Introduced. Student exposure to Mathematica Wolfram products increased dramatically as a result.

2010 - Meisel, Anne and Steinhauer start physics students on N-body calculations as special projects using NDSolve[ ]. Offered special Mathematica workshops for students and faculty. Several Math faculty agree to allow students to submit problem solutions in Mathematica as part of their regular course work.

2012 - Senior physics majors request some instruction in Mathematica be provided prior to the junior level. Anne works with Ed Pogozelski, instructor of the Math Methods (sophomore-level physics course) by attending class, presenting a Mathematica introduction and to supply Mathematica versions of problems previously programmed in MathCad®.  Students beginning to use Mathematica in summer research projects.

2013 - Mathematica used in advanced electricty and magnetism course.   Mathematica exclusively used in Math Methods course, MathCad® license dropped. The book “Astrophysics Through Computation with Mathematica Support” appears. Mathematica used by Anne in numerical analysis course. Course materials for astrophysics, advanced mechanics, statistical mechanics, advanced electricity and magnetism, math methods, and numerical analysis continue to be developed.

2013 - The latest shoot out. In a numerical analysis class, challenge problems will be assigned to students to be solved with Maple®, MatLab®, and Mathematica. It is expected that Mathematica will win widespread student approval.

Main Advantages of Using Mathematica in Undergraduate Classes

1) Cross platform development and notebook exchange is superior to other packages.

2) Mathematica online documentation is superior to other packages. (The trick to quick access that you have to teach students (and some faculty) how to formulate key words).

3) Sophisticated Mathematica prototyping is much faster and more compact than that in the usually available programming languages and systems. Prototyping and run speed is much faster than can be achieved with the usual write, submit, compile, correct cycle such as experienced with say R and students appreciate that.

4) Although texts using specific systems tend to become obsolete quickly, updating of “out-of-date” Mathematica notebooks written in earlier versions usually requires much less effort than similar updates in standard programming languages. These can be used in tutorial sessions whether by instructors or by peer tutors.

5) The ability to document each notebook in great and accurate detail makes it possible to create “templates” from which students (and other faculty) can learn to program similar types of problems or to play “what-if” with the parameter space. This enchances the learning value of the notebooks considerably and makes an extended development effort rewarding. Often in such tempate notebooks we leave out extensive documentation and plot labels and require students to provide such information. This seems to aid their understanding of the subject material.

6) The library depth of Mathematica is truly astounding and increases dramatically with each new version as the structure becomes more and more inclusive while the evolution of other packages is progressively fragmented.

7) Forty years ago computer science gurus promised 5th generation programming languages to lighten one’s computational load. It is our opinion that Mathematica is the only system that approaches that promise. Happy 25th Birthday Mathematica !

8) Currently we are working on using Grid Mathematica to develop a campus wide Mathematica  cluster mainly for undergraduate research. Eventually we hope it will be a gateway for student and faculty access to available supercomputers.

The Problem That Started It All at SUNY-Geneseo (Static)

Starting with Recursion Formulae Obtained from Watson, 1869 Mathematica Produced Analytical Functions.  ArcTan Functions are for Elliptical Orbits, ArcTanh Functions are for Hyperboilc orbits.

uu[τ]=Sum[30*n*(((-1))^n)*(τ^n)/(3+2*(n-1)),{n,1,∞}]

u[τ]=Sum[-30*n*(τ^n)/(3+2*(n-1)),{n,1,∞}]

vv[τ]=Sum[(10+2*n)*(((-1))^((n-1)))*(τ^n)/(1+2*n),{n,0,∞}]

v[τ]=Sum[(10+2*n)*(τ^n)/(1+2*n),{n,0,∞}]

aa=uu[τ]/vv[τ]

a=u[τ]/v[τ]

Producing Expansions Far Beyond Those Shown by Watson

cc=Series[a,{τ,0,10}]/τ

dd=Series[aa,{τ,0,10}]/τ

Operational Formulae Needed for Elliptical/Hyperbolic Orbit Computations

ww[τ]=2*(Sqrt[τ])*vv[τ]

w[τ]=2*(Sqrt[τ])*v[τ]

bb=ww[τ]/(Sqrt[20*Abs[aa]])

b=w[τ]/(Sqrt[20*Abs[a]])

Simplify[bb]

Simplify[b]

Some New Mathematica Technology Samples

Disclaimer: 
Note that these examples are not engineered to be as 
elegant as the usual “DEMO” examples found in the 
Demonstration Project Section of the Wolfram Web Site.

 Instead they are designed to be working templates that can 
 be readily understood by upper division undergraduate science 
  students just learning to make better use of various advanced
  Mathematica features.

Astrophysics Applications :The Book

;

Synchrotron Radiation

Technology Used: Manipulate, Units, Sliders (Run)

Clear["Global`*"]

Initialization using the Units functions, but modified so references to online values are eliminated during the talk.

(*c1=QuantityMagnitude[UnitConvert[Quantity[ "SpeedOfLight"]]]*)

c1=299792458

(*qq=QuantityMagnitude[UnitConvert[Quantity["ElectronCharge"]]]*)

qq=1.60217656500000003606379190082454356534`7.359612335017222*^-19

(*me=QuantityMagnitude[UnitConvert[Quantity["ElectronMass"]]]*)

me=.10938291707429769542081766`7.0617630921243855*^-31

(*hh=QuantityMagnitude[UnitConvert[Quantity["PlanckConstant"]]]*)

hh=6.6260695699999995915614325125943122898228578`7.057827997854494*^-34

(*kk=QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"]]]*)

hh=1.3806488000000000456839657751372779`5.725109871795768*^-23

dPdΩ[θ_,β_]:=((qq^2)/(4 π (c1^3)))(((Sin[θ °])^2)/(((1-β Cos[θ °]))^6));

Linear Acceleration

The pattern displayed is in a plane cut that includes the velocity vector (the straight line).  In this case, the acceleration vector is along the same line and the observer is viewing perpendicular to the plane along y axis.

Manipulate[Show[ParametricPlot[{dPdΩ[tt,β1]Cos[tt °],dPdΩ[tt,β1]Sin[tt °]},{tt,-90,90},PlotRange->All,Axes->False],ParametricPlot[{dPdΩ[tt,β1]Cos[tt °],0},{tt,0,90},PlotRange->All,Axes->False]],{β1,0,(1-(10^(-12)))}];

;

General  Acceleration  Direction (includes Transverse Case)

dPdΩc[θ1_,ϕ1_,inc_,β_]:=(((qq^2)(vdot^2))/(4 π (c1^3))) ((1/(((1-β Cos[θ1 °]))^4))+((2  (Sin[θ1 °]Cos[ϕ1 °]Sin[inc °]+Cos[θ1 °]Cos[inc °])β Cos[inc °])/(((1-β Cos[θ1 °]))^5))-(((1-(β^2)) (((Sin[θ1 °]Cos[ϕ1 °]Sin[inc °]+Cos[θ1 °]Cos[inc °]))^2) )/(((1-β Cos[θ1 °]))^6)));

Here is the relativistic radiation pattern for the case where the velocity and acceleration are at arbitrary angles with respect to each other. The pattern displayed is the section that is in a plane containing the velocity. It is assumed that v is in the z-direction and vdot is in the x-z plane with the k vector pointing in the {x, y, z} = {Sin θ Cos ϕ, Sin θ Sin ϕ, Cos θ } direction.  The acceleration in the x-z plane makes an angle “inc =ip” with respect to the velocity. The full expression is extremely complicated and includes the previous case as inc =0.

The polar angle of the observer is ϕϕ and the ip is the inclination angle between the vdot and v.  Setting inc = 0 gives the collinear situation (above) while inc = 90 is the perpendicular situation with the acceleration in the line of sight.  The default is ip= 90 so that the non-collinear case comes up first. If you want to see the collinear case again set the slider to ip=0.

(*ip=90;*)

DynamicModule[{},{Slider[Dynamic[ϕϕ],{0,90}],"ϕϕ"->Dynamic[ϕϕ]}]

DynamicModule[{},{Slider[Dynamic[ip],{0,90}],"ip"->Dynamic[ip]}]

vdot=1;

pattern=Manipulate[Show[ParametricPlot[{dPdΩc[tt,ϕϕ,ip,β1]Cos[tt °],dPdΩc[tt,ϕϕ,ip,β1]Sin[tt °] Cos[ϕϕ °] },{tt,-180,180},Axes->False,PlotRange->All],ParametricPlot[{dPdΩc[tt,ϕϕ,ip,β1]Sin[tt °] Cos[ϕϕ °],0},{tt,0,90},Axes->False]],{β1,0,(1-(10^(-12)))}];

;

Microwave Galactic Center Direction Analysis (WMAP)

Technology Used: Image Processing, Fourier Transforms, Wavelet Packet Transforms (Static)

WMAP Fluctuations Within 22.5 degrees of the Galactic Center

;

WMAP Power Spectrum for the Region Within 22.5 degrees of the Galactic Center

;

WMAP Wavelet Packet Spectrum for the Region Within 22.5 degrees of the Galactic Center

;

Gravitational N-Body Problem

Technology Used: NDSolve (Symplectic Partitioned Runge-Kutta - SPRK Method), Variational Methods Package, Graphics 3D, ListAnimate,Trekker

A Research (Student Project) Study of Dynamics of the Orion Nebula Trapezium System

Here is the nebula via Hubble Space Telescope (Trapezium = (θ^1) Ori is circled)

;

Here is a simulation of the 3D motion of the components as seen from earth over 1 million years in the past to the present. Red =  (θ^1) Ori A, Green =  (θ^1) Ori B, White =  (θ^1) Ori C, and Blue =  (θ^1)  Ori D. (Slider left goes into past)

;

Here is a simulation of the 3D motion of the components as seen from earth over the present to 1 million years into the future. (Slider right goes into future.)

;

Some Classical Mechanics Applications

Hohmann planetary transfer orbits

Technology Used: Manipulate, Sliders, Graphics Constructions, NDSolve (Run)

C. Freeman

Clear["Global`*"]

r1=1.5;r2=2.75;
T1=r1^(1.5);T2=r2^(1.5);w1=2 Pi/T1;w2=2 Pi/T2;k=4*Pi^2;μ=1.0;

at=(r1+r2)/2.;Tt=at^(1.5);phi=w1*tlaunch;

Opp=Pi-w2 Tt/2.0;ϵ=(r2-r1)/(r2+r1);α=at*(1-ϵ^2);ell=Sqrt[α μ k];topp=If[r1>r2,(Opp+2 Pi)/(w2-w1),(Opp-2 Pi)/(w2-w1)];

tmax=Tt/2+tlaunch;tlaunch=topp;

s=NDSolve[{q'[t]==(ell/(μ (α^2)))(((1+ϵ Cos[q[t]-phi]))^2),q[tlaunch]==phi},q,{t,0,tmax}];

θ[t_]:=q[t]/.Flatten[s]

r[qq_]:=α*(1+ϵ Cos[qq-phi])^(-1)

PolarPlot[r[qq],{qq,0,2*Pi}];

X[t_]:=If[t<tlaunch,r1*Cos[w1*t],r[θ[t]]*Cos[θ[t]]]
Y[t_]:=If[t<tlaunch,r1*Sin[w1*t],r[θ[t]]*Sin[θ[t]]]

g1=PolarPlot[{r1,r2},{gg,0,2 Pi},PlotStyle->{Blue,Red}];

g2[t_]:=ParametricPlot[{{u r1 Cos[w1 t],u r1 Sin[w1 t]},{u r2 Cos[w1 t+Opp],u r2 Sin[w1 t+Opp]}},{u,0,1},PlotStyle->{Purple,Purple}]

g3[t_]:=ParametricPlot[{X[tt],Y[tt]},{tt,-0.1,t},PlotStyle->{Darker[Green],Dashed,Thin}]

Manipulate[Show[ListPlot[{{{X[t],Y[t]}},{{r1 Cos[w1 t],r1 Sin[w1 t]}},{{r2 Cos[w2 t+0*Opp],r2 Sin[w2 t+0*Opp]}}},PlotRange->{{-3,3},{-3,3}},AspectRatio->Automatic,PlotMarkers->{●},PlotStyle->{Darker[Green],Blue,Red}],g3[t],g2[t],g1],{t,0,tmax,tmax/1000.0}]

Tracking a Foucalt Pendulum

Technology Used: Manipulate, Sliders, Graphics Constructions, NDSolve (Run)

C. Freeman

Clear["Global`*"]

r0=0.9;w=0.3;α=10;

tmax= π/w

sol=NDSolve[{x''[t]==0(w^2) x[t]+2 w y'[t]-(α^2) x[t],y''[t]==-2 w x'[t]+0(w^2) y[t]-(α^2) y[t],x[0]==r0,x'[0]==0,y[0]==0,y'[0]==-w r0},{x,y},{t,0,tmax}];

X[t_]:=x[t]/.Flatten[sol]
Y[t_]:=y[t]/.Flatten[sol]

Rot[t_]:={{Cos[w t],-Sin[w t]},{Sin[w t], Cos[w t]}}

Xirf[t_]:=First[Rot[t].{X[t],Y[t]}]
Yirf[t_]:=Last[Rot[t].{X[t],Y[t]}]

Xt[t_,tf_]:=First[Rot[tf].{X[t],Y[t]}]
Yt[t_,tf_]:=Last[Rot[tf].{X[t],Y[t]}]

gcirc=Graphics[Circle[{0,0},1]];

gline[t_]:=Graphics[{Gray,Line[{{-Cos[w t],-Sin[w t]},{Cos[w t],Sin[w t]}}],Line[{{Sin[w t],-Cos [w t]},{-Sin[w t],Cos[w t]}}]}]

gt[tf_]:=ParametricPlot[{Xt[t,tf],Yt[t,tf]},{t,0,tf}]

girf[tf_]:=ParametricPlot[{Xirf[t],Yirf[t]},{t,0,tf},PlotStyle->{Black,Dashed}]

gg1[tf_]:=ParametricPlot[{X[t],Y[t]},{t,0,tf}];

gpuck[tf_]:=Graphics[Disk[{X[tf],Y[tf]},0.03]]
gpuck2[tf_]:=Graphics[Disk[{Xirf[tf],Yirf[tf]},0.03]]

gred[t_]:=Graphics[{Red,Disk[{Cos[w t],Sin[w t]},0.05]}]
gblue[t_]:=Graphics[{Blue,Disk[{Cos[w t+Pi/2],Sin[w t+Pi/2]},0.05]}]
ggreen[t_]:=Graphics[{Green,Disk[{Cos[w t+Pi],Sin[w t+Pi]},0.05]}]
gpurple[t_]:=Graphics[{Purple,Disk[{Cos[w t+3Pi/2],Sin[w t+3Pi/2]},0.05]}]

ani1[tf_]:=Show[gg1[tf],gcirc,gred[0],gblue[0],ggreen[0],gpurple[0],gpuck[tf],gline[0],PlotRange->{{-1.2,1.2},{-1.2,1.2}},Axes->None]

ani2[tf_]:=Show[gt[tf],girf[tf],gcirc,gred[tf],gblue[tf],ggreen[tf],gpurple[tf],gpuck2[tf],gline[tf],PlotRange->{{-1.2,1.2},{-1.2,1.2}},Axes->None]

Manipulate[ani1[tf],{tf,0.01,tmax}]

Manipulate[GraphicsGrid[{{ani1[tf],ani2[tf]}}],{tf,0.01,tmax}]

Precessional Motion

Technology Used: Manipulate, Sliders, Graphics Constructions, NDSolve (Run)

C. Freeman

Clear["Global`*"]

I1=2;I3=4;wmag=1.5;β=30 Degree;w3=wmag Cos[β];Ω=(I3-I1)*w3/I1;

Lmag=Norm[{I1*wmag*Sin[β],0,I3*w3}];ϕdot=Lmag/I1;

α=ArcCos[(I1*(wmag*Sin[β])^2+w3*I3*w3)/(wmag*Lmag)];

ww[t_]:=wmag*{Sin[α]Cos[ϕdot t],Sin[α]Sin[ϕdot t],Cos[α]};

x3[t_]:=wmag*{Sin[β+α]Cos[ϕdot t],Sin[β+α]Sin[ϕdot t],Cos[β+α]}

wg[t_]:=Graphics3D[Arrow[{{0,0,0},ww[t]}],PlotRange->{{-1.5,1.5},{-1.5,1.5},{-2,2}}]

x3g[t_]:=Graphics3D[Arrow[{{0,0,0},x3[t]}]]

Lg=Graphics3D[Arrow[{{0,0,0},{0,0,wmag}}]];

cone=Graphics3D[{Opacity[0.8],Cone[{wmag*{0,0,Cos[α]},{0,0,0}},wmag*Sin[α]]}];

cone2[t_]:=Graphics3D[{Opacity[0.5],Cone[{wmag*{Sin[β+α]Cos[ϕdot t],Sin[β+α]Sin[ϕdot t],Cos[β+α]},{0,0,0}},wmag*Tan[β]]}]

Manipulate[Show[wg[t],x3g[t],Lg,cone,cone2[t]],{t,0,2π/ϕdot}];

;

Motion of a Tossed Tennis Racket

Technology Used: Manipulate, Sliders, Graphics Constructions, NDSolve (Run)

C. Freeman

Clear["Global`*"]

I1=1;I2=1.594;I3=2.559;tmax=100;mass=1.25;

a=Sqrt[(I2+I3-I1)/mass];b=Sqrt[(I1+I3-I2)/mass];c=Sqrt[(I1+I2-I3)/mass];

wcomp={0,1,0.02};

wvec=wcomp Norm[wcomp];

sol=NDSolve[{I1 w1'[t]==(I2-I3)w2[t]w3[t],I2 w2'[t]==(I3-I1)w3[t]w1[t],I3 w3'[t]==(I1-I2)w1[t]w2[t],w1[0]==wvec[[1]],w2[0]==wvec[[2]],w3[0]==wvec[[3]]},{w1,w2,w3},{t,0,tmax}];

W1[t_]:=w1[t]/.Flatten[sol]
W2[t_]:=w2[t]/.Flatten[sol]
W3[t_]:=w3[t]/.Flatten[sol]

Wbody[t_]:={W1[t],W2[t],W3[t]}

Lbody[t_]:={{I1,0,0},{0,I2,0},{0,0,I3}}.Wbody[t]/Norm[{{I1,0,0},{0,I2,0},{0,0,I3}}.Wbody[t]]

gWbody[t_]:=Arrow[{{0,0,0},Wbody[t]}]

e1={1,0,0};e2={0,1,0};e3={0,0,1};

body=GraphicsGroup[{Arrow[{{0,0,0},e1}],Arrow[{{0,0,0},e2}],Arrow[{{0,0,0},e3}]}];

cube[t_]:=Cuboid[{-a/2,-b/2,-c/2},{a/2,b/2,c/2}]

Manipulate[Graphics3D[{Opacity[0.9],cube[t],Thick,Red,Arrow[{{0,0,0},Wbody[t]}],Thick,Blue,Arrow[{{0,0,0},Lbody[t]}],Thin,Black,Dashed,body},PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1.1,1.1}}],{t,0,tmax}];

;

In the above simulation we are in the body-fixed frame.  The blue arrow is the angular momentum vector, and the red arrow is the angular velocity vector.

sol2=NDSolve[{W1[t]==ϕ'[t]Sin[θ[t]]Sin[ψ[t]]+θ'[t]Cos[ψ[t]],W2[t]==ϕ'[t]Sin[θ[t]]Cos[ψ[t]]-θ'[t]Sin[ψ[t]],W3[t]==ϕ'[t]Cos[θ[t]]+ψ'[t],θ[0]==0.01,ϕ[0]==0.01,ψ[0]==0.01},{ϕ,θ,ψ},{t,0,tmax}];

ϕout[t_]:=ϕ[t]/.Flatten[sol2]
θout[t_]:=θ[t]/.Flatten[sol2]
ψout[t_]:=ψ[t]/.Flatten[sol2]

λ1[t_]:={{Cos[ϕout[t]],Sin[ϕout[t]],0},{-Sin[ϕout[t]],Cos[ϕout[t]],0},{0,0,1}}
λ2[t_]:={{1,0,0},{0,Cos[θout[t]],Sin[θout[t]]},{0,-Sin[θout[t]],Cos[θout[t]]}}
λ3[t_]:={{Cos[ψout[t]],Sin[ψout[t]],0},{-Sin[ψout[t]],Cos[ψout[t]],0},{0,0,1}}

λ[t_]:=λ3[t].λ2[t].λ1[t]

E1[t_]:=Transpose[λ[t]].e1;E2[t_]:=Transpose[λ[t]].e2;E3[t_]:=Transpose[λ[t]].e3;

bodyaxes[t_]:=Show[Graphics3D[{Dashed,Arrow[{{0,0,0},E1[t]}]}],Graphics3D[{Dashed,Arrow[{{0,0,0},E2[t]}]}],Graphics3D[{Dashed,Arrow[{{0,0,0},E3[t]}]}]];

v1[t_]:=Transpose[λ[t]].{a/2,b/2,c/2}
v2[t_]:=Transpose[λ[t]].{a/2,-b/2,c/2}
v3[t_]:=Transpose[λ[t]].{a/2,-b/2,-c/2}
v4[t_]:=Transpose[λ[t]].{a/2,b/2,-c/2}
opa=0.9;

cu[t_]:=GraphicsGroup[{Opacity[opa],White,Polygon[{v1[t],v2[t],v3[t],v4[t]}],Red,Polygon[{v1[t],v2[t],-v4[t],-v3[t]}],Blue,Polygon[{v2[t],v3[t],-v1[t],-v4[t]}],Yellow,Polygon[{-v1[t],-v2[t],-v3[t],-v4[t]}],Orange,Polygon[{-v1[t],-v2[t],v4[t],v3[t]}],Green,Polygon[{-v2[t],-v3[t],v1[t],v4[t]}]}]

Manipulate[Graphics3D[{cu[t],Thick,Red,Arrow[{{0,0,0},Transpose[λ[t]].Wbody[t]}],Thick,Blue,Arrow[{{0,0,0},Transpose[λ[t]].Lbody[t]}],Thin,Black,Dashed,Arrow[{{0,0,0},E1[t]}],Arrow[{{0,0,0},E2[t]}],Arrow[{{0,0,0},E3[t]}]},PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1.1,1.1}}],{t,0,tmax,0.01}];

;

In this simulation we are in the inertial coordinate system.  Again, the blue arrow is the angular momentum vector (constant in the inertial frame) and the red arrow is the angular velocity vector.  Notice how the block “flips” when the angular velocity vector is closely aligned with the 2-axis (intermediate moment of inertia).

Simple Projectile Motion 155mm Howitzer Firing an HE M107 Shell

Technology Used:  Manipulate, Sliders, Vector Analysis, Graphics3D, ParametricPlot3D. (Run)

K. Anne

Clear["Global`*"]
Ignore[Set::shape];
Ignore[CoordinateTransformData::bdpt];
Ignore[Solve::ratnz];
Ignore[Solve::fulldim];

Manipulate[
ShellMass=43;               (* HE M107 Shell in kg *)
TargetRange =10000; (* Maximum Range 14,600 m *)
(* 
θ=45Degree;                  (* Elevation from -2° to 63° *)
ϕ=23 Degree;                  (* Traverse from -25° to 25° *)
*)
MuzzleVelocity=564;  (* Muzzle velocity of 155 mm Howitzer in m/s *)
BarrelLength=3.564;
Acceleration=((MuzzleVelocity^2)/(2 BarrelLength));
(Subscript[T, b])=BarrelLength/MuzzleVelocity;
{(Subscript[x, 0]),(Subscript[y, 0]),(Subscript[z, 0])}={0,0,zinit};
{(Subscript[a, x]),(Subscript[a, y]),(Subscript[a, z])}={0,0,-9.8};
{(Subscript[Vx, 0]),(Subscript[Vy, 0]),(Subscript[Vz, 0])}=CoordinateTransformData["Spherical"->"Cartesian","Mapping",{MuzzleVelocity,90 Degree - θ,ϕ}];
{(Subscript[Ax, 0]),(Subscript[Ay, 0]),(Subscript[Az, 0])}=CoordinateTransformData["Spherical"->"Cartesian","Mapping",{Acceleration,90 Degree - θ,ϕ}];
Ax[t_]:=If[t<(Subscript[T, b]),(Subscript[Ax, 0]),(Subscript[a, x])];
Ay[t_]:=If[t<(Subscript[T, b]),(Subscript[Ay, 0]),(Subscript[a, y])];
Az[t_]:=If[t<(Subscript[T, b]),(Subscript[Az, 0]),(Subscript[a, z])];
Vx[t_]:=(Subscript[Vx, 0])+Ax[t]*t;
Vy[t_]:=(Subscript[Vy, 0])+Ay[t]*t;
Vz[t_]:=(Subscript[Vz, 0])+Az[t]*t;
x[t_]:=(Subscript[x, 0])+Vx[t]*t+Ax[t]*(t^2);
y[t_]:=(Subscript[y, 0])+Vy[t]*t+Ay[t]*(t^2);
z[t_]:=(Subscript[z, 0])+Vz[t]*t+Az[t]*(t^2);
solutions=Solve[z[t]==ztarget,t];
timpact=Max[Select[t /. solutions,#>0 &]];
s=Table[Graphics3D[Cylinder[{{x[timpact],y[timpact],z[timpact]+20},{x[timpact],y[timpact],z[timpact]-20}},i]],{i,0,2000,500}];
p=ParametricPlot3D[{x[t],y[t],z[t]},{t,0,100},PlotRange->{{-1000,20000},{-8000,8000},{-10000,10000}},AspectRatio->1];
zp=Graphics3D[Polygon[{{-1000,-8000,ztarget},{-1000,8000,ztarget},{20000,8000,ztarget},{20000,-8000,ztarget},{-1000,-8000,ztarget}}]];
Show[{p,zp,s}],
{{θ,30 Degree, "Elevation"},-2  Degree,63 Degree,1 Degree},
{{ϕ,0 Degree, "Rotation"},-25 Degree,25 Degree,1 Degree},
{{zinit,0,"Altitude"},-10000,10000},
{{ztarget,-10000,"Target Altitude"},-10000,zinit}
,ControlPlacement->Left]

New Insights to Old Problems (2 Examples)

Monte Carlo Solution of the Birthday Coincidence Problem

Technology Used:  RandomReal, Histogram, ListPlot (Static)

K. Anne

How many matches of two or more persons with birthdays on the same date in different class sizes ?

;

What is the probability of  two or more persons with birthdays on the same date in different class sizes ?

;

This graph shows an interesting phenomenon that is not evident at lower values of sample size. Notice that the maximum probability of pairs occurs at size ~ 360 (365 ?), triplets at twice that, ~ 720 (730 ?), and quadruplets at three times that,  ~ 1100 (1095 ?). The maximum probabilities are a factor of two lower each time, 0.18, 0.09, 0.045... It is likely such behavior would only be discovered so easily by Monte Carlo techniques.

Analysis of Ultrasonic Two Circular-Aperture Source Diffraction

Technology Used:  NonlinearModelFit, LinearModelFit, Fourier Analysis, KaiserWindow

C.Freeman, J. Shortino, D. Meisel

Clear["Global`*"]

Original Data

data={{2,5.5696`},{3,5.1076`},{4,5.5696`},{5,7.29`},{6,7.6176`},{7,6.0516`},{8,6.4516`},{9,5.76`},{10,4},{11,3.61`},{12,3.0976`},{13,1.8496`},{14,0.9216`},{15,0.49`},{16,0.25`},{17,0.4356`},{18,2.0164`},{19,2.9584`},{20,4},{21,5.6644`},{22,7.29`},{23,9.2416`},{24,8.8804`},{25,9.1204`},{26,10.89`},{27,10.24`},{28,9},{29,8.7616`},{30,8.41`},{31,6.8644`},{32,3.61`},{33,1.7956`},{34,0.6084`},{35,0.36`},{36,1.3456`},{37,3.0976`},{38,5.3824`},{39,7.29`},{40,10.3684`},{41,13.8384`},{42,16.6464`},{43,19.36`},{44,19.7136`},{45,19.0096`},{46,17.64`},{47,14.44`},{48,9.7344`},{49,6.1504`},{50,3.24`},{51,1.44`},{52,0.4096`},{53,0.5776`},{54,2.3104`},{55,5.3824`},{56,8.2944`},{57,11.8336`},{58,13.8384`},{59,17.9776`},{60,20.7936`},{61,19.0096`},{62,17.3056`},{63,14.1376`},{64,11.56`},{65,8.2944`},{66,4.84`},{67,2.56`},{68,0.8464`},{69,0.1936`},{70,1.2544`},{71,3.5344`},{72,5.76`},{73,9.4864`},{74,12.3904`},{75,14.44`},{76,15.6816`},{77,15.0544`},{78,15.3664`},{79,14.44`},{80,13.5424`},{81,10.4976`},{82,7.6176`},{83,4},{84,2.0736`},{85,0.64`},{86,0.16`},{87,0.4624`},{88,1.7424`},{89,3.5344`},{90,5.3824`},{91,6.3504`},{92,6.76`},{93,8.5264`},{94,10.24`},{95,9.7344`},{96,9.7344`},{97,8.7616`},{98,7.1824`},{99,5.1984`},{100,3.8416`},{101,2.9584`},{102,1.96`},{103,0.8464`},{104,0.3136`},{105,0.0784`},{106,0.2704`},{107,0.7744`},{108,0.8464`},{109,1.2544`},{110,2.8224`},{111,4},{112,4.84`},{113,5.5696`},{114,5.0176`},{115,4.6656`},{116,5.5696`},{117,5.5696`},{118,5.3824`},{119,4.3264`},{120,3.0976`}};

Previous numerical experiments have shown that precision values of the following two parameters are needed for good convergence of NonlinearModelFit[ ]. The value for L was obtained by direct measurement and the precise ultrasonic wavelength was calculated from a frequency meter reading and the speed of sound in air, itself calculated from the ambient air temperature.

L=113.0;λ=1.05512;

Traditional Analysis

Fit of Linearized Diffraction Equation with Allowance for a Possible Linear Trend in the Data due to Non-perpendicularity of Measuring Rail.

fit=NonlinearModelFit[data,Imax*(((Sin[((x-x0)/(Sqrt[(((x-x0))^2)+(L^2)]))*((π*a)/λ)]*Cos[((x-x0)/(Sqrt[(((x-x0))^2)+(L^2)]))*((π*d)/λ)]/(((x-x0)/(Sqrt[(((x-x0))^2)+(L^2)]))*((π*a)/λ))))^2)*((L^2)/((L^2)+(((x-x0))^2)))+I0+b x,{{a,1.5},{d,8.0},{x0,64.945},{I0,0},{b,0},{Imax,25.0}},x];

Plot of Fit and Original Points

Show[Plot[fit[x],{x,0,120},PlotStyle->{Red}],ListPlot[data]]

Goodness of Fit

(Sqrt[fit["RSquared"]])

Parameter and Error List

fit["ParameterTable"]

fitcoeff=fit["BestFitParameters"]

Parameter and Error List

(Notice that the parameters for the linear trend, I0 and b, are above their errors so that the inclusion of those terms must be considered satisfactory and that the detector track indeed was tilted to the plane of the two sources.   With the trend removed, the residuals are now well-balanced (= “flat”), but still the residuals still show signs of remaining quasi-periodicity. )

resid=fit["FitResiduals"];

residplot=ListPlot[resid,Joined->True,ImageSize->Scaled[0.8]]

Classical Fourier Analysis

Because the residuals do not look random, examination of power spectra of the data and residuals may be useful to clarify certain issues.

Power Spectrum of Raw Points with the Trend Removed

cst=fitcoeff[[4,2]];slope=fitcoeff[[5,2]];

data1=Table[data[[ii,2]]-cst-slope data[[ii,1]],{ii,1,119}];

avdata1=Mean[data1];

Mean Subtracted to Reduce Zero-Order Frequency Component

data0=Table[data1[[ii]]-avdata1,{ii,1,119}];

psdata1=((Abs[Fourier[data0]])^2);

Plotting Log 10 of the Power Spectrum for Clarity

rawplot=ListPlot[Log[10,psdata1],Joined->True,PlotRange->{{0,61},{-4,3}},PlotStyle->{Dashed,Blue},AxesLabel->{"Frequency Index","Log 10 Power"},ImageSize->Scaled[0.9]]

First one needs to recognize that the actual frequency is given in terms of the x-axis values (=frequency index) as
		freq((cm^(-1)) ) = 0.5 (frequency index-1)

Secondly, in this raw log 10 spectrum we see that the zero frequency component is zero (log = -∞) as it should be since the trend was removed. The dominant frequency is at 4 (cm^(-1)) (the index position is 8-1=7, equal to the number of narrow peaks visible in the data plot). The next most powerful frequency is 0.5 (cm^(-1))(the index position is 2-1=1, equal to the number of broad peaks visible). But there are other peaks visible, the most intense ones being at the first harmonics of these two fundamental frequencies. These harmonics are signatures of sample aliasing. The remaining part of the spectrum oscillates around a descending linear trend line in the log plot. This indicates a power law type of noise.

Extended Fourier Analysis

Exploring a Windowed Fourier Transform Power Spectrum

As noted above there are certain features of the raw Fourier spectrum that may be artifacts of the discrete sampling used to obtain the original data. These are considered distortions and to minimize these much effort has gone into designing “window” functions that multiply the data in the “time” domain to remedy such artifacts in the frequency domain.  One personal favorite is the Kaiser window based on Bessel functions. Since Bessel functions are often used to design realistic electronic filters, Kaiser windows simulate the action of “realizable” filter banks.

Create Three Kaiser Window Functions

data2=Table[Table[KaiserWindow[ii/120.,3.+( jj-1)]data0[[ii]],{ii,1,119}],{jj,1,3}];

Find the Power Spectra of the Windowed FT

psdata2=((Abs[Fourier[data2]])^2);

smoothplot=ListPlot[Log[10,psdata2],Joined->True,PlotRange->{{0,61},{-4,3}},AxesLabel->{"Frequency Index","Log 10 Power"},PlotStyle->{Blue,Red,Green},ImageSize->Scaled[0.9],PlotLegends->Placed[LineLegend[{{Dashed,Blue},Blue,Red,Green},{"No filter","α=3","α=4","α=5"},LegendLayout->"Row"],{0.55,0.98} ]];

Show[rawplot,smoothplot]

The action of the filter with α =4 or 5 is practically identical in making the fundamental peaks “clean” and the rest of the spectrum smoother. The noise section is still present but of reduced amplitude and with some obvious phase shifts.  The origin of the noise is still not certain but since it is “spatial” noise, it most likely arises from self-reflections within the measurement area itself or perhaps interference from other experimental setups in the room at the same time.

Frequency Analysis of the Residuals

If interference from reflections and other experimental setups in the same room, we examine the spectral properties of the residuals. Here again is a plot of the residuals with the linear trend removed.

residplot

Power Spectrum of the Residuals

pspect0=((Abs[Fourier[resid]])^2);

As suspected, the residuals from the non-linear fit produce the power law spectrum at all frequencies.

ListPlot[Log[10,pspect0],Joined->True,PlotRange->{{0,61},{-3,3}},AxesLabel->{"Frequency Index","Log 10 Power"},ImageSize->Scaled[0.8]]

Converting the frequencies to octaves and the power spectra to dB of intensity we find the slope of the noise in db per octave. Since the data are detrended the residuals have zero power at zero frequency so we omit the intercept and that gives the highest fit correlation.

resdata=Table[{Log[8,0.5(ii-1)],10 Log[10,(Sqrt[pspect0])][[ii]]},{ii,2,61}];

pnoise=LinearModelFit[resdata,x,x,IncludeConstantBasis->False]

pnoise["ParameterTable"]

(Sqrt[pnoise["RSquared"]])

Since the slope is nearly -3 dB/octave, this particular noise is called pink noise because it increases toward lower frequencies as 1/f.  (See Wikipedia http://en.wikipedia.org/wiki/Colors_of_noise)

Autocorrelation of the Residuals

The InverseFourier of the power spectrum produces the autocorrelation function. The autocorrelation is a function of the “lag” distance (here in cm) and is a measure of the interaction distances that produces the pink noise.

auto=((Abs[InverseFourier[pspect0]])^2);

ListPlot[auto,Joined->True,PlotRange->{{0,61},{0,78}},AxesLabel->{"Lag in cm","Autocorrelation"},ImageSize->Scaled[0.8]]

The autocorrelation function away from zero lag decreases until lags of ~10 cm are reached. Thus 10 cm is the approximate correlation length and that includes the distance between the two source centers and the outside edge to edge distance across the two sources. The next peak  is 20 cm which is the distance of travel going from one source to the other and back again. This diagram seems sufficient to indicate that the noise orginates in close range interactions between the two sources. Such interactions are called near-field diffraction effects. They cannot be described by the standard diffraction formula as it represents far-field diffraction, not near field.

Lessons Learned from Our Experiences

1) Inertia among faculty against switching programming languages is fierce with experimentalists often more biased than theoreticians.

2) Adoption of a new, uniform programming system is perceived as requiring considerable “personal” retooling on the part of instructors/researchers to keep themselves ahead of their students.

3) Perceived difficulty and expense of converting legacy code (usually acquired during graduate school and often used in subsequent instructional situations), is another excuse given for not adopting a superior programming system.

4) Today's undergraduate students have no such prejudices and are much more willing to switch programming systems if they see the move is to their advantage.

5) Student pressure for the adoption of a superior programming system is often able to overcome items (1-3) even with availablity of formal programming courses. (Availability of Wolfram/Alpha and a campus visit by the Wolfram college representative were believed to be instrumental in forming student opinion our case).

6) Support of student/faculty peer efforts to adapt to programming system change requires constant attention and patience on the part of those with more experience.  

7) Version updating of out-of-date notebooks from books previously published is an excellent way to develop Mathematica agility among faculty and students. This is consistent with our template approach to Mathematica instruction.

8) Using Mathematica enables undergraduate students to embark and achieve considerable progress (before they graduate) on serious computational research projects such as light propagation around a spinning black hole, the N-body problem under a galactic potential, the Boltzmann transport equation treatment of fluids, and relativistic cosmology.