Changeset 2869

Timestamp:

May 15, 2006, 5:34:11 PM (18 years ago)

Author:

sexton

Message:

updates to slump/slide formulation and document

Files:

• documentation/experimentation/smf.tex (modified) (1 diff)
• inundation/pyvolution/smf.py (modified) (5 diffs)
• production/onslow_2006/data.tex (modified) (1 diff)
• production/pt_hedland_2006/data.tex (modified) (1 diff)
• production/sydney_2006/find_max.py (modified) (4 diffs)

documentation/experimentation/smf.tex

@@ -120 +120 @@
 formulation, the values of interest shown in Figure 2 of [3] would
 be ($\kappa', \Delta x) = (1,2), (1,4), (1.2, 13.48)$ and shown in
-Figure \ref{fig:eta_vary}.
+Figure \ref{fig:eta_vary}. Note, this has not been scaled by $\eta_{\rm min}$.
 
 

inundation/pyvolution/smf.py

@@ -184 +184 @@
                   gravity=9.8, gamma=1.85, \
                   massco=1, dragco=1, frictionco=0, \
-                  dx=None, kappa=3.0, kappad=0.8, zsmall=0.01, \
+                  dx=None, kappa=3.0, kappad=1.0, zsmall=0.01, \
                   domain=None,
                   verbose=False):
@@ -248 +248 @@
         print '\t2D amplitude: ', a2D
         print '\t3D amplitude: ', a3D
+        print '\tDelta x ', dx
+        print '\tsmall ', zsmall
+        print '\tKappa d ', kappad
 
     #keep an eye on some of the assumptions built into the maths
@@ -343 +346 @@
         kappa = self.kappa
         kappad = self.kappad
+        amin = self.find_min(x0,wa,kappad,kappa)
 
         #double Gaussian calculation assumes water displacement is oriented
@@ -355 +359 @@
 
         z = zeros(N, Float)
-
+        maxz = 0.0
+        print 'here i am about to calculate z with dx = ', dx
         for i in range(N):
             try:
-                z[i] =  -am / ((cosh(kappa*(yr[i]-y0)/(wi+wa)))**2) \
+                z[i] =  -am / (amin*(cosh(kappa*(yr[i]-y0)/(wi+wa)))**2) \
                             * (exp(-((xr[i]-x0)/wa)**2) - \
                                 kappad*exp(-((xr[i]-dx-x0)/wa)**2))
+                if z[i] > maxz: maxz = z[i]
             except OverflowError:
                 pass
-
+        
+        print 'maximum of Double Gaussian function', maxz
         return z
 
@@ -386 +393 @@
 
         self.dx = 2.0 * (c * sqrt(-log((zsmall/a),e))) / 5.0
+
+    def find_min(self, x0, wa, kappad, kappa):
+        """Determine eta_min to scale eta(x,y)
+        """
+
+        from math import exp, cosh
+
+        step = 0.001
+        x = x0
+        deriv = 10.0
+        tol = 0.001
+        count_max = 100000
+        c = 0
+
+        while c < count_max and deriv > 0:
+            deriv = (x-x0)*exp(-((x-x0)/wa)**2.0) - \
+                    kappad*(x-dx-x0)*exp(-((x-dx-x0)/wa)**2.0)
+            if deriv < 0: xstar = x
+            x -= step
+            c += 1
+        
+        etastar = -am * (exp(-((xstar-x0)/wa)**2) - \
+                   kappad*exp(-((xstar-dx-x0)/wa)**2))
+
+        return etastar

production/onslow_2006/data.tex

@@ -49 +49 @@
 \begin{figure}[hbt]
 
-  \centerline{ \includegraphics[width=75mm, height=75mm]{onslow_data_extent.png}}
+  \centerline{ \includegraphics[width=100mm, height=75mm]{onslow_data_extent.png}}
 
   \caption{Data extent for Onslow scenario. Offshore data shown in blue and onshore data

production/pt_hedland_2006/data.tex

@@ -55 +55 @@
 \begin{figure}[hbt]
 
-  \centerline{ \includegraphics[width=75mm, height=75mm]{pt_hedland_data_extent.png}}
+  \centerline{ \includegraphics[width=100mm, height=75mm]{pt_hedland_data_extent.png}}
 
   \caption{Data extent for Pt Hedland scenario. Offshore data shown in blue and onshore data

production/sydney_2006/find_max.py

@@ -5 +5 @@
 
 from math import exp, cosh
+import Numeric
 
 # Grilli and Watts 2005 example
@@ -19 +20 @@
     
 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
+    numerator = exp(-((x-x0)/lam0)**2.0) - kappad*exp(-((x-dx-x0)/lam0)**2.0)
+    denominator = cosh(kappa*(y-y0)/(w+lam0))**2.0
     return -numerator/denominator
 
 def dfuncdx(x,kappad):
-    dfdx = (x-x0)**3.0*exp(-(x-x0)**2.0/lam0)-kappad*(x-dx-x0)**3.0*exp(-(x-dx-x0)**2.0/lam0)
+    dfdx = (x-x0)*exp(-((x-x0)/lam0)**2.0)-kappad*(x-dx-x0)*exp(-((x-dx-x0)/lam0)**2.0)
     return dfdx
 
@@ -34 +35 @@
 direction = 'positive'
 while c < 100000 and deriv > 0:
-    deriv = dfuncdx(x,0.83)
+    deriv = dfuncdx(x,1.0)
     if deriv < 0: xstar = x
     if direction == 'positive': x += step
@@ -40 +41 @@
     c += 1
 
-print 'location of maximum of surface elevation function: xstar = %f' % (xstar)
+print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0)
 const = 1.0 #a3D = ? #sydney = 86? and kappad=1
-print 'eta(xstar): %f' % (const*func(xstar,y0,0.83))
 
+x = arange(-20,25,0.001)
+y = arange(-10,10,0.1)
+X,Y = meshgrid(x,y)
 
-
-
+test = 0.0
+for xi in x:
+    testi = func(xi,0.0,1.0)
+    if direction == 'positive':
+        if testi > test:
+            test = testi
+            xstar = xi
+    else:
+        if testi < test:
+            test = testi
+            xstar = xi
+    
+print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)