Damped driven pendulum DE

solve.py

from numpy import sin, cos
import scipy
from scipy.integrate import odeint
from scipy import arange, array

def func(y, t, L, S, g, c, scale):
    theta = y[0]
    d_theta_dt = y[1]
    return (d_theta_dt, -g/L*sin(theta)-scale*S(t)/L*cos(theta)-c*d_theta_dt) 

def solve_de(L, S, g, t, c, scale = 1.0):
    u0 = array([0.,0.])
    u = odeint(func, u0, t, args = (L,S,g,c,scale))
    return u

def radial_force(theta, m, g, S, t, scale=1.0):
    return(m*g*cos(theta)+scale*m*S(t)*sin(theta))

if __name__ == '__main__':
    from matplotlib import *
    from pylab import *
    import read_rscth
    import interpolate
    f = open('rscth.out', 'r')
    print('reading data')
    data = read_rscth.read(f)
    print('interpolate')
    S = interpolate.interpolate(data[:,0], data[:,1])
    m = 10.0
    L = 1.00
    g = 9.81
    c = 0.05
    t = arange(0.0,50.0,0.0001)
    with open('time.csv', 'wb') as f:
        for time in t:
            f.write('%s\n' % (time,))
    with open('acc.csv', 'wb') as f:
        acc = S(t)
        print('Maximum acceleration (cm/s^2):%s' % (scipy.amax(scipy.absolute(acc)),))
        for a in acc:
            f.write('%s\n' % (a,))

    u = solve_de(L, S, g, t, c, scale=0.01)
    force = radial_force(u[:,0], m, g, S, t, scale=0.01)
    with open('force_vs_t.csv', 'wb') as f:
        for time, r in zip(t,force):
            f.write('%s,%s\n' % (time,r))
    print('Maximum cable force (N):%s' % (scipy.amax(scipy.absolute(force)),))
    '''
    figure(1)
    plot(t,u[:,0],label='theta')
    plot(t,u[:,1],label='d(theta)/dt')
    xlabel('Time (s)')
    title('Damped driven pendulum')
    legend()
    figure(2)
    plot(u[:,0],u[:,1])
    title('phase-space')
    xlabel('Theta (rad)')
    ylabel('d(theta)/dt')
    figure(3)
    plot(t, force)
    title('Force in cable')
    xlabel('Time (s)')
    ylabel('Force (N)')
    figure(4)
    title('Force in cable')
    xlabel('Position (rad)')
    ylabel('Force (N)')
    plot(u[:,0], force)
    figure(5)
    title('Time Series of Seismic Acceleration')
    xlabel('Time (s)')
    ylabel('Horizontal Acceleration (cm/s^2)')
    plot(t, acc)
    show()
    '''