One dimension kinematics

movement.py

#!/usr/bin/env python 
'''one dimension kinematics
to be run from:
$ ipython -pylab
[1]: run movement.py
'''
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import sys

def fun(c0=8,c1=6,c2=-1,c3=0,c4=0,c5=0):
    #lambda t: 8 -6.*t+10*t**2-t**3
    return lambda t: c0+c1*t+c2*t**2+c3*t**3+c4*t**4+c5*t**5

def dynamics(x,Dt,m=1):
    '''For the worldline defined by the ndarray x, calculates
    Dt, 
    and ndarrays: 
      v, a, T, V, E and S
    for each segment'''
    
    x=np.asarray(x)
    v=(x[1:]-x[:-1])/Dt
    a=(v[1:]-v[:-1])/Dt
    T=0.5*m*v**2
    V=0.5*m*a[-1]*(x[1:]+x[:-1])
    E=T+V
    SS=(T-V)*Dt
    return v,a,T,V,E,SS

def vplot(t,x,v,step=1,dot='ro'):
    plt.plot(t,x)
    plt.plot(t,x,dot)
    for i in range(0,v.size,step):
        plt.annotate(str(v[i]),(t[i+1],x[i+1]))

    plt.ylabel('x (m)',size=20)
    plt.xlabel('t (s)',size=20)
    
def tlist(tmin=0,tmax=10,Dt=1):
    tmin=0
    tmax=10.
    return np.arange(tmin,tmax,Dt),Dt

if __name__ == "__main__":
    x=fun()
    (t,Dt)=tlist(Dt=1)
    (v,a,T,V,E,SS)=dynamics(x(t),Dt)
    vplot(t,x(t),v)
    plt.title(r'$x(t)=8+6 t-t^2$')
    print "v",v
    print "a",a
    plt.show()
    s=raw_input('continue')
    if s=='s':  sys.exit()
    (t,Dt)=tlist(Dt=0.1)
    (v,a,T,V,E,SS)=dynamics(x(t),Dt)
    vplot(t,x(t),v,10)
    plt.title(r'$x(t)=8+6 t-t^2$')
    print "v",v[0]
    print "a",a[0]
    s=raw_input('continue')
    if s=='s':  sys.exit()
    (t,Dt)=tlist(Dt=0.01)
    (v,a,T,V,E,SS)=dynamics(x(t),Dt)
    vplot(t,x(t),v,100,dot='r.')
    plt.title(r'$x(t)=8+6 t-t^2$')
    print "v",v[0]
    print "a",a[:10]
    #print "E",E[:10]
    plt.plot(t[:-1],v,'b.')
    plt.plot(t[:-2],a,'g.')

    s=raw_input('continue')
    if s=='s':  sys.exit()
    p=[-0.32517483,0.8240676,2.10067016,-0.89087995,0.12194056,-0.00544872]
    x=fun(p[0],p[1],p[2],p[3],p[4],p[5])
    (t,Dt)=tlist(Dt=0.01)
    (v,a,T,V,E,SS)=dynamics(x(t),Dt)
    vplot(t,x(t),v,100,dot='r.')
    print "v",v[0]
    print "a",a[0:10]
    #print "E",E[:10]
    plt.plot(t[:-1],v,'b.')
    plt.plot(t[:-2],a,'g.')
    plt.hlines(0,0,10)