Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 4319)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 4320)
@@ -1060,4 +1060,6 @@
         comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv'
         fid_compare = open(comparefile, 'w')
+        file0 = file_loc[j]+'gauges_t0.csv'
+        fid_0 = open(file0, 'w')
         ##### loop over each gauge #####
         for k in gauge_index:
@@ -1102,4 +1104,7 @@
                     s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel, z, uh, vh)
                     fid_out.write(s)
+                    if t == 0:
+                        s = '%.2f, %.2f, %.2f\n' %(g[0], g[1], w)
+                        fid_0.write(s)
                     if t/60.0 <= 13920: tindex = i
                     if w > max_stage: max_stage = w
@@ -1157,4 +1162,5 @@
         depth_axis = axis([time_min/60.0, time_max/60.0, -0.1, max(max_depths)*1.1])
         stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
+        #stage_axis = axis([50.0, time_max/60.0, -3.0, 2.0])
         vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
         mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1])
@@ -1195,5 +1201,5 @@
                         axis(depth_axis)
                     if which_quantity == 'stage':
-                        if elevations[0,k,j] < 0:
+                        if elevations[0,k,j] <= 0:
                             plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j])
                             axis(stage_axis)
@@ -1273,4 +1279,5 @@
                     if title_on == True:
                         title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc2))
+                        #title('Gauge %s (MOST elevation %.2f, ANUGA elevation %.2f)' %(gaugeloc2, elevations[10,k,0], elevations[10,k,1] ))
 
                     savefig(graphname) # save figures with sww file
Index: /anuga_core/source/anuga/shallow_water/eqf.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/eqf.py	(revision 4319)
+++ /anuga_core/source/anuga/shallow_water/eqf.py	(revision 4320)
@@ -26,9 +26,9 @@
 """
 
-def earthquake_tsunami(length, width, strike, dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\
-                  domain=None,
-                  verbose=False):
+def earthquake_tsunami(length, width, strike, depth, \
+                       dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\
+                       domain=None, verbose=False):
 
-    from math import sin, radians, sqrt
+    from math import sin, radians
 
     if domain is not None:
@@ -37,5 +37,5 @@
         x0 = x0 - xllcorner  # fault origin (relative)
         y0 = y0 - yllcorner
-
+    
     #a few temporary print statements
     if verbose is True:
@@ -50,8 +50,15 @@
         print '\ty0: ', y0
 
+    # warning test
+    test = width*1000.0*sin(radians(dip)) - depth
 
-    return Okada(length=length, width=width, dip=dip, \
-                 x0=x0, y0=y0, strike=strike, depth=depth, \
-                 slip=slip, rake=rake)
+    if verbose is True:
+        if test > 0.0:
+            print 'Earthquake source not located below seafloor'
+            print 'Please check depth'
+
+    return Okada_func(length=length, width=width, dip=dip, \
+                      x0=x0, y0=y0, strike=strike, depth=depth, \
+                      slip=slip, rake=rake, test=test)
 
 
@@ -80,8 +87,8 @@
 """
 
-class Okada:
+class Okada_func:
 
     def __init__(self, length, width, dip, x0, y0, strike, \
-                 depth, slip, rake):
+                 depth, slip, rake, test):
         self.dip = dip
         self.length = length
@@ -93,8 +100,9 @@
         self.slip = slip
         self.rake = rake
+        self.test = test
 
 
     def __call__(self, x, y):
-        """Make Okada a callable object.
+        """Make Okada_func a callable object.
 
         If called as a function, this object returns z values representing
@@ -105,5 +113,5 @@
         from math import sin, cos, radians, exp, cosh
         from Numeric import zeros, Float
-        from okada import okada
+        from okada import okadatest
 
         #ensure vectors x and y have the same length
Index: /anuga_core/source/anuga/shallow_water/eqf_v2.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/eqf_v2.py	(revision 4320)
+++ /anuga_core/source/anuga/shallow_water/eqf_v2.py	(revision 4320)
@@ -0,0 +1,540 @@
+#
+# earthquake_tsunami function
+#
+
+"""This function returns a callable object representing an initial water
+   displacement generated by a submarine earthqauke.
+   
+Using input parameters:
+
+Required
+ length  along-stike length of rupture area
+ width   down-dip width of rupture area
+ strike  azimuth (degrees, measured from north) of fault axis
+ dip     angle of fault dip in degrees w.r.t. horizontal
+ depth   depth to base of rupture area
+ 
+Optional
+ x0      x origin (0)
+ y0      y origin (0)
+ slip    metres of fault slip (1)
+ rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
+
+The returned object is a callable okada function that represents
+the initial water displacement generated by a submarine earthuake.
+
+"""
+
+def earthquake_tsunami(length, width, strike, depth, \
+                       dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\
+                       domain=None, verbose=False):
+
+    from math import sin, radians
+
+    if domain is not None:
+        xllcorner = domain.geo_reference.get_xllcorner()
+        yllcorner = domain.geo_reference.get_yllcorner()
+        x0 = x0 - xllcorner  # fault origin (relative)
+        y0 = y0 - yllcorner
+
+    #a few temporary print statements
+    if verbose is True:
+        print '\nThe Earthquake ...'
+        print '\tLength: ', length
+        print '\tDepth: ', depth
+        print '\tStrike: ', strike
+        print '\tWidth: ', width
+        print '\tDip: ', dip
+        print '\tSlip: ', slip
+        print '\tx0: ', x0
+        print '\ty0: ', y0
+
+    # warning state
+    test = width*1000.0*sin(radians(dip)) - depth
+    
+    if test > 0.0:
+        msg = 'Earthquake source not located below seafloor - check depth'
+        raise Exception, msg
+
+    return Okada_func(length=length, width=width, dip=dip, \
+                      x0=x0, y0=y0, strike=strike, depth=depth, \
+                      slip=slip, rake=rake)
+
+#
+# Okada class
+#
+
+"""This is a callable class representing the initial water displacment 
+   generated by an earthquake.
+
+Using input parameters:
+
+Required
+ length  along-stike length of rupture area
+ width   down-dip width of rupture area
+ strike  azimuth (degrees, measured from north) of fault axis
+ dip     angle of fault dip in degrees w.r.t. horizontal
+ depth   depth to base of rupture area
+ 
+Optional
+ x0      x origin (0)
+ y0      y origin (0)
+ slip    metres of fault slip (1)
+ rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
+
+"""
+
+class Okada_func:
+
+    def __init__(self, length, width, dip, x0, y0, strike, \
+                 depth, slip, rake):
+        self.dip = dip
+        self.length = length
+        self.width = width
+        self.x0 = x0
+        self.y0 = y0
+        self.strike = strike
+        self.depth = depth
+        self.slip = slip
+        self.rake = rake
+
+
+    def __call__(self, x, y):
+        """Make Okada_func 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
+        #from okada import okadatest
+
+        #ensure vectors x and y have the same length
+        N = len(x)
+        assert N == len(y)
+
+        depth = self.depth
+        dip = self.dip
+        length = self.length
+        width = self.width
+        x0 = self.x0
+        y0 = self.y0
+        strike = self.strike
+        dip = self.dip
+        rake = self.rake
+        slip = self.slip
+        
+        #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(strike))
+        sina = sin(radians(strike))
+
+        xr = ( (x-x0) * sina + (y-y0) * cosa)
+        yr = (-(x-x0) * cosa + (y-y0) * sina)
+
+        # works on nautilus when have already done
+        # f2py -c okada.f -m okada
+        #z1 = okada(xr,yr,depth,length,width,dip,rake,slip)
+
+        z2 = zeros(N, Float)
+        alp = 0.5
+        disl3 = 0.0
+        zero = 0.0
+        disl1 = slip*cos(radians(rake))
+        disl2 = slip*sin(radians(rake))
+        cd = cos(radians(dip))
+        sd = sin(radians(dip))
+
+        for i in range(N-1):
+            self.SRECTF(alp,xr[i]*.001,yr[i]*.001,depth*.001,zero,length,\
+                        zero,width,sd,cd,disl1,disl2,disl3)
+            
+            z2[i] = self.U3
+
+        return z2
+
+    def SRECTF(self,ALP,X,Y,DEP,AL1,AL2,AW1,AW2,SD,CD,DISL1,DISL2,DISL3):
+        """ SURFACE DISPLACEMENT,STRAIN,TILT DUE TO RECTANGULAR FAULT
+            IN A HALF-SPACE. CODED BY  Y.OKADA ... JAN 1985                
+             
+                                                                          
+         INPUT                                                            
+           ALP     : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)   
+           X,Y     : COORDINATE OF STATION                                
+           DEP     : SOURCE DEPTH                                         
+           AL1,AL2 : FAULT LENGTH RANGE                                   
+           AW1,AW2 : FAULT WIDTH RANGE                                    
+           SD,CD   : SIN,COS OF DIP-ANGLE                                 
+                  (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)
+           DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION      
+                                                                              
+         OUTPUT                                                           
+           U1, U2, U3      : DISPLACEMENT ( UNIT= UNIT OF DISL     )      
+           U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF DISL /          
+           U31,U32         : TILT                 UNIT OF X,Y,,,AW )      
+                                                                              
+         SUBROUTINE USED...SRECTG                                         
+        """                                                                    
+        
+        from Numeric import zeros, Float
+        U = zeros(9, Float)
+        DU = zeros(9, Float)
+        
+        F0 = 0.0
+        F1 = 1.0
+                                                                    
+        P = Y*CD + DEP*SD                                                 
+        Q = Y*SD - DEP*CD
+
+        KVEC = [1,2]
+        JVEC = [1,2]
+        for K in KVEC:
+            if K == 1: ET=P-AW1
+            if K == 2: ET=P-AW2
+            for J in JVEC:
+                if J == 1: XI=X-AL1
+                if J == 2: XI=X-AL2
+                JK=J+K
+                if JK <> 3:
+                    SIGN= F1
+                else:
+                    SIGN=-F1
+
+                self.SRECTG(ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3)
+            
+                DU[0] = self.DU1
+                DU[1] = self.DU2
+                DU[2] = self.DU3
+                DU[3] = self.DU11
+                DU[4] = self.DU12
+                DU[5] = self.DU21
+                DU[6] = self.DU22
+                DU[7] = self.DU31
+                DU[8] = self.DU32
+            
+                for i in range(len(U)):
+                    U[i]=U[i]+SIGN*DU[i]
+            
+        U1 =U[0]                                                          
+        U2 =U[1]                                                          
+        U3 =U[2]                                                          
+        U11=U[3]                                                          
+        U12=U[4]                                                          
+        U21=U[5]                                                          
+        U22=U[6]                                                          
+        U31=U[7]                                                          
+        U32=U[8]
+
+        self.U3 = U3
+        
+    def SRECTG(self,ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3):
+        """
+            C                                                                       
+            C*********************************************************************  
+            C*****                                                           *****  
+            C*****  INDEFINITE INTEGRAL OF SURFACE DISPLACEMENT,STRAIN,TILT  *****  
+            C*****  DUE TO RECTANGULAR FAULT IN A HALF-SPACE                 *****  
+            C*****                          CODED BY  Y.OKADA ... JAN 1985   *****  
+            C*****                                                           *****  
+            C*********************************************************************  
+            C                                                                       
+            C***** INPUT                                                            
+            C*****   ALP     : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)   
+            C*****   XI,ET,Q : FAULT COORDINATE                                     
+            C*****   SD,CD   : SIN,COS OF DIP-ANGLE                                 
+            C*****          (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)
+            C*****   DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION      
+            C                                                                       
+            C***** OUTPUT                                                           
+            C*****   U1, U2, U3      : DISPLACEMENT ( UNIT= UNIT OF DISL    )       
+            C*****   U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF DISL /          
+            C*****   U31,U32         : TILT                 UNIT OF XI,ET,Q )       
+            C
+        """
+        
+        from math import sqrt, atan, log, radians
+        F0 = 0.0
+        F1 = 1.0
+        F2 = 2.0
+        PI2=6.283185307179586
+
+        XI2=XI*XI                                                         
+        ET2=ET*ET                                                         
+        Q2=Q*Q                                                            
+        R2=XI2+ET2+Q2                                                     
+        R =sqrt(R2)                                                      
+        R3=R*R2                                                           
+        D =ET*SD-Q*CD                                                     
+        Y =ET*CD+Q*SD                                                     
+        RET=R+ET
+        if RET < F0: RET=F0
+        RD =R+D                                                           
+        RRD=F1/(R*RD)                                                     
+
+        if Q <> F0:
+            TT = atan( radians( XI*ET/(Q*R) ))
+        else:
+            TT = F0
+
+        if RET <> F0:
+            RE = F1/RET
+            DLE= log(RET)
+        else:
+            RE = F0
+            DLE=-log(R-ET)
+                           
+        RRX=F1/(R*(R+XI))                                                 
+        RRE=RE/R                                                          
+        AXI=(F2*R+XI)*RRX*RRX/R                                           
+        AET=(F2*R+ET)*RRE*RRE/R
+
+        if CD == 0:
+            #C==============================                                         
+            #C=====   INCLINED FAULT   =====                                         
+            #C==============================
+            RD2=RD*RD
+            A1=-ALP/F2*XI*Q/RD2                                               
+            A3= ALP/F2*( ET/RD + Y*Q/RD2 - DLE )                              
+            A4=-ALP*Q/RD                                                      
+            A5=-ALP*XI*SD/RD                                                  
+            B1= ALP/F2*  Q  /RD2*(F2*XI2*RRD - F1)                            
+            B2= ALP/F2*XI*SD/RD2*(F2*Q2 *RRD - F1)                            
+            C1= ALP*XI*Q*RRD/RD                                               
+            C3= ALP*SD/RD*(XI2*RRD - F1)  
+        else:
+            #C==============================                                         
+            #C=====   VERTICAL FAULT   =====                                         
+            #C==============================
+            TD=SD/CD                                                          
+            X =sqrt(XI2+Q2)
+            if XI == F0:
+                A5=F0
+            else:                                                    
+                A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / (XI*(R+X)*CD) ))   
+            A4= ALP/CD*(  log(RD) - SD*DLE )                                  
+            A3= ALP*(Y/RD/CD - DLE) + TD*A4                                   
+            A1=-ALP/CD*XI/RD        - TD*A5                                   
+            C1= ALP/CD*XI*(RRD - SD*RRE)                                      
+            C3= ALP/CD*(Q*RRE - Y*RRD)                                        
+            B1= ALP/CD*(XI2*RRD - F1)/RD - TD*C3                              
+            B2= ALP/CD*XI*Y*RRD/RD       - TD*C1
+                                        
+        A2=-ALP*DLE - A3                                                  
+        B3=-ALP*XI*RRE - B2                                               
+        B4=-ALP*( CD/R + Q*SD*RRE ) - B1                                  
+        C2= ALP*(-SD/R + Q*CD*RRE ) - C3                                  
+                                                                 
+        U1 =F0                                                            
+        U2 =F0                                                            
+        U3 =F0                                                            
+        U11=F0                                                            
+        U12=F0                                                            
+        U21=F0                                                            
+        U22=F0                                                            
+        U31=F0                                                            
+        U32=F0
+        
+        if DISL1 <> F0:
+            #C======================================                                 
+            #C=====  STRIKE-SLIP CONTRIBUTION  =====                                 
+            #C======================================
+            UN=DISL1/PI2                                                      
+            REQ=RRE*Q                                                         
+            U1 =U1 - UN*( REQ*XI +   TT    + A1*SD )                          
+            U2 =U2 - UN*( REQ*Y  + Q*CD*RE + A2*SD )                          
+            U3 =U3 - UN*( REQ*D  + Q*SD*RE + A4*SD )                          
+            U11=U11+ UN*( XI2*Q*AET - B1*SD )                                 
+            U12=U12+ UN*( XI2*XI*( D/(ET2+Q2)/R3 - AET*SD ) - B2*SD )         
+            U21=U21+ UN*( XI*Q/R3*CD + (XI*Q2*AET - B2)*SD )                  
+            U22=U22+ UN*( Y *Q/R3*CD + (Q*SD*(Q2*AET-F2*RRE) -(XI2+ET2)/R3*CD - B4)*SD )           
+            U31=U31+ UN*(-XI*Q2*AET*CD + (XI*Q/R3 - C1)*SD )                  
+            U32=U32+ UN*( D*Q/R3*CD + (XI2*Q*AET*CD - SD/R + Y*Q/R3 - C2)*SD )
+
+        if DISL2 <> F0:
+            #C===================================                                    
+            #C=====  DIP-SLIP CONTRIBUTION  =====                                    
+            #C=================================== 
+            UN=DISL2/PI2                                                      
+            SDCD=SD*CD                                                        
+            U1 =U1 - UN*( Q/R             - A3*SDCD )                         
+            U2 =U2 - UN*( Y*Q*RRX + CD*TT - A1*SDCD )                         
+            U3 =U3 - UN*( D*Q*RRX + SD*TT - A5*SDCD )                         
+            U11=U11+ UN*( XI*Q/R3            + B3*SDCD )                      
+            U12=U12+ UN*( Y *Q/R3 - SD/R     + B1*SDCD )                      
+            U21=U21+ UN*( Y *Q/R3 + Q*CD*RRE + B1*SDCD )                      
+            U22=U22+ UN*( Y*Y*Q*AXI - (F2*Y*RRX + XI*CD*RRE)*SD + B2*SDCD )   
+            U31=U31+ UN*( D *Q/R3 + Q*SD*RRE + C3*SDCD )                      
+            U32=U32+ UN*( Y*D*Q*AXI - (F2*D*RRX + XI*SD*RRE)*SD + C1*SDCD ) 
+            
+        if DISL3 <> F0:
+            #C========================================
+            #C=====  TENSILE-FAULT CONTRIBUTION  =====                               
+            #C======================================== 
+            UN=DISL3/PI2                                                      
+            SDSD=SD*SD                                                        
+            U1 =U1 + UN*( Q2*RRE                       - A3*SDSD )            
+            U2 =U2 + UN*(-D*Q*RRX - SD*(XI*Q*RRE - TT) - A1*SDSD )            
+            U3 =U3 + UN*( Y*Q*RRX + CD*(XI*Q*RRE - TT) - A5*SDSD )            
+            U11=U11- UN*( XI*Q2*AET             + B3*SDSD )                   
+            U12=U12- UN*(-D*Q/R3 - XI2*Q*AET*SD + B1*SDSD )                   
+            U21=U21- UN*( Q2*(CD/R3 + Q*AET*SD) + B1*SDSD )                   
+            U22=U22- UN*((Y*CD-D*SD)*Q2*AXI - F2*Q*SD*CD*RRX - (XI*Q2*AET - B2)*SDSD )                   
+            U31=U31- UN*( Q2*(SD/R3 - Q*AET*CD) + C3*SDSD )                   
+            U32=U32- UN*((Y*SD+D*CD)*Q2*AXI + XI*Q2*AET*SD*CD - (F2*Q*RRX - C1)*SDSD )
+
+        self.DU1 = U1
+        self.DU2 = U2
+        self.DU3 = U3
+
+        self.DU11 = U11
+        self.DU12 = U12
+
+        self.DU21 = U21
+        self.DU22 = U22
+        
+        self.DU31 = U31
+        self.DU32 = U32
+        
+    def spoint(self,alp,x,y,dep,sd,cd,pot1,pot2,pot3):
+        """ Calculate surface displacement, strain, tilt due to buried point
+        source in a semiinfinite medium. Y. Okada Jan 1985
+
+        Input:
+                                                            
+        ALP   : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)     
+        X,Y   : COORDINATE OF STATION                                  
+        DEP     : SOURCE DEPTH                                         
+        SD,CD : SIN,COS OF DIP-ANGLE                                   
+               (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)
+        POT1,POT2,POT3 : STRIKE-, DIP- AND TENSILE-POTENCY             
+            POTENCY=(  MOMENT OF DOUBLE-COUPLE  )/MYU    FOR POT1,2    
+            POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMDA  FOR POT3
+            
+        Output:
+
+        U1, U2, U3      : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) /      
+                        :                     (UNIT OF X,Y,D)**2  )    
+        U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF POTENCY) /      
+        U31,U32         : TILT                (UNIT OF X,Y,D)**3  )
+        """
+
+        from math import sqrt
+
+        F0 = 0.0
+        F1 = 1.0
+        F2 = 2.0
+        F3 = 3.0
+        F4 = 4.0
+        F5 = 5.0
+        F8 = 8.0
+        F9 = 9.0
+        PI2 = 6.283185307179586
+
+        D =DEP                                                            
+        P =Y*CD + D*SD                                                    
+        Q =Y*SD - D*CD                                                    
+        S =P*SD + Q*CD                                                    
+        X2=X*X                                                            
+        Y2=Y*Y                                                            
+        XY=X*Y                                                            
+        D2=D*D                                                            
+        R2=X2 + Y2 + D2                                                   
+        R =sqrt(R2)                                                       
+        R3=R *R2
+        R5=R3*R
+        QR=F3*Q/R
+        XR =F5*X2/R2
+        YR =F5*Y2/R2
+        XYR=F5*XY/R2
+        DR =F5*D /R2
+        RD =R + D
+        R12=F1/(R*RD*RD)
+        R32=R12*(F2*R + D)/ R2
+        R33=R12*(F3*R + D)/(R2*RD)
+        R53=R12*(F8*R2 + F9*R*D + F3*D2)/(R2*R2*RD)
+        R54=R12*(F5*R2 + F4*R*D +    D2)/R3*R12                           
+
+        A1= ALP*Y*(R12-X2*R33)                                            
+        A2= ALP*X*(R12-Y2*R33)                                            
+        A3= ALP*X/R3 - A2                                                 
+        A4=-ALP*XY*R32                                                    
+        A5= ALP*( F1/(R*RD) - X2*R32 )                                    
+        B1= ALP*(-F3*XY*R33      + F3*X2*XY*R54)                          
+        B2= ALP*( F1/R3 - F3*R12 + F3*X2*Y2*R54)                          
+        B3= ALP*( F1/R3 - F3*X2/R5) - B2                                  
+        B4=-ALP*F3*XY/R5 - B1                                             
+        C1=-ALP*Y*(R32 - X2*R53)                                          
+        C2=-ALP*X*(R32 - Y2*R53)                                          
+        C3=-ALP*F3*X*D/R5 - C2                                            
+
+        U1 =F0                                                            
+        U2 =F0                                                            
+        U3 =F0                                                            
+        U11=F0                                                            
+        U12=F0                                                            
+        U21=F0                                                            
+        U22=F0                                                            
+        U31=F0                                                            
+        U32=F0
+        
+        #======================================                                 
+        #=====  STRIKE-SLIP CONTRIBUTION  =====                                 
+        #======================================                                 
+
+        if POT1 <> F0:
+            UN=POT1/PI2                                                       
+            QRX=QR*X                                                          
+            FX=F3*X/R5*SD                                                     
+            U1 =U1 - UN*( QRX*X + A1*SD )                                     
+            U2 =U2 - UN*( QRX*Y + A2*SD )                                     
+            U3 =U3 - UN*( QRX*D + A4*SD )
+              
+            U11=U11- UN*( QRX* (F2-XR)        + B1*SD )                       
+            U12=U12- UN*(-QRX*XYR      + FX*X + B2*SD )                       
+            U21=U21- UN*( QR*Y*(F1-XR)        + B2*SD )                       
+            U22=U22- UN*( QRX *(F1-YR) + FX*Y + B4*SD )                       
+            U31=U31- UN*( QR*D*(F1-XR)        + C1*SD )                       
+            U32=U32- UN*(-QRX*DR*Y     + FX*D + C2*SD )                       
+
+        #======================================                                 
+        #=====    DIP-SLIP CONTRIBUTION   =====                                 
+        #======================================                                 
+
+        if POT2 <> F0:
+            UN=POT2/PI2                                                       
+            SDCD=SD*CD                                                        
+            QRP=QR*P                                                          
+            FS=F3*S/R5                                                        
+            U1 =U1 - UN*( QRP*X - A3*SDCD )                                   
+            U2 =U2 - UN*( QRP*Y - A1*SDCD )                                   
+            U3 =U3 - UN*( QRP*D - A5*SDCD )                                   
+            U11=U11- UN*( QRP*(F1-XR)        - B3*SDCD )                      
+            U12=U12- UN*(-QRP*XYR     + FS*X - B1*SDCD )                      
+            U21=U21- UN*(-QRP*XYR            - B1*SDCD )                      
+            U22=U22- UN*( QRP*(F1-YR) + FS*Y - B2*SDCD )                      
+            U31=U31- UN*(-QRP*DR*X           - C3*SDCD )                      
+            U32=U32- UN*(-QRP*DR*Y    + FS*D - C1*SDCD )
+            
+        #========================================                               
+        #=====  TENSILE-FAULT CONTRIBUTION  =====                               
+        #========================================                                                                    
+
+        if POT3 <> F0:
+            UN=POT3/PI2                                                       
+            SDSD=SD*SD                                                        
+            QRQ=QR*Q                                                          
+            FQ=F2*QR*SD                                                       
+            U1 =U1 + UN*( QRQ*X - A3*SDSD )                                   
+            U2 =U2 + UN*( QRQ*Y - A1*SDSD )                                   
+            U3 =U3 + UN*( QRQ*D - A5*SDSD )                                   
+            U11=U11+ UN*( QRQ*(F1-XR)        - B3*SDSD )                      
+            U12=U12+ UN*(-QRQ*XYR     + FQ*X - B1*SDSD )                      
+            U21=U21+ UN*(-QRQ*XYR            - B1*SDSD )                      
+            U22=U22+ UN*( QRQ*(F1-YR) + FQ*Y - B2*SDSD )                      
+            U31=U31+ UN*(-QRQ*DR*X           - C3*SDSD )                      
+            U32=U32+ UN*(-QRQ*DR*Y    + FQ*D - C1*SDSD )
+
Index: /anuga_core/source/anuga/shallow_water/okada.f
===================================================================
--- /anuga_core/source/anuga/shallow_water/okada.f	(revision 4319)
+++ /anuga_core/source/anuga/shallow_water/okada.f	(revision 4320)
@@ -1,3 +1,3 @@
-      subroutine okada(x,y,z,n,depth,xlength,xwidth,dip,rake,slip)
+      subroutine okadatest(x,y,z,n,depth,xlength,xwidth,dip,rake,slip)
 Cf2py   intent(out) z
 Cf2py   intent(in) x,y,n,depth,dip,rake,slip,xlength,xwidth
Index: /anuga_work/production/pt_hedland_2006/run_building_inundation.py
===================================================================
--- /anuga_work/production/pt_hedland_2006/run_building_inundation.py	(revision 4319)
+++ /anuga_work/production/pt_hedland_2006/run_building_inundation.py	(revision 4320)
@@ -34,4 +34,5 @@
 #add_depth_and_momentum2csv(swwfile, buildings_filename, buildings_filename_out)
 inundation_damage(swwfile, buildings_filename, buildings_filename_out)
+#inundation_damage(swwfile, buildings_filename, buildings_filename_out)
 
 print '\n Augmented building file written to %s \n' \
Index: /anuga_work/production/sydney_2006/report/interpretation.tex
===================================================================
--- /anuga_work/production/sydney_2006/report/interpretation.tex	(revision 4319)
+++ /anuga_work/production/sydney_2006/report/interpretation.tex	(revision 4320)
@@ -11,5 +11,5 @@
 \begin{tabular}{cc}
 \includegraphics[scale=0.4]{../report_figures/overallmodelv3.jpg} &
-\includegraphics[scale=0.4]{../report_figures/detailedmodel.jpg} &
+\includegraphics[scale=0.4]{../report_figures/detailedmodel.jpg} \\
 (a) & (b) \\
 \end{tabular}
