#!usr/bin/env python
"""This script shows the progress of the height of the water level 
of a cylindrical bucket with a small hole at the bottom.

Soil Physics Theroy - 6583
Oklahoma State University

Andres Patrignani - October 2015"""

# Import modules
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import mpld3
mpld3.enable_notebook()


# Define parameters and constants
Q_in           = 7                          # cm^3/s
nbuckets       = 1                          # number of buckets
radius_bucket  = 2.5                        # cm
area_bucket    = np.pi*radius_bucket**2     # cm^2
height_bucket  = 30                         # cm
radius_hole    = 0.1                        # cm
area_hole      = np.pi*radius_hole**2       # cm^2
g              = 981                        # cm/s^2
density_water  = 1                          # [g/cm^3]   rho
dyn_viscosity  = 0.01                       # g/(cm s)
dt             = 1                          # s
total_time     = 250                        # s
time_vector    = np.arange(0,total_time,dt) # s

# Pre-allocate variables
L = len(time_vector)
height = np.ones((L,nbuckets))

# Reynolds number (Re) and coefficient of discharge (CD) lookup table
Re_table = np.log10([1, 10, 100, 1000, 10000, 100000]); # dimensionless
CD_table = [0.04, 0.28, 0.73, 0.74, 0.64, 0.61];        # dimensionless

### EXPLICIT METHOD (Forward Euler method) ###

# Initial conditions (t=0)
height[0,:] = 0 # Bucket is empty at t=0

# Subsequent time steps (t=1 to end)
for i in range(1, L):
    # Calculate flux (q_in) into the bucket for current time step
    # Inflow (Q_in) is constant, so this step is easy
    q_in = Q_in / area_bucket # volume per unit area per unit time
    
    # Calculate flux out (q_out) of the bucket for current time step
    # Outflow (Q_out) is dependent on the current water into the bucket
    # so we need to do some extra calculations to find q_out.
    velocity = sqrt(2 * g * height[i-1,0])
    Re = 2 * radius_hole * density_water / dyn_viscosity * velocity
    Re = max(Re, 1)
    CD = np.interp(log10(Re), Re_table, CD_table)
    Q_out = CD * area_hole * velocity
    q_out = Q_out / area_bucket # volume per unit area per unit time
    
    # Calculate height of water into the bucket (mass balance)
    height[i,0] = min(height[i-1,0] + dt*(q_in - q_out), height_bucket)
    
    
# Plot results (Explore the chart using the toolbar at the bottom left corner of the plot)
plt.title('Explicit method',fontsize=24)
plt.plot(time_vector,height,'-')

plt.xlabel('Time (s)',fontsize=20)
plt.xticks(fontsize=16)

plt.ylabel('Height (cm)',fontsize=20)
plt.yticks(fontsize=16)
plt.ylim([0,height_bucket+5])
plt.grid(color='lightgray')

plt.show()