kinema.py

#!/usr/bin/env python3

from math import acosh, atan, cos, cosh, log, sqrt, tan, tanh

viscosity = {'air': 1.9e-5,
     'water': 8.9e-4,
     'olive oil': 0.081}

density = {'air': 1.275,
     'olive oil': 860.,
     'water': 1000.,
     'granite': 2600,
     'iron': 7874.,
     'lead': 11340}

def reynolds(ρ, l, v, η):
    '''Parameters:
    ρ: density (kg/m^3)
    l: cross-sectional length (m)
    v: velocity (m/s)
    η: dynamic viscosity (Pa s)
    
    Returns:
    Reynolds number (dimensionless)'''
    
    return (ρ*l*v) / η

def quadratic_drag(ρ, η, ρ_ball, R, v_0, gravity=9.80665):
    '''Parameters:
    ρ: density of medium (kg/m^3)
    η: dynamic viscosity of medium (Pa s)
    ρ_ball: density of object (kg/m^3)
    R: radius of object (m)
    v_0: initial velocity of object (m/s)
    gravity: gravity (m/s)
    
    Returns:
    v: velocity as a function of time (m/s)
    h: height as a function of time (m)
    t_up: time to reach maximum height (s)
    t_down: time to reach ground from maximum height (s)
    v_term: terminal velocity (m/s)
    v_hit: hit velocity (m/s)'''
    
    g = (1 - ρ/ρ_ball) * gravity
    γ = sqrt(15) * sqrt(ρ/(ρ_ball*g*R)) / 10
    
    v_up = lambda t: 1/γ * tan(-γ*g*t + atan(γ*v_0))
    t_up = 1/(γ*g) * atan(γ*v_0)
    h_max = 1/(γ**2*g) * log(sqrt(1 + (γ*v_0)**2))
    h_up = lambda t: 1/(γ**2*g) * log(sqrt(1 + (γ*v_0)**2) * cos(-γ*g*t + atan(γ*v_0)))
    
    v_down = lambda t: -1/γ * tanh(γ*g*(t-t_up))
    v_term = -1/γ
    h_down = lambda t: 1/(γ**2*g) * log(sqrt(1+(γ*v_0)**2) / cosh(γ*g*(t-t_up)))
    t_down = 1/(γ*g) * acosh(sqrt(1 + (γ*v_0)**2))
    
    v_hit = v_down(t_up + t_down) # = -v_0/sqrt(1+(γ*v_0)**2)
    
    v = lambda t: v_up(t) if t <= t_up else v_down(t)
    h = lambda t: h_up(t) if t <= t_up else h_down(t)
    
    return (v, h, t_up, t_down, v_term, v_hit)

# Sources and attribution:
# http://blog.kjempekjekt.com/2014/03/24/spillkata-2-kinema/
# http://www.math.upatras.gr/~weele/files/What%20goes%20up%20must%20come%20down%20(paper%2027).pdf