Damped driven pendulum

Fp_dynamic.py

from numpy import sqrt, cos, sin, exp, arange, amax, absolute
from math import pow
from matplotlib import *
from pylab import *

def theta_func(L, g, m, x):
    tmp = (179.5*pow(g,3./2.)*L*sin((sqrt(g)*x)/sqrt(L))*sin(x*(3.132-sqrt(g)/sqrt(L)))-\
           179.5*pow(g,3./2.)*L*cos((sqrt(g)*x)/sqrt(L))*cos(x*(3.132-sqrt(g)/sqrt(L)))-\
           5514.8*pow(L,5./2.)*sin((sqrt(g)*x)/sqrt(L))*sin(x*(3.132-sqrt(g)/sqrt(L)))+\
           5514.8*pow(L,5./2.)*cos((sqrt(g)*x)/sqrt(L))*cos(x*(3.132-sqrt(g)/sqrt(L)))+\
           562.194*g*pow(L,3./2.)*sin((sqrt(g)*x)/sqrt(L))*sin(x*(3.132-sqrt(g)/sqrt(L)))-\
           562.194*g*pow(L,3./2.)*cos((sqrt(g)*x)/sqrt(L))*cos(x*(3.132-sqrt(g)/sqrt(L)))-\
           1760.79*sqrt(g)*pow(L,2)*sin((sqrt(g)*x)/sqrt(L))*sin(x*(3.132-sqrt(g)/sqrt(L)))+\
           1760.79*sqrt(g)*pow(L,2)*cos((sqrt(g)*x)/sqrt(L))*cos(x*(3.132-sqrt(g)/sqrt(L)))-\
           179.5*pow(g,3./2.)*L*sin(x*(sqrt(g)/sqrt(L)+3.132))*sin((sqrt(g)*x)/sqrt(L))+\
           359.*pow(g,3./2.)*L*cos((sqrt(g)*x)/sqrt(L))-\
           179.5*pow(g,3./2.)*L*cos(x*(sqrt(g)/sqrt(L)+3.132))*cos((sqrt(g)*x)/sqrt(L))-\
           5514.8*pow(L,5./2.)*sin(x*(sqrt(g)/sqrt(L)+3.132))*sin((sqrt(g)*x)/sqrt(L))-\
           5514.8*pow(L,5./2.)*cos(x*(sqrt(g)/sqrt(L)+3.132))*cos((sqrt(g)*x)/sqrt(L))+\
           562.194*g*pow(L,3./2.)*sin(x*(sqrt(g)/sqrt(L)+3.132))*sin((sqrt(g)*x)/sqrt(L))+\
           562.194*g*pow(L,3./2.)*cos(x*(sqrt(g)/sqrt(L)+3.132))*cos((sqrt(g)*x)/sqrt(L))+\
           1760.79*sqrt(g)*pow(L,2)*sin(x*(sqrt(g)/sqrt(L)+3.132))*sin((sqrt(g)*x)/sqrt(L))+\
           1760.79*sqrt(g)*pow(L,2)*cos(x*(sqrt(g)/sqrt(L)+3.132))*cos((sqrt(g)*x)/sqrt(L))-\
           3521.58*sqrt(g)*pow(L,2)*cos((sqrt(g)*x)/sqrt(L)))/\
        (sqrt(g)*m*(3.132*sqrt(L)-sqrt(g))*(sqrt(g)+3.132*sqrt(L))*(9.80942*L-g))

    return tmp
L = 1.0
g = 9.81
m = 10.0
t = arange(0.,20.,0.001)
theta = theta_func(L,g,m,t)
force = m*g*cos(theta)-359.0*sin(3.132*theta)
print('Maximum cable force (N):%s' % (amax(absolute(force)),))
figure(1)
plot(t,theta,label='theta')
xlabel('Time (s)')
title('Damped driven pendulum')
legend()
figure(2)
plot(t,force,label='force')
xlabel('Time (s)')
title('Damped driven pendulum')
legend()
show()