Changeset 8466 for trunk/anuga_core/source/anuga

Timestamp:

Jul 11, 2012, 9:42:10 PM (12 years ago)

Author:

davies

Message:

Adding new 'okada_tsunami' routines, which are well tested
and can treat both simple and complex faults

Location:

trunk/anuga_core/source/anuga/shallow_water

Files:

• eqf.py (modified) (1 diff)
• eqf_v2.py (modified) (6 diffs)
• okada_tsunami.py (added)
• okada_tsunami_fortran.f (added)
• okada_tsunami_octave_95.txt (added)
• test_okada_tsunami.py (added)

trunk/anuga_core/source/anuga/shallow_water/eqf.py

@@ -2 +2 @@
 # earthquake_tsunami function
 #
-
+#import okada
 """This function returns a callable object representing an initial water
    displacement generated by a submarine earthqauke.

trunk/anuga_core/source/anuga/shallow_water/eqf_v2.py

@@ -37 +37 @@
     from math import sin, radians
 
+    raise Exception, 'WARNING: (GD, 11/07/2012) there seem to be bugs in eqf_v2.py\
+    , I suggest you use okada_tsunami'
+
     if domain is not None:
         xllcorner = domain.geo_reference.get_xllcorner()
@@ -156 +159 @@
 
         for i in range(N-1):
-            self.SRECTF(alp, xr[i]*.001, yr[i]*.001, depth*.001, zero, length,\
+            # GD CHANGE
+            self.SRECTF(alp, xr[i]*.001, yr[i]*.001, depth, zero, length,\
                         zero, width, sd, cd, disl1, disl2, disl3)
+            #self.SRECTF(alp, xr[i]*.001, yr[i]*.001, depth*0.001, zero, length,\
+            #            zero, width, sd, cd, disl1, disl2, disl3)
             
             z2[i] = self.U3
@@ -281 +287 @@
 
         if Q != F0:
-            TT = atan( radians( XI*ET/(Q*R) ))
+            #TT = atan( radians( XI*ET/(Q*R) ))
+            # GD CHANGE
+            TT = atan( XI*ET/(Q*R) )
         else:
             TT = F0
@@ -299 +307 @@
         if CD == 0:
             #C==============================     
-            #C=====   INCLINED FAULT   =====     
+            #C=====   VERTICAL FAULT   =====     
             #C==============================
             RD2=RD*RD
@@ -312 +320 @@
         else:
             #C==============================   
-            #C=====   VERTICAL FAULT   =====     
+            #C=====   INCLINED FAULT   =====     
             #C==============================
             TD=SD/CD                                                          
@@ -319 +327 @@
                 A5=F0
             else:                                                    
-                A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / \
+                #A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / \
+                #                    (XI*(R+X)*CD) ))   
+                # GD CHANGE
+                A5= ALP*F2/CD*atan( ((ET*(X+Q*CD)+X*(R+X)*SD) / \
                                     (XI*(R+X)*CD) ))   
             A4= ALP/CD*(  log(RD) - SD*DLE )