#
# 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.

"""

import numpy as num


def earthquake_tsunami(length, width, strike, depth, \
                       dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\
                       domain=None, verbose=False):
    """ return a function representing a tsunami event
    """

    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()
        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
    test = width*1000.0*sin(radians(dip)) - depth*1000
    
    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
        #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 = num.zeros(N, num.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):
            # 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

        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                                         
        """                                                                    
        
        U = num.zeros(9, num.float)
        DU = num.zeros(9, num.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) ))
            # GD CHANGE
            TT = atan( 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=====   VERTICAL 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=====   INCLINED 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) ))   
                # 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 )                                  
            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*R2
        QR = F3*Q/R5
        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 )

