[http://bit.ly/aaptw15titus](http://bit.ly/aaptw15titus)

Aaron Titus

atitus@highpoint.edu, [@physicslogos](http://twitter.com/physicslogos)

High Point University



In [1]:

    
# I keep this as a cell in my title slide so I can rerun 
# it easily if I make changes, but it's low enough it won't
# be visible in presentation mode.
%run talktools

Links

Presentation. http://nbviewer.ipython.org/github/atitus/PHY3110/blob/master/presentation/AAPTw15/AAPTw15.ipynb
Course on Trinket. https://trinket.io/atitus/courses/phy-3110-advanced-classical-mechanics#/
Github respository. https://github.com/atitus/PHY3110
Instructions for installing ipython notebook. https://trinket.io/atitus/courses/phy-3110-advanced-classical-mechanics#/tech-tools/installing-python-and-ivisual
iVisual. https://pypi.python.org/pypi/IVisual
Plot.ly. http://nbviewer.ipython.org/gist/msund/61cdbd5b22c103fffb84

Outline

**An experiment** to flip an advanced classical mechanics course.
**Trinket** for displaying course content.
**iPython and iVisual** for computational models and numerical calculations.

Goals of the Course

My explicit goal is to help students adopt a Growth Mindset. To achieve this goal, I

gave students **a problem each day** to solve during class-time in collaboration with others.
incorporated a few **experiments**.
moved derivations, class notes, and examples to **Trinket**.
used **iPython Notebook** (often with iVisual) to create computational models, visualize motion, solve differential equations numerically, and graph data for lab reports.
applied standards based grading (**SBG**).

I refer to the daily problems as "modified problem based learning (PBL)" because the problems are typically doable in one hour and are indicative of homework problems and exam problems. (Typically PBL refers to more comprehensive problems that may be ill-defined and may require research, estimation, collaboration, and no clear answer.)

Trinket



In [2]:

    
website('https://trinket.io/atitus/courses/phy-3110-advanced-classical-mechanics#/', 'Trinket Course')









    Out[2]:





Trinket Course

Central Force



In [24]:

    
website('https://atitus.trinket.io/phy-3110-advanced-classical-mechanics#/chapter-04/example-particle-barbell-sticky-collision','Sticky Collision of a Particle and Barbell')









    Out[24]:





Sticky Collision of a Particle and Barbell



In [4]:

    
from __future__ import division, print_function
from ivisual import *
from math import *



In [5]:

    
scene=canvas(title="Binary Stars")

#constants
au=1.5e11
R=7e8
G=6.67e-11

#stars
m1=2e30
m2=2*m1
r1=vector(au,0,0)
r2=-m1/m2*r1
speed=1.1*3e4
v1=vector(0,speed,0)
v2=-m1/m2*v1

#time
t=0
dt=1e5

#data lists
Ulist=[]
Klist=[]
Elist=[]
rlist=[]

#ivisual objects
star1=sphere(radius=10*R, pos=r1, color=color.red, make_trail=True, retain=400)
star2=sphere(radius=10*R, pos=r2, color=color.yellow, make_trail=True, retain=400)
CM=sphere(radius=star1.radius/2, pos=(m1*r1+m2*r2)/(m1+m2), color=color.white)

while t<1000*dt:
    rate(100)
    r=r1-r2
    rmag=mag(r)
    runit=norm(r)
    Fgrav=G*m1*m2/rmag/rmag
    F1=Fgrav*(-runit)
    v1=v1+F1/m1*dt
    r1=r1+v1*dt
    F2=-F1
    v2=v2+F2/m2*dt
    r2=r2+v2*dt
    
    U=-G*m1*m2/rmag
    K=1/2*m1*mag(v1)*mag(v1)+1/2*m2*mag(v2)*mag(v2)
    E=U+K
    
    Ulist.append(U)
    Klist.append(K)
    Elist.append(E)
    rlist.append(rmag)
    
    star1.pos=r1
    star2.pos=r2
    
    t=t+dt



In [6]:

    
%matplotlib inline



In [7]:

    
import matplotlib.pyplot as plt



In [8]:

    
#energy graphs
plt.title('energy vs. distance')
plt.ylabel('E (J)')
plt.xlabel('r (m)')
plt.plot(rlist,Ulist,'m.', rlist,Klist,'y.', rlist,Elist,'c.')
plt.show()



In [9]:

    
scene2=canvas(title="One Star Model")

#reset variables
r1=vector(au,0,0)
r2=-m1/m2*r1
v1=vector(0,speed,0)
v2=-m1/m2*v1

#one star model
mu=m1*m2/(m1+m2)
L=m1*cross(r1,v1)+m2*cross(r2,v2)
Lmag=mag(L)
r=r1-r2
rmag=mag(r)
rdot=0
phi=0
phidot=Lmag/rmag/rmag/mu
v=vector(rdot,rmag*phidot)

#time
t=0

#data lists
Ulist=[]
Klist=[]
Elist=[]
rlist=[]

#ivisual objects
starmu=sphere(radius=20*R, pos=r, color=color.orange, make_trail=True, retain=400)

while t<1000*dt:
    rate(100)
    
    #update rmag
    rdotdot=(-G*m1*m2/(rmag*rmag)+Lmag*Lmag/(mu*rmag*rmag*rmag))/mu
    rdot=rdot+rdotdot*dt
    rmag=rmag+rdot*dt

    #update phi
    phidot=Lmag/rmag/rmag/mu
    phi=phi+phidot*dt
    
    #calculate velocity and position
    v=vector(rdot,rmag*phidot)
    r=vector(rmag*cos(phi),rmag*sin(phi),0)
    
    U=-G*m1*m2/rmag + Lmag*Lmag/(2*mu*rmag*rmag)
    K=1/2*mu*rdot*rdot
    E=U+K
    
    Ulist.append(U)
    Klist.append(K)
    Elist.append(E)
    rlist.append(rmag)

    #ivisual objects
    starmu.pos=r
    
    t=t+dt



In [10]:

    
#energy graphs
plt.title('energy vs. distance')
plt.xlabel('r (m)')
plt.ylabel('E (J)')
plt.plot(rlist,Ulist,'m.', rlist,Klist,'y.', rlist,Elist,'c.')
plt.show()



In [12]:

    
scene3=canvas(title="One Star Model with Two Stars Calculated")

#reset variables
r1=vector(au,0,0)
r2=-m1/m2*r1
v1=vector(0,speed,0)
v2=-m1/m2*v1

#one star model
mu=m1*m2/(m1+m2)
L=m1*cross(r1,v1)+m2*cross(r2,v2)
Lmag=mag(L)
r=r1-r2
rmag=mag(r)
rdot=0
phi=0
phidot=Lmag/rmag/rmag/mu
v=vector(rdot,rmag*phidot)

#time
t=0

#data lists
Ulist=[]
Klist=[]
Elist=[]
rlist=[]

#ivisual objects
starmu=sphere(radius=20*R, pos=r, color=color.orange, make_trail=True, retain=400)
star1=sphere(radius=10*R, pos=r1, color=color.red, make_trail=True, retain=400)
star2=sphere(radius=10*R, pos=r2, color=color.yellow, make_trail=True, retain=400)
CM=sphere(radius=star1.radius/2, pos=(m1*r1+m2*r2)/(m1+m2), color=color.white)


while t<1000*dt:
    rate(100)
    
    #update rmag
    rdotdot=(-G*m1*m2/(rmag*rmag)+Lmag*Lmag/(mu*rmag*rmag*rmag))/mu
    rdot=rdot+rdotdot*dt
    rmag=rmag+rdot*dt

    #update phi
    phidot=Lmag/rmag/rmag/mu
    phi=phi+phidot*dt
    
    #calculate velocity and position
    v=vector(rdot,rmag*phidot)
    r=vector(rmag*cos(phi),rmag*sin(phi),0)
    
    U=-G*m1*m2/rmag + Lmag*Lmag/(2*mu*rmag*rmag)
    K=1/2*mu*rdot*rdot
    E=U+K
    
    Ulist.append(U)
    Klist.append(K)
    Elist.append(E)
    rlist.append(rmag)

    #ivisual objects
    starmu.pos=r
    star1.pos=m2/(m1+m2)*r
    star2.pos=-m1/(m1+m2)*r
    
    t=t+dt

LAB -- Rocket Accelerated Disk



In [13]:

    
from IPython.display import HTML
from base64 import b64encode
video = open("video_disk_accelerated_rocket.mov", "rb").read()
video_encoded = b64encode(video)
video_tag = '<video controls alt="test" src="data:video/x-m4v;base64,{0}">'.format(video_encoded)
HTML(data=video_tag)









    Out[13]:



In [14]:

    
#t (s) and omega (deg/s)

data = """
0.033	7.9
0.067	22.8
0.1	45.8
0.133	76.1
0.167	109.7
0.2	147.7
0.234	193.2
0.267	232.3
0.3	254.2
0.334	274.7
0.367	297.4
0.4	318.1
0.434	335.2
0.467	349.1
0.501	371.4
0.534	400.3
0.567	421.5
0.601	437.5
0.634	460.4
0.667	472.7
0.701	491.8
0.734	518.9
0.767	533.3
0.801	550.5
0.834	561.8
0.868	570.0
0.901	578.5
0.934	589.7
0.968	593.6
1.001	589.8
1.034	606.7
1.068	624.6
1.101	638.9
1.134	666.8
1.168	698.0
1.201	735.5
1.235	767.4
1.268	784.3
1.301	800.4
1.335	815.7
1.368	828.9
1.401	827.2
1.435	813.1
1.468	811.8
1.502	806.2
1.535	816.9
1.568	845.8
1.602	883.6
1.635	938.0
1.668	981.7
1.702	1012.1
1.735	1024.3
1.768	1027.8
1.802	1006.6
1.835	978.3
1.869	959.8
1.902	957.7
1.935	992.7
1.969	1052.1
2.002	1121.2
2.035	1169.7
2.069	1183.1
2.102	1157.3
2.135	1112.8
2.169	1065.2
2.202	1027.6
2.236	1016.7
2.269	1039.4
2.302	1088.5
2.336	1145.0
2.369	1181.2
2.402	1180.0
2.436	1154.2
2.469	1104.6
2.503	1048.4
2.536	1014.9
2.569	1017.1
2.603	1050.6
2.636	1094.9
2.669	1144.3
2.703	1178.7
2.736	1173.7
2.769	1142.8
2.803	1087.3
2.836	1045.6
2.87	1012.4
2.903	1005.5
2.936	1047.9
2.97	1096.2
3.003	1150.2
3.036	1171.1
3.07	1155.3
3.103	1130.6
3.136	1080.1
3.17	1029.0
3.203	1003.7
3.237	1004.3
""".split('\n')  # split this string on the "newline" character.

print(len(data))



In [15]:

    
#
# Here we'll take the list of strings defined above and break it into actual numbers in 
# mks units.
#

tlist = []
omegalist = []
for s in data:
    if s:
        t,omega = s.split()     # break string in two
        t=float(t)          # convert time to float
        omega=float(omega)/180*pi #convert omega to float
        tlist.append(t)
        omegalist.append(omega)
        
#print "tlist=",tlist
#print "omegalist=",omegalist



In [16]:

    
from numpy import *
omegatheor=1.465*sin(19.21*array(tlist)+5.905)+19.12

fig1 = plt.figure()
plt.title('omega vs t')
plt.xlabel('t (s)')
plt.ylabel('omega (rad/s)')
plt.plot(tlist, omegalist,'b.',tlist,omegatheor,'m-')
plt.show









    Out[16]:





<function matplotlib.pyplot.show>



In [17]:

    
import plotly.plotly as py
import plotly.tools as tls
from plotly.graph_objs import *
py.sign_in("hpuphysics", "fwxi6w237k")
py.iplot_mpl(fig1, strip_style = True)









    Out[17]:

Computer Model



In [18]:

    
from IPython.display import Image
Image(filename='disk-final.png')









    Out[18]:

$$ \frac{d\omega_{A,z}}{dt} = \ddot{\phi} = \frac{ F_{thrust}R -(m_1+m_2)gd\cos(\phi)}{(I_{disk}+I_1+I_2)} $$$$ d = \left(\frac{m_1-m_2}{m_1+m_2}\right) R \\ F_{thrust} =-\dot{m}v_{exhaust} $$

Computer Model After Burnout



In [19]:

    
#disk
M=2.94
R=0.6048/2
m1=0.52
m2=0.15*m1
d=(m1-m2)/(m1+m2)*R
I=R*R*(1/2*M+m1+m2)
print("d = ",d)

#initial angle with respect to -y axis
phi_deg=180
phi=phi_deg*pi/180 #phi in rad

#initial angular velocity
phidot=19  #rad/s

#gravitational field in N/kg
g=9.8

t=0
dt=0.005

#position vector
r=d*vector(sin(phi),-cos(phi),0)

#period
T=0
ti=0
atpeak=False

tlist=[]
omegalist=[]
philist=[]

while t<1.25:   
    omegai=phidot
    
    #update phidotdot, phidot, phi
    phidotdot=-(m1+m2)*g*d/I*cos(phi)
    phidot=phidot+phidotdot*dt
    phi=phi+phidot*dt
    
    #update r
    r=R*vector(sin(phi),-cos(phi),0)
    
    #find period
    if(phidotdot>0 and (phidot-omegai)<0.001 and phidot>19 and atpeak==False):
        ti=t
#        print("t_i = ",t)
        atpeak=True
    if(phidotdot>0 and (phidot-omegai)<0.001 and phidot>19 and atpeak==True):
        T=t-ti
        ti=t
#        print("T = ",T)
    
    
    tlist.append(t)
    omegalist.append(phidot)
    philist.append(phi)
    
    t=t+dt

plt.title('omega vs. t')
plt.xlabel('t (s)')
plt.ylabel('omega (rad/s)')
plt.plot(tlist,omegalist,'m.')
plt.show()









    



d =  0.223513043478

Computer Model with Rocket Burning



In [20]:

    
from ivisual import *

g=9.8
mdisk=2.94 #kg
Rdisk=0.6048/2 #m
m0engine=520/1000 #kg
mfuel=0.85*m0engine #kg
vex=30 #m/s
mf=m0engine-mfuel
m=m0engine
m1=m0engine

r=0.3 #m
mdot=-mfuel/2 #2 s burn time

t=0
dt=0.001
theta=0
omega=0
Idisk=1/2*mdisk*Rdisk*Rdisk
Ieng1=m1*r*r
Ieng2=m*r*r
I=Idisk+Ieng1+Ieng2

r1=r*vector(cos(theta+pi),sin(theta+pi),0)
r2=r*vector(cos(theta),sin(theta),0)

tlist=[]
thetalist=[]
omegalist=[]
mlist=[]
taulist=[]

tau=vector(0,0,0)

while t<3:
    if(m>mf):
        m=m+mdot*dt
    Fgrav1=m1*vector(0,-g,0)
    Fgrav2=m*vector(0,-g,0)
    thrust=-mdot*vex
    
    taugrav1=cross(r1,Fgrav1)
    taugrav2=cross(r2,Fgrav2)
    if(m>mf):
        tauthrust=vector(0,0,thrust*r)
    else:
        tauthrust=vector(0,0,0)
    tau=taugrav1+taugrav2+tauthrust
    
    Ieng2=m*r*r
    I=Idisk+Ieng1+Ieng2
    
    omega=omega+tau.z/I*dt
    theta=theta+omega*dt
    
    r1=r*vector(cos(theta+pi),sin(theta+pi),0)
    r2=r*vector(cos(theta),sin(theta),0)
    
    t=t+dt
    
    tlist.append(t)
    thetalist.append(theta)
    omegalist.append(omega)
    mlist.append(m)
    taulist.append(tau.z)



In [21]:

    
plt.title('omega vs. t')
plt.xlabel('t (s)')
plt.ylabel('omega (rad/s)')
plt.plot(tlist,omegalist,'m.')
plt.show()

iPython for Homework

A simple pendulum of length $b$ and mass $m$ is suspended in a gravitational field $g$ from a point on the circumference of a thin low-mass disc of radius $R$ that rotates counterclockwise with a constant angular velocity $\omega$ about is central axis as shown below.

Find a differential equation for $\theta$, the angle the pendulum makes with the vertical. Call downward $\theta=0$ and angles to the right as positive.
Write an iPython Notebook to model the motion of the pendulum. Graph $\theta(t)$, $x(t)$, and $y(t)$ for the pendulum, where $x$ and $y$ are the Cartesian positions of the mass $m$. You can select reasonable values for the constants $b$, $R$, $m$, $\omega$, $\theta_0$ and $\dot{\theta}_0$.



In [22]:

    
Image(filename='fig-1.png')









    Out[22]:



In [23]:

    
scene10=canvas("Pendulum hanging from rotating pivot")

omega=2*pi/2
g=9.8
b=1
R=b/5

phi=0
theta=30*pi/180
thetadot=0

t=0
dt=0.005

rpivot=R*vector(cos(phi),sin(phi),0)
bvec=b*vector(sin(theta),-cos(theta),0)
r=bvec+rpivot

wheel=cylinder(pos=vector(0,0,-R/10), axis=vector(0,0,R/10), radius=R, color=color.red)
pivot=sphere(pos=rpivot, radius=R/20, color=color.yellow)
pendulum=sphere(pos=r, radius=R/4, color=color.yellow, make_trail=True, retain=2000)
string=cylinder(pos=pivot.pos, axis=bvec, radius=R/20, color=color.yellow)

while t<25:
    rate(100)
    
    thetadotdot=-omega**2*R/b*cos(theta+omega*t)-g/b*sin(theta)
    thetadot=thetadot+thetadotdot*dt
    theta=theta+thetadot*dt
    
    phi=phi+omega*dt
    
    rpivot=R*vector(cos(phi),sin(phi),0)
    bvec=b*vector(sin(theta),-cos(theta),0)
    r=bvec+rpivot
    
    pivot.pos=rpivot
    pendulum.pos=r
    string.pos=rpivot
    string.axis=r-rpivot
    
    t=t+dt



In [ ]: