"""
Determining volume conservation for SMF slump formulation.
Jane Sexton, GA 2006
"""

from math import sqrt, atan, degrees, exp, cosh, pi, radians
from Numeric import ones
from pylab import *
import scipy
from scipy.special import erf

ion()
hold(False)


x = arange(-20,25,0.1)
y = arange(-10,10,0.1)
X,Y = meshgrid(x,y)

# Grilli and Watts 2005 example
dx = 2.0
x0 = 10.0
y0 = 0.0
w = 2.0
kappa = 3
g = 9.8
d = 0.259
T = 0.052
theta = 15.0 #acos(radians(-d/T))  
lam0 = 5.0 #sqrt(g*d)
    
def func(x,y,kappad):
    
    numerator = exp(-((x-x0)/lam0)**2) - kappad*exp(-((x-dx-x0)/lam0)**2)
    denominator = cosh(kappa*(y-y0)/(w+lam0))**2
    return -numerator/denominator

eta1 = func(X,0,0.83)
eta2 = func(X,0,1)

#eta = func(X,Y,kappad)
#im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')
#plot(X/lam0,eta)
#axis([0, 5, -2, 2])

figure(1)
plot(X/lam0,eta1,'b.-',X/lam0,eta2,'g.-')
title('Surface elevation of SMF: kappa = 0.83 and 1')
xlabel('x/lambda_0')
ylabel('eta(x,y=0)')
axis([0, 5, -2, 2])

def set_integral_zero(x,dx):
    numerator = erf( (x-x0)/sqrt(lam0) )
    denominator = erf( (x-dx-x0)/sqrt(lam0) )
    return numerator/denominator

def set_integral_zero_2(x,x0):
    numerator = erf( (x-x0)/sqrt(lam0) )
    xg = d/tan(radians(theta)) + T/radians(sin(theta))
    #print 'xg', xg
    denominator = erf( (x-xg-2*x0)/sqrt(lam0) )
    return numerator/denominator

n = 100
kappad_vec1 = []
dx_vec1 = []
step = 0.001
dx = -step #+ 0.016
for i in range(n):
    dx += step
    kappad1 = float(max(set_integral_zero(x,dx)))
    dx_vec1.append(dx)
    kappad_vec1.append(kappad1)
    #print 'dx: %0.5f kappad: %f' % (dx, kappad1)

figure(2)
plot(dx_vec1,kappad_vec,'.-')
title('Volume conservation')
xlabel('delta x')
ylabel('kappa_d')
savefig('Volume_conservation')
#axis([0,0.1,0,80])

n = 60
kappad_vec2 = []
x0_vec2 = []
step = 0.1
x0 = -step
for i in range(n):
    x0 += step
    kappad2 = float(max(set_integral_zero_2(x,x0)))
    x0_vec2.append(x0)
    kappad_vec2.append(kappad2)
    print 'dx: %0.5f kappad: %f' % (x0, kappad2)


figure(3)
plot(x0_vec2,kappad_vec2,'.-')
title('Volume conservation - using xg')
xlabel('x0')
ylabel('kappa_d')
savefig('Volume_conservation_xg')