#!/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()