#!/usr/bin/env python
'''
vhull.py
Eric Ayars
June 2016
A solution to the vertical motion of a V-hulled boat.
'''
from pylab import *
from scipy.integrate import odeint
# definitions
w = 1.0 # Well, we need to pick something. Might as well
# normalize the frequency. :-)
delta = 0.2 # not sure what this should be, play with it!
def vhull_ODE(s,t):
# derivs function for v-hull ODE, without viscous damping
g0 = s[1]
g1 = -w*(1+s[0]/2)*s[0]
return array([g0,g1])
def vhull_ODE2(s,t):
# derivs function for v-hull ODE, with viscous damping added
g0 = s[1]
g1 = -w*(1+s[0])*s[0] - delta*s[1]
return array([g0,g1])
# Now set up the solutions.
time = linspace(0,20,500) # 500 points on interval [0,20]
# And now plot some things!
# Here's a low-amplitude one, should be approximately SHO.
A = 0.1 # Fairly small amplitude
yi = array([A,0.0])
answer = odeint(vhull_ODE, yi, time)
y = answer[:,0] # position is just the first row
plot(time,y, 'b-', label="Actual solution")
plot(time,A*cos(w*time), 'r-', label="SHO approximation")
A = 0.5 # Not small amplitude
yi = array([A,0.0])
answer = odeint(vhull_ODE, yi, time)
y = answer[:,0] # position is just the first row
plot(time,y, 'b-')
plot(time,A*cos(w*time), 'r-')
'''
# For damped motion, just use vhull_ODE2() instead of vhull_ODE.
A = 0.5
yi = array([A,0.0])
answer = odeint(vhull_ODE2, yi, time)
y = answer[:,0] # position is just the first row
plot(time,y, 'b-', label="Actual solution")
plot(time,A*cos(w*time), 'r-', label="SHO approximation")
'''
title("V-Hull Boat Solution")
xlabel("time")
ylabel("Position")
legend(loc="lower right")
show()