# tipping_bucket_model
# theorize for a bit about the output of a rain gauge
# Matt Ball, October 2009
import math
import sys
import numpy as np
from pylab import *

def main():
    flip_vol, catch = getinputs()
    rain_res_sim, flips_sim, period_sim, fshape, cshape = simevents(flip_vol, catch)
    plot_sims(rain_res_sim, flips_sim, period_sim, fshape, cshape)




def getinputs():                          # properties of the bucket
    flip_vol = 4                          # typical tipping volume, mL
    catch    = math.pi*(1000/2)**2/100    # typical cathcment area, cm^2
    return flip_vol, catch




def simevents(flip_vol, catch):
    
    flip_vol_rng = np.array(range(flip_vol*0.2, flip_vol*2.0, flip_vol*0.2), dtype = float)
    catch_rng    = np.array(range(catch*0.2,   catch*2.0,    catch*0.2),     dtype = float)

    fshape = flip_vol_rng.shape   # length of flip_vol_rng, called later
    cshape = catch_rng.shape

    rain_res_sim = np.zeros([fshape[0], cshape[0]], dtype = float)   # must preallocate
    flips_sim    = np.zeros([fshape[0], cshape[0]], dtype = float)   
    period_sim   = np.zeros([fshape[0], cshape[0]], dtype = float)

    for i in range(0,fshape[0]):
        for k in range(0,cshape[0]):
            rain_res, flips, period = simrain(catch_rng[k], flip_vol_rng[i])

            rain_res_sim[i, k] = rain_res  # pull those values into separate arrays
            flips_sim[i,k] = flips
            period_sim[i,k] = period

    return rain_res_sim, flips_sim, period_sim, fshape, cshape




def simrain(catch, flip_vol):          # properties of a hypothetical storm
    rain_res = 2.54*flip_vol/catch     # smallest amount of rain that can be detected, inches

    precip   = 2.3*2.54                # cm rainfall, 2.3 is that 9/22/09 storm
    time     = 1*3600                  # seconds it took precip to fall
    in_vol   = precip*catch            # incoming volume, mL
    flips    = in_vol/flip_vol         # number of tips the bucket will see
    period   = time/flips              # seconds between bucket tips

    return rain_res, flips, period




def plot_sims(rain_res_sim, flips_sim, period_sim, fshape, cshape):
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm
    import matplotlib.pyplot as plt
    clf()
    
    fig = plt.figure()
    ax = Axes3D(fig)
    
    X = range(0,cshape[0])
    Y = range(0,fshape[0])
    X, Y = np.meshgrid(X, Y)
    #R = np.sqrt(X**2 + Y**2)
    #Z = np.sin(R)
    ax.plot_surface(X, Y, period_sim, rstride=1, cstride=1, cmap=cm.jet)
    plt.xlabel('catchment_size (cm^2)')
    plt.ylabel('flip_volume (mL)')
    plt.title('rain_resolution (in)')
    plt.show()

    
    
if __name__ == '__main__': main()