Figures out the mass of your narrowboat and how much ballast you need to add to achieve a desired result.

ballast.py

from __future__ import division
import numpy as np

# A boat is a dict with keys
# mass, center_of_mass, cm_depth, front_depth, back_depth, pitch, length, width
# units are radians (angles), meters (depths and lengths), tons (mass)
# Boats are idealized as rectangular prisms

water_density = 1

kg_to_pound = 2.2
ton_to_pound = kg_to_pound*1000
inch_to_meter = .0254
meter_to_inch = 1/inch_to_meter

def buoyancy_constraint(center_of_mass, cm_depth, length, pitch):
    """
    Gives zero for the correct center of mass.
    
    This is the torque about the center of mass due to buoyancy, found by \int_0^l P(x)(x-cm)dx
    where P is pressure. Since the hydrostatic force is normal to the baseplate, there's no need
    for a cos(theta) term. The pressure at position x is equal to d_cm  - (x-cm)\sin \theta where
    d_cm is the depth at the cm and \theta is the pitch angle.
    """
    l = length-center_of_mass
    return (l**2-center_of_mass**2)/2.*cm_depth - (l**3+center_of_mass**3)/3.*np.sin(pitch)

def mass_from_displacement(center_of_mass, cm_depth, length, pitch, width, water_density=water_density):
    """
    Computes the mass from the cm, the depth at cm, pitch, length and width
    of the boat, and the density of water.
    """
    back_depth = cm_depth + center_of_mass*np.sin(pitch)
    front_depth = back_depth - length*np.sin(pitch)
    return (front_depth+back_depth)/2.*length*width*np.cos(pitch)*water_density
    
def init_boat(front_depth, back_depth, length, width):
    "Computes the mass, cm and attitude of the boat from readily-measurable data."
    if np.abs(front_depth-back_depth) > length:
        raise ValueError, 'Difference between front and back depths exceeds length.'
    
    # Compute the pitch angle first.
    pitch = np.arcsin((back_depth-front_depth)/length)
    
    # Find the CM by satisfying the buoyancy constraint.
    from scipy import optimize
    f = lambda cm, bd=back_depth, l=length, p=pitch: buoyancy_constraint(cm, bd - cm*np.sin(p), l, p)
    center_of_mass = optimize.bisect(f, 0, length)
    cm_depth = (front_depth*(center_of_mass)+back_depth*(length-center_of_mass))/length
    
    # Everything else is easy once the CM is known...
    mass = mass_from_displacement(center_of_mass, cm_depth, length, pitch, width)
    
    return {'mass': mass,
            'center_of_mass': center_of_mass,
            'cm_depth': cm_depth,
            'front_depth': front_depth,
            'back_depth': back_depth,
            'pitch': pitch, 
            'length': length,
            'width': width}
            
def plot_boat(boat):
    "Plots the attitude of the boat to scale, with the waterline shown."
    import pylab as pl
    pl.plot([-boat['center_of_mass']*np.cos(boat['pitch']),
                (boat['length']-boat['center_of_mass'])*np.cos(boat['pitch'])],
            [-boat['back_depth'],-boat['front_depth']],'r-')
    pl.plot([-boat['center_of_mass'],boat['length']-boat['center_of_mass']],[0,0],'b-')
    pl.axis('image')
    pl.axis('off')

def ballast(boat, new_mass, new_center_of_mass):
    """
    Adds a new mass distribution with given mass and center of mass to the boat.
    Returns a new boat with the requested mass added.
    """
    # Get new mass and CM. Remember, this isn't like adding two RV's, 
    # it's like mixing two distributions.
    mass = boat['mass']+new_mass
    center_of_mass = (boat['center_of_mass']*boat['mass']+new_center_of_mass*new_mass)/mass
    
    # These don't change.
    width=boat['width']
    length=boat['length']
    
    # A function of depth at cm and pitch angle
    # which is only zero when the mass is correct
    # and the buoyancy constraint is satisfied.
    def f(x, cm=center_of_mass, m=mass, w=width, l=length):
        cm_depth=x[0]
        pitch=x[1]
        md = mass_from_displacement(cm, cm_depth, l, pitch, w)
        bc = buoyancy_constraint(cm, cm_depth, l, pitch)
        return (md-m)**2+bc**2
    
    # Find the depth at cm and the pitch angle.    
    from scipy import optimize
    x = optimize.fmin(f, np.array([boat['cm_depth'], boat['pitch']]))#, bounds=[(0,None),(-np.pi/2.,np.pi/2.)])
    cm_depth, pitch = x
    
    # Everything else is easy...
    back_depth = cm_depth+center_of_mass*np.sin(pitch)
    front_depth = back_depth-length*np.sin(pitch)

    return {'mass': mass,
            'center_of_mass': center_of_mass,
            'cm_depth': cm_depth,
            'front_depth': front_depth,
            'back_depth': back_depth,
            'pitch': pitch, 
            'length': length,
            'width': width}
            
def find_ballast_amount(f, m0, ballast_center_of_mass):
    """
    Returns a ballast mass that would bring f(ballast_mass, ballast_center_of_mass) to zero.
    f must be provided as an input and should presumably close on other characteristics
    of the boat.
    """
    from scipy import optimize
    m = optimize.fmin(f, m0, (ballast_center_of_mass,))
    return np.asscalar(m)
    
if __name__ == '__main__':
    
    length = (57*12-22-22/2.-92/2.)*inch_to_meter
    width = (6*12+10)*inch_to_meter
    ballast_cm = length-92/2*inch_to_meter
    front_depth = (9+.75)*inch_to_meter
    back_depth = 23*inch_to_meter
    desired_front_drop = 5
    desired_front_depth = desired_front_drop*inch_to_meter+front_depth
    
    b = init_boat(front_depth, back_depth, length, width)
    def f(m, cm, b=b):
        return np.abs(ballast(b, m, cm)['front_depth']-desired_front_depth)
    m_ballast = find_ballast_amount(f, 2, ballast_cm)
    b2 = ballast(b, m_ballast, ballast_cm)
    
    print 'Adding %f pounds of ballast changes front depth from %f to %f inches, back depth from %f to %f inches.'%(m_ballast*ton_to_pound, b['front_depth']*meter_to_inch, b2['front_depth']*meter_to_inch, b['back_depth']*meter_to_inch, b2['back_depth']*meter_to_inch)