Index: /inundation/ga/storm_surge/pyvolution/smf.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/smf.py	(revision 1578)
+++ /inundation/ga/storm_surge/pyvolution/smf.py	(revision 1578)
@@ -0,0 +1,359 @@
+#
+# slide_tsunami function
+#
+
+"""This function returns a callable object representing an initial water
+   displacement generated by a submarine sediment slide.
+
+Using input parameters:
+
+Required
+ length  downslope slide length
+ depth   water depth to slide centre of mass
+ slope   bathymetric slope
+
+Optional
+ x0      x origin (0)
+ y0      y origin (0)
+ alpha   angular orientation of slide in xy plane (0)
+ w       slide width (0.25*b)
+ T       slide thickness (0.01*b)
+ g       acceleration due to gravity (9.8)
+ gamma   specific density of sediments (1.85)
+ Cm      added mass coefficient (1)
+ Cd      drag coefficient (1)
+ Cn      friction coefficient (0)
+ psi     (0)
+ dx      offset of second Gaussian (0.2*width of first Gaussian)
+ kappa   multiplier for sech^2 function (3.0)
+ kappad  multiplier for second Gaussian function (0.8)
+ zsmall  an amount near to zero (0.01)
+
+The following parameters are calculated:
+
+ a0      initial acceleration
+ ut      theoretical terminal velocity
+ s0      charactistic distance of motion
+ t0      characteristic time of motion
+ w       initial wavelength of tsunami
+ a2D     2D initial amplitude of tsunami
+ a3D     3D initial amplitude of tsunami
+
+The returned object is a callable double Gaussian function that represents
+the initial water displacement generated by a submarine sediment slide.
+
+Adrian Hitchman
+Geoscience Australia, June 2005
+"""
+
+def slide_tsunami(length, depth, slope, width=None, thickness=None, \
+                  x0=0.0, y0=0.0, alpha=0.0, \
+                  gravity=9.8, gamma=1.85, \
+                  massco=1, dragco=1, frictionco=0, psi=0, \
+                  dx=None, kappa=3.0, kappad=0.8, zsmall=0.01):
+
+    from math import sin, tan, radians, pi, sqrt, exp
+
+    #if width not provided, set to typical value
+    if width is None:
+        width = 0.25 * length
+
+    #if thickness not provided, set to typical value
+    if thickness is None:
+        thickness = 0.01 * length
+
+    #calculate some parameters of the slide
+
+    sint = sin(radians(slope))
+    tant = tan(radians(slope))
+    tanp = tan(radians(psi))
+
+    a0 = gravity * sint * ((gamma-1)/(gamma+massco)) * (1-(tanp/tant))
+    ut = sqrt((gravity*depth) * (length*sint/depth) \
+                    * (pi*(gamma-1)/(2*dragco)) * (1-(tanp/tant)))
+    s0 = ut**2 / a0
+    t0 = ut / a0
+
+    #calculate some parameters of the water displacement produced by the slide
+
+    w = t0 * sqrt(gravity*depth)
+    a2D = s0 * (0.0574 - (0.0431*sint)) \
+             * (thickness/length) \
+             * ((length*sint/depth)**1.25) \
+             * (1 - exp(-2.2*(gamma-1)))
+    a3D = a2D / (1 + (15.5*sqrt(depth/(length*sint))))
+
+    #a few temporary print statements
+    print '\nThe slide ...'
+    print '\tLength: ', length
+    print '\tDepth: ', depth
+    print '\tSlope: ', slope
+    print '\tWidth: ', width
+    print '\tThickness: ', thickness
+    print '\tx0: ', x0
+    print '\ty0: ', y0
+    print '\tAlpha: ', alpha
+    print '\tAcceleration: ', a0
+    print '\tTerminal velocity: ', ut
+    print '\tChar time: ', t0
+    print '\tChar distance: ', s0
+    print '\nThe tsunami ...'
+    print '\tWavelength: ', w
+    print '\t2D amplitude: ', a2D
+    print '\t3D amplitude: ', a3D
+
+    #keep an eye on some of the assumptions built into the maths
+
+    if ((slope < 5) or (slope > 30)):
+        print 'WARNING: slope out of range (5 - 30 degrees) ', slope
+    if ((depth/length < 0.06) or (depth/length > 1.5)):
+        print  'WARNING: d/b out of range (0.06 - 1.5) ', depth/length
+    if ((thickness/length < 0.008) or (thickness/length > 0.2)):
+        print 'WARNING: T/b out of range (0.008 - 0.2) ', thickness/length
+    if ((gamma < 1.46) or (gamma > 2.93)):
+        print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
+
+    return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
+                           x0=x0, y0=y0, alpha=alpha, \
+                           dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall)
+
+#
+# slump_tsunami function
+#
+
+"""This function returns a callable object representing an initial water
+   displacement generated by a submarine sediment slump.
+
+Using input parameters:
+
+Required
+ length  downslope slump length
+ depth   water depth to slump centre of mass
+ slope   bathymetric slope
+
+Optional
+ x0      x origin (0)
+ y0      y origin (0)
+ alpha   angular orientation of slide in xy plane (0)
+ w       slump width (1.0*b)
+ T       slump thickness (0.1*b)
+ R       slump radius of curvature (b^2/(8*T))
+ del_phi slump angular displacement (0.48)
+ g       acceleration due to gravity (9.8)
+ gamma   specific density of sediments (1.85)
+ Cm      added mass coefficient (1)
+ Cd      drag coefficient (1)
+ Cn      friction coefficient (0)
+ dx      offset of second Gaussian (0.2*width of first Gaussian)
+ kappa   multiplier for sech^2 function (3.0)
+ kappad  multiplier for second Gaussian function (0.8)
+ zsmall  an amount near to zero (0.01)
+
+The following parameters are calculated:
+
+ a0      initial acceleration
+ um      maximum velocity
+ s0      charactistic distance of motion
+ t0      characteristic time of motion
+ w       initial wavelength of tsunami
+ a2D     2D initial amplitude of tsunami
+ a3D     3D initial amplitude of tsunami
+
+The returned object is a callable double Gaussian function that represents
+the initial water displacement generated by a submarine sediment slump.
+
+Adrian Hitchman
+Geoscience Australia, June 2005
+"""
+
+def slump_tsunami(length, depth, slope, width=None, thickness=None, \
+                  radius=None, dphi=0.48, x0=0.0, y0=0.0, alpha=0.0, \
+                  gravity=9.8, gamma=1.85, \
+                  massco=1, dragco=1, frictionco=0, \
+                  dx=None, kappa=3.0, kappad=0.8, zsmall=0.01):
+
+    from math import sin, radians, sqrt
+
+    #if width not provided, set to typical value
+    if width is None:
+        width = length
+
+    #if thickness not provided, set to typical value
+    if thickness is None:
+        thickness = 0.1 * length
+
+    #if radius not provided, set to typical value
+    if radius is None:
+        radius = length**2 / (8.0 * thickness)
+
+    #calculate some parameters of the slump
+
+    sint = sin(radians(slope))
+
+    s0 = radius * dphi / 2
+    t0 = sqrt((radius*(gamma+massco)) / (gravity*(gamma-1)))
+    a0 = s0 / t0**2
+    um = s0 / t0
+
+    #calculate some parameters of the water displacement produced by the slump
+
+    w = t0 * sqrt(gravity*depth)
+    a2D = s0 * (0.131/sint) \
+             * (thickness/length) \
+             * (length*sint/depth)**1.25 \
+             * (length/radius)**0.63 * dphi**0.39 \
+             * (1.47 - (0.35*(gamma-1))) * (gamma-1)
+    a3D = a2D / (1 + (2.06*sqrt(depth/length)))
+
+    #a few temporary print statements
+    print '\nThe slump ...'
+    print '\tLength: ', length
+    print '\tDepth: ', depth
+    print '\tSlope: ', slope
+    print '\tWidth: ', width
+    print '\tThickness: ', thickness
+    print '\tRadius: ', radius
+    print '\tDphi: ', dphi
+    print '\tx0: ', x0
+    print '\ty0: ', y0
+    print '\tAlpha: ', alpha
+    print '\tAcceleration: ', a0
+    print '\tMaximum velocity: ', um
+    print '\tChar time: ', t0
+    print '\tChar distance: ', s0
+    print '\nThe tsunami ...'
+    print '\tWavelength: ', w
+    print '\t2D amplitude: ', a2D
+    print '\t3D amplitude: ', a3D
+
+    #keep an eye on some of the assumptions built into the maths
+
+    if ((slope < 10) or (slope > 30)):
+        print 'WARNING: slope out of range (10 - 30 degrees) ', slope
+    if ((depth/length < 0.34) or (depth/length > 0.5)):
+        print  'WARNING: d/b out of range (0.34 - 0.5) ', depth/length
+    if ((thickness/length < 0.10) or (thickness/length > 0.15)):
+        print 'WARNING: T/b out of range (0.10 - 0.15) ', thickness/length
+    if ((radius/length < 1.0) or (radius/length > 2.0)):
+        print 'WARNING: R/b out of range (1 - 2) ', radius/length
+    if ((dphi < 0.10) or (dphi > 0.52)):
+        print 'WARNING: del_phi out of range (0.10 - 0.52) ', dphi
+    if ((gamma < 1.46) or (gamma > 2.93)):
+        print 'WARNING: gamma out of range (1.46 - 2.93) ', gamma
+
+    return Double_gaussian(a3D=a3D, wavelength=w, width=width, \
+                           x0=x0, y0=y0, alpha=alpha, \
+                           dx=dx, kappa=kappa, kappad=kappad, zsmall=zsmall)
+
+#
+# Double_gaussian class
+#
+
+"""This is a callable class representing the initial water displacment 
+   generated by a submarine mass failure (smf -> sediment slide or slump).
+
+Using input parameters:
+
+Required
+ w       initial wavelength of tsunami
+ a3D     3D initial amplitude of tsunami
+ width   width of smf
+
+Optional
+ x0      x origin of smf
+ y0      y origin of smf
+ alpha   angular orientation of smf in xy plane (0)
+ dx      offset of second Gaussian (0.2*width of first Gaussian)
+ kappa   multiplier for sech^2 function (3.0)
+ kappad  multiplier for second Gaussian function (0.8)
+ zsmall  an amount near to zero (0.01)
+
+Adrian Hitchman
+Geoscience Australia, June 2005
+"""
+
+class Double_gaussian:
+
+    def __init__(self, a3D, wavelength, width, x0, y0, alpha, \
+                 dx, kappa, kappad, zsmall):
+        self.a3D = a3D
+        self.wavelength = wavelength
+        self.width = width
+        self.x0 = x0
+        self.y0 = y0
+        self.alpha = alpha
+        self.kappa = kappa
+        self.kappad = kappad
+
+        if dx is None:
+            self.determineDX(zsmall=zsmall)
+        else:
+            self.dx = dx
+
+    def __call__(self, x, y):
+        """Make Double_gaussian a callable object.
+
+        If called as a function, this object returns z values representing
+        the initial 3D distribution of water heights at the points (x,y)
+        produced by a submarine mass failure.
+        """
+
+        from math import sin, cos, radians, exp, cosh
+        from Numeric import zeros, Float
+
+        #ensure vectors x and y have the same length
+        N = len(x)
+        assert N == len(y)
+
+        am = self.a3D
+        wa = self.wavelength
+        wi = self.width
+        x0 = self.x0
+        y0 = self.y0
+        alpha = self.alpha
+        dx = self.dx
+        kappa = self.kappa
+        kappad = self.kappad
+
+        #double Gaussian calculation assumes water displacement is oriented
+        #E-W, so, for displacement at some angle alpha clockwise from the E-W
+        #direction, rotate (x,y) coordinates anti-clockwise by alpha
+
+        cosa = cos(radians(alpha))
+        sina = sin(radians(alpha))
+
+        xr = ((x-x0) * cosa - (y-y0) * sina) + x0
+        yr = ((x-x0) * sina + (y-y0) * cosa) + y0
+
+        z = zeros(N, Float)
+
+        for i in range(N):
+            try:
+                z[i] =  -am / ((cosh(kappa*(yr[i]-y0)/(wi+wa)))**2) \
+                            * (exp(-((xr[i]-x0)/wa)**2) - \
+                                kappad*exp(-((xr[i]-dx-x0)/wa)**2))
+            except OverflowError:
+                pass
+
+        return z
+
+    def determineDX(self, zsmall):
+        """Determine a suitable offset for the second Gaussian function.
+
+        A suitable offset for the second Gaussian function is taken to
+        be some fraction of the 'width' of the first Gaussian function.
+
+        The 'width' of the first Gaussian is obtained from the range of
+        the x coordinates over which the function takes values from
+        'near zero', through 1, and back to 'near zero'.
+
+        The parameter zsmall passed to this function specifies how much
+        'near zero' is.
+        """
+
+        from math import sqrt, log, e
+
+        a = self.a3D
+        c = self.wavelength
+
+        self.dx = 2.0 * (c * sqrt(-log((zsmall/a),e))) / 5.0
