Visual Python in a Computational Physics Course

1
Visual Python in a Computational Physics Course

• Richard P. Olenick
• David Case
• Eric Pepin
• William Spearman

2
Outline

• Physical Modeling Demos with VPython
• How Does a Slider Differ from a Fastball?
• Boids
• Stability of L4 and L5 Lagrange Points
• How Lightning Forms
• Background and Mathematical Tools
• Euler Methods, Runge-Kutta, Poincare Maps
• Advantages of VPython
• Demos of Students Semester Projects
• 3D Dynamic Lattice
• Gravitational Shaping by Galactic Encounters
• Sand Avalanches
• Percolation Growth of Galaxies
• Non-relativistic Electron Scattering

3
Elements of the Computational Physics Course

• Emphasizes the modeling of physical systems.
• Lectures are kept to a minimum and students work
  on programs during classes.
• Instructor provides in-class debugging.
• Students have different initial conditions to run
  and their results are shown to the entire class
  and discussed.

4
Magnus Effect and Baseball Pitches

• Application of the Magnus effect and empirical
  viscous drag
• Prospective
• Fastballs, curveballs, sliders, and screwballs
  thrown by both right-handed and left-handed
  pitchers can be simulated and viewed from the
  prospective of the pitcher, catcher, or umpire.

5
(No Transcript)
6
def go() global omega, w, ball, g, dt, B,
phi, rot, speed initial velocity of pitch
for throwing angle ball.vvector(-smath.sqrt(
1-math.cos(theta)math.cos(theta)-math.cos(beta)m
ath.cos(beta)),smath.cos(theta),smath.cos(beta))
while( ball.xgt-.5 and ball.ygt0) while ball
in the air and not passed home plate
c.interact() rate(20speed)
ball.rotate(angleomegadt,axisrot,originball.po
s)animates the rotation of the ball
ball.v_mid ball.v (ball.v
(-F(ball.v.mag)ball.v.magball.v
Bcross(w,ball.v)-g)dt/2.)dt/2 calculates new
pos and vel ball.v ball.v
(-F(ball.v_mid.mag)ball.v_mid.magball.v_mid
Bcross(w,ball.v_mid)-g)dt ball.pos
ball.pos ball.vdt ball.path.append(pos
ball.pos) marks path of pitch if a
strike, make ball red if ((ball.pos.xlt0
and ball.pos.xgt-.216)and(ball.pos.ylt1.24 and
ball.pos.ygt.52)and(ball.pos.zlt.216 and
ball.pos.zgt-.216))
ball2.colorcolor.red
7
Boids

• Models the distributive behavior that is observed
  in nature among groups such as swarms, flocks,
  and herds with an application to ferromagnetism.
• Speed, cohesion, separation, alignment and goal
  seeking can be controlled in the simulation. The
  average separation is measured as the boids
  flock.

8
for i in range(numboids) v1
rule_1(i)D1 velocity correction based on
rule 1 (cohesion) v2 rule_2(i)D2
velocity correction based on rule 2
(separation) v3 rule_3(i)D3
velocity correction based on rule 3 (alignment)
v4 tend_to_place(i)D4 velocity
correction based on common goal point
boidsi.velocity Speedness(boidsi.velocityv
1v2v3v4)/abs(boidsi.velocityv1v2v3v4)
boidsi.pos (boidsi.velocity)dt
directedi.pos boidsi.pos
directedi.axis (0.5boidsi.velocity/abs(boid
si.velocity)) total boidsi.pos
9
Stability of Lagrange Points

• The stability of particle orbits near the L4 and
  L5 Lagrange points is simulated.
• The effects of mass ratio, ?, of the two major
  bodies as well as the initial offset positions,
  ?x and ?y can be studied.
• Simulations use the Runge-Kutta-Fehlberg
  algorithm.

10
if graphics 1 draws the scene
objectposvector(x,y,0) use angular
velocity to find how much the objects are
rotated objectposrotate(objectpos,
angleomegat, axis0,0,1)
body2posvector(1-mu, 0,0) draws the
scene   body2posrotate(body2pos,
angleomegat,axis0,0,1)
body1posvector(-mu, 0,0)
body1posrotate(body1pos, angleomegat,axis0,0,
1) object.posobjectpos
body2.posbody2pos body1.posbody1pos
object.trail.append(posobject.pos)
11
How Lightning FormsFractal Growth Processes

• Diffusion-limited growth is simulated.
• Box method used to calculate the fractal
  dimension.
• The effects of growth time and fractal
  dimension can be studied.

12
def walk ()     Checks to see if particle
became attached from the movement of another
particle,     If not, moves the particles,    
Checks for attachment     Only checks if
particle is within the radius of the existing
structure def move(point)     Move the
particle randomly,     The grid is cyclic -
going off one side puts it on the other def
check(point)     Checks to see if the particle
is touching any particle in the aggregate   
Changes the radius if necessary def fractDim ()
    Uses boxes to calculate box-dimension    
First box has smallest dimensions that enclose
the aggregate     Subsequent boxes are 1/8 the
size     Counts the number of boxes necessary to
cover the aggregate at each size     Plots the
common log of box diameter versus the common log
of the count     Calculates best fit linear
regression line (R2gt0.999)     Negative of the
slope is the box-dimension
13
Some Mathematical Methods

• Euler Methods
• Improved Euler

14

• Runge-Kutta 4th order

15
Poincare Maps

• Choose any three from x, y, vx, and vy and
  calculate the fourth using the integral of the
  system, such as x, y, and vy.
• Follow a solution in the three-dimensional space
  of these three coordinates. If the surface is the
  plane x 0, then we are looking at points that
  cross in the y-vy plane.
• Now let a solution cross this plane at the point
  P, where x is increasing or vx gt 0. The next
  crossing occurs at Q, where vx lt 0. The next
  cross occurs when x is increasing, at R.

16
Poincaré map for the Hénon-Heiles problem with E
0.1053, x(0)0.4, y(0) 0, vx(0) 0, vy(0)
0.225
Poincaré map for the Hénon-Heiles problem with E
0.1503, x(0)0.5, y(0) 0, vx(0) 0, vy(0)
0.225
17
Advantages of VPython

• The semester-long experiment using VPython in a
  computational physics course proved highly
  successful. The notable successes included
• Shorter learning curve for programming.
• 3D visualization.
• Correspondence between vector laws and object
  oriented programming.
• More in-depth student projects.
• Student enjoyment with developing simulations.
• Carry over of modeling into other courses.

18
Gravitational Shaping by Galactic Encounters
As a galaxy encounters another galaxy
(represented as a point mass), the orbits of
individual stars are perturbed.
19
AvalanchesAn Example of Self-Organized
Criticality
Avalanches occur when a critical threshold is
exceeded.
20
def fall(x,y,D)     global grid, sand    
gridxD3-1y1-D/301 adds one to
receiving square     sgridxy1.pop()
index of sand being moved    
gridxD3-1y1-D/31.append(s) adds index
to receiving square     gridxy0-1
removes one from falling square    
sands.pos(xD3-1,y1-D/3,gridxD3-1y1-D/3
0) relocates sand
21
Percolation Growth of Galaxies
In this percolation model of galaxies, long range
spiral structure is achieved through local
interactions.
22
def growth()     Each inhabited sector has a p
chance of giving birth to a new star, in which
case it is destroyed                             
def rotate()     Rotates each sector
according to the angular velocity def spiral()
    Determines if stars are born or have died in
each sector, updates image to reflect current
population and position
23
3D Dynamic Lattice
A 3D lattice of particle connected by springs
allows for an investigation of normal modes.
24
sum the forces acting on the current
particle Forces (Force(i,j,k,k,j,idim)
Force(i,j,k,k,j,iaug)
Force(i,j,k,k,jdim,i)
Force(i,j,k,k,jaug,i)
Force(i,j,k,kdim,j,i)
Force(i,j,k,kaug,j,i)) Forces
-(masskji.velocityb) if dropping
add gravity and handle collisions Forces
mgvector(0,1,0) if masskji.pos.y lt
floor if masskji.velocity.y lt 0
masskji.velocity.y
-masskji.velocity.y if masskji.fixed
do no move if particle is fixed in place
Forces vector(0,0,0) masskji.velocity
vector(0,0,0) cur_Forces Forces set
variables for passing to rkf4 var0masskji
.pos var1masskji.velocity use RKF4 to
evaluate differential equation
25
Non-relativistic Electron Scattering
This project involves a Monte Carlo simulation of
non-relatavistic electron scattering in a thin
film. The backscattering, transmission, and
absorption can be studied for various metals.
26
calcluate new positions for particles
for n in range(0, maxn)
c1.interact() theta.append(0) add
the n1 element phi.append(0) add
the n1 element s.append(0)
initialize position element to 0
T.append(0) R.append(vector(0,0,0))  
Rtheta random() Rphi
random() Rs random()  
p_costheta 1-(2betaRtheta /
(1beta-Rtheta)) p_theta
acos(1-(2betaRtheta / (1beta-Rtheta)))
p_theta 1 - p_theta p_phi
2piRphi   costhetanew
cos(thetan)cos(p_theta)sin(thetan)sin(p_the
ta)cos(p_phi)   thetan1
acos(costhetanew) top (
sin(thetan)sin(p_theta)sin(p_phi) )
bottom ( cos(p_theta)-(costhetanew)cos(t
hetan))   if bottom 0
bottom 1 top
sin(p_theta)sin(p_phi) phin1
atan( top / bottom ) phin sn
(-1.02beta(beta1)A(Tn)2 /
(Z(Z1)rho)) log(Rs) Tn1
Tn - abs(dTds(Tn))sn Rn1
Rn sn(sin(thetan1)cos(phin1)
vector(1,0,0)
sin(thetan1)sin(phin1)vector(0,1,0)

cos(thetan1)vector(0,0,1))
27
Follow-Up
For additional information please contact
olenick_at_udallas.edu.