#
# 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 (km)
 width   down-dip width of rupture area (km)
 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 (km)
 ns      the really used number of rectangular sources
 NSMAX   the upper limit of ns
Optional
 x0      x origin (0)
 y0      y origin (0)
 z0      z origin
 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.

"""
from __future__ import print_function
from __future__ import division

from builtins import str
from builtins import range
from builtins import object
from past.utils import old_div
from future.utils import raise_
import numpy as num
import anuga.utilities.log as log


def earthquake_tsunami(ns, NSMAX, length, width, strike, depth,
                       dip, xi, yi, z0, slip, rake,
                       domain=None, verbose=False):

    from anuga.abstract_2d_finite_volumes.quantity import Quantity
    from math import sin, radians

    #zrec0 = Quantity(domain)
    zrec0 = Quantity(domain)
    zrec0.set_values(z0)
    zrec = zrec0.get_vertex_values(xy=True)

    x0 = num.zeros(ns, float)
    y0 = num.zeros(ns, float)
    if ns == 1:
        x0[0] = xi
        y0[0] = yi
    else:
        x0 = xi
        y0 = yi
    if domain is not None:
        xllcorner = domain.geo_reference.get_xllcorner()
        yllcorner = domain.geo_reference.get_yllcorner()
        for i in range(0, ns):
            x0[i] = x0[i] - xllcorner  # fault origin (relative)
            y0[i] = y0[i] - yllcorner
    # a few temporary print statements
    if verbose is True:
        log.critical('\nThe Earthquake ...')
        log.critical('\tns: %s' % str(ns))
        log.critical('\tNSMAX: %s' % str(NSMAX))
        log.critical('\tLength: %s' % str(length))
        log.critical('\tDepth: %s' % str(depth))
        log.critical('\tStrike: %s' % str(strike))
        log.critical('\tWidth: %s' % str(width))
        log.critical('\tDip: %s' % str(dip))
        log.critical('\tSlip: %s' % str(slip))
        log.critical('\tx0: %s' % str(x0))
        log.critical('\ty0: %s' % str(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(ns=ns, NSMAX=NSMAX, length=length, width=width, dip=dip,
                      x0=x0, y0=y0, strike=strike, depth=depth,
                      slip=slip, rake=rake, zrec=zrec)

#
# 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(object):

    def __init__(self, ns, NSMAX, length, width, dip, x0, y0, strike,
                 depth, slip, rake, zrec):
        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
        self.ns = ns
        self.zrec = zrec

    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,z)
        produced by a submarine mass failure.
        """

        from math import sin, cos, radians, exp, cosh
        # ensure vectors x and y have the same length
        N = len(x)
        assert N == len(y)
        # nrec=N*N
        depth = self.depth
        dip = self.dip
        length = self.length
        width = self.width
        y0 = self.x0
        x0 = self.y0
        strike = self.strike
        dip = self.dip
        rake = self.rake
        dislocation = self.slip
        ns = self.ns
        zrec = self.zrec
        # initialization
        disp0 = num.zeros(3, float)
        strain0 = num.zeros(6, float)
        tilt0 = num.zeros(2, float)
        dislocations = num.zeros(ns, float)
        depths = num.zeros(ns, float)
        strikes = num.zeros(ns, float)
        lengths = num.zeros(ns, float)
        slips = num.zeros(ns, float)
        rakes = num.zeros(ns, float)
        widths = num.zeros(ns, float)
        dips = num.zeros(ns, float)
        strikes = num.zeros(ns, float)
        strikes = num.zeros(ns, float)
        strain = num.zeros((N, 6), float)
        disp = num.zeros((N, 3), float)
        tilt = num.zeros((N, 2), float)
        xs = num.zeros(ns, float)
        ys = num.zeros(ns, float)
        z = []
        if ns == 1:
            dislocations[0] = dislocation
            depths[0] = depth
            strikes[0] = strike
            lengths[0] = length
            rakes[0] = rake
            widths[0] = width
            dips[0] = dip
            try:
                xs[0] = x0
                ys[0] = y0
            except:
                xs[0] = x0[0]
                ys[0] = y0[0]
        else:
            dislocations = dislocation
            strikes = strike
            lengths = length
            rakes = rake
            widths = width
            dips = dip
            xs = x0
            ys = y0
            depths = depth
        # 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
        ALPHA = 0.5
        POT3 = 0.0
        POT4 = 0.0
        DISL3 = 0.0
        AL1 = 0.0
        AW2 = 0.0
        # rad=1.745329252e-2
        # zrec = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])   #get Zrec... has to be constant!!!
        #zrec=zreci.get_values(interpolation_points=[[x, y]])
        # zrec=0
        # Z=-zrec
        eps = 1.0e-6

#
        for irec in range(0, N):
            xrec = y
            yrec = x
            for i in range(0, len(zrec[0])):
                if (num.any(zrec[0][i] == yrec) and num.any(zrec[1][i] == xrec)):
                    Z = zrec[2][i]
                    Z = 0.001*Z
                    break
                else:
                    continue
            # zrec=zreci.get_values(interpolation_points=[[xrec[irec],yrec[irec]]],location='edges')

            for ist in range(0, ns):
                #st = radians(strikes[ist])
                st = radians(strikes[ist])
                csst = cos(st)
                ssst = sin(st)
                cs2st = cos(2.0*st)
                ss2st = sin(2.0*st)
#
                #di = radians(dips[ist])
                di = radians(dips[ist])
                csdi = cos(di)
                ssdi = sin(di)
#
                #ra= radians(rakes[ist])
                ra = radians(rakes[ist])
                csra = cos(ra)
                ssra = sin(ra)
                # transform from Aki's to Okada's system
                # first attribute x to north and y to east to match OKADA axis

                # X=0.001*(x[irec]-xs[ist])*ssst+0.001*(y[irec]-ys[ist])*csst
                # Y=-0.001*(x[irec]-xs[ist])*csst+0.001*(y[irec]-ys[ist])*ssst

                X = 0.001*((xrec[irec]-xs[ist])*csst+(yrec[irec]-ys[ist])*ssst)
                Y = 0.001*((xrec[irec]-xs[ist])*ssst-(yrec[irec]-ys[ist])*csst)
                DEPTH = depths[ist]
                DIP = dips[ist]
                if lengths[ist] == 0 and widths[ist] == 0:
                    #
                    #	      point source
                    #
                    POT1 = dislocations[ist]*csra
                    POT2 = dislocations[ist]*ssra
                    IRET = 1
                    self.DC3D0(ALPHA, X, Y, Z, DEPTH, DIP,
                               POT1, POT2, POT3, POT4)
                    UX = self.UX
                    UY = self.UY
                    UZ = self.UZ
                    UXX = self.UXX
                    UYX = self.UYX
                    UZX = self.UZX
                    UXY = self.UXY
                    UYY = self.UYY
                    UZY = self.UZY
                    UXZ = self.UXZ
                    UYZ = self.UYZ
                    UZZ = self.UZZ
                    IRET = self.IRET
                    if IRET == 1:
                        log.critical('There is a problem in Okada subroutine!')
                        break
                else:
                    #	        finite source
                    AL2 = lengths[ist]
                    AW1 = -widths[ist]
                    DISL1 = dislocations[ist]*csra
                    DISL2 = dislocations[ist]*ssra
                    if lengths[ist] == 0:
                        AL2 = widths[ist]*eps
                        DISL1 = old_div(DISL1, AL2)
                        DISL2 = old_div(DISL2, AL2)
                    elif widths[ist] == 0.0:
                        AW1 = -lengths[ist]*eps
                        DISL1 = old_div(DISL1, (-AW1))
                        DISL2 = old_div(DISL2, (-AW1))
                    IRET = 1
                    self.DC3D(ALPHA, X, Y, Z, DEPTH, DIP, AL1, AL2, AW1, AW2,
                              DISL1, DISL2, DISL3)
                    UX = self.UX
                    UY = self.UY
                    UZ = self.UZ
                    UXX = self.UXX
                    UYX = self.UYX
                    UZX = self.UZX
                    UXY = self.UXY
                    UYY = self.UYY
                    UZY = self.UZY
                    UXZ = self.UXZ
                    UYZ = self.UYZ
                    UZZ = self.UZZ
                    IRET = self.IRET
                    # if X==-6 and Y==-1:
                    # print UY
                    # print 'hello'
                    if IRET == 1:
                        print(' There is a problem in Okada subroutine!')
                        break
#
#	    transform from Okada's to Aki's system
#
                disp0[0] = UX*csst+UY*ssst
                disp0[1] = UX*ssst-UY*csst
                disp0[2] = -UZ

                tilt0[0] = -(UXZ*csst+UYZ*ssst)
                tilt0[1] = -(UXZ*ssst-UYZ*csst)
#
                strain0[0] = (UXX*csst*csst+UYY*ssst*ssst
                              + 0.5*(UXY+UYX)*ss2st)
                strain0[1] = (UXX*ssst*ssst+UYY*csst*csst
                              - 0.5*(UXY+UYX)*ss2st)
                strain0[2] = (UZZ)
                strain0[3] = (0.5*(UXX-UYY)*ss2st
                              - 0.5*(UXY+UYX)*cs2st)
                strain0[4] = (-0.5*(UZX+UXZ)*ssst
                              + 0.5*(UYZ+UZY)*csst)
                strain0[5] = (-0.5*(UZX+UXZ)*csst
                              - 0.5*(UYZ+UZY)*ssst)
                for j in range(0, 3):
                    disp[irec][j] = disp[irec][j] + disp0[j]
                for j in range(0, 2):
                    tilt[irec][j] = tilt[irec][j] + tilt0[j]
                for j in range(0, 6):
                    strain[irec][j] = strain[irec][j] + strain0[j]
            z.append(-disp[irec][2])
        return z

        # works on nautilus when have already done
        # f2py -c okada.f -m okada
        #z1 = okada(xr,yr,depth,length,width,dip,rake,slip)

    def DC3D0(self, ALPHA, X, Y, Z, DEPTH, DIP, POT1, POT2, POT3, POT4):
        """********************************************************************
        *****                                                          *****
        #*****    DISPLACEMENT AND STRAIN AT DEPTH                      *****
        #*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
        #*****                         CODED BY  Y.OKADA ... SEP.1991   *****
        #*****                         REVISED   Y.OKADA ... NOV.1991   *****
        #*****                                                          *****
        #********************************************************************
        #***** INPUT
        #*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)
        #*****   X,Y,Z : COORDINATE OF OBSERVING POINT
        #*****   DEPTH : SOURCE DEPTH
        #*****   DIP   : DIP-ANGLE (DEGREE)
        #*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
        #*****       POTENCY=(  MOMENT OF DOUBLE-COUPLE  )/MYU     FOR POT1,2
        #*****       POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMBDA  FOR POT3
        #*****       POTENCY=(INTENSITY OF LINEAR DIPOLE )/MYU     FOR POT4
        #***** OUTPUT
        #*****   UX, UY, UZ  : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) /
        #*****               :                     (UNIT OF X,Y,Z,DEPTH)**2  )
        #*****   UXX,UYX,UZX : X-DERIVATIVE ( UNIT= UNIT OF POTENCY) /
        #*****   UXY,UYY,UZY : Y-DERIVATIVE        (UNIT OF X,Y,Z,DEPTH)**3  )
        #*****   UXZ,UYZ,UZZ : Z-DERIVATIVE
        #*****   IRET        : RETURN CODE  ( =0....NORMAL,   =1....SINGULAR )"""

#        COMMON /C1/DUMMY(8),R,dumm(15)
#          DIMENSION  U(12),DUA(12),DUB(12),DUC(12)

        F0 = 0.0
        U = num.zeros((12, 1), float)
        DUA = num.zeros((12, 1), float)
        DUB = num.zeros((12, 1), float)
        DUC = num.zeros((12, 1), float)

        if Z > 0:
            print('(''0** POSITIVE Z WAS GIVEN IN SUB-DC3D0'')')
        for I in range(0, 12):
            U[I] = F0
            DUA[I] = F0
            DUB[I] = F0
            DUC[I] = F0

        AALPHA = ALPHA
        DDIP = DIP
        self.DCC0N0(AALPHA, DDIP)
# ======================================
# =====  REAL-SOURCE CONTRIBUTION  =====
# ======================================
        XX = X
        YY = Y
        ZZ = Z
        DD = DEPTH+Z
        self.DCCON1(XX, YY, DD)
        R = self.R
        if R == F0:
            UX = F0
            UY = F0
            UZ = F0
            UXX = F0
            UYX = F0
            UZX = F0
            UXY = F0
            UYY = F0
            UZY = F0
            UXZ = F0
            UYZ = F0
            UZZ = F0
            IRET = 1
            return
        else:
            PP1 = POT1
            PP2 = POT2
            PP3 = POT3
            PP4 = POT4

            self.UA0(XX, YY, DD, PP1, PP2, PP3, PP4)
            DUA = self.DUA

            for I in range(0, 12):
                if I < 10:
                    U[I] = U[I]-DUA[I]
                if I >= 10:
                    U[I] = U[I]+DUA[I]
# =======================================
# =====  IMAGE-SOURCE CONTRIBUTION  =====
# =======================================
            DD = DEPTH-Z
            self.DCCON1(XX, YY, DD)
            self.UA0(XX, YY, DD, PP1, PP2, PP3, PP4)
            DUA = self.DUA
            self.UB0(XX, YY, DD, ZZ, PP1, PP2, PP3, PP4)
            DUB = self.DUB
            self.UC0(XX, YY, DD, ZZ, PP1, PP2, PP3, PP4)
            DUC = self.DUC

# -----
            for I in range(0, 12):
                DU = DUA[I]+DUB[I]+ZZ*DUC[I]
                if I >= 9:
                    DU = DU+DUC[I-9]
                U[I] = U[I]+DU
# =====
            self.UX = U[0]
            self.UY = U[1]
            self.UZ = U[2]
            self.UXX = U[3]
            self.UYX = U[4]
            self.UZX = U[5]
            self.UXY = U[6]
            self.UYY = U[7]
            self.UZY = U[8]
            self.UXZ = U[9]
            self.UYZ = U[10]
            self.UZZ = U[11]
            self.IRET = 0

    def UA0(self, X, Y, D, POT1, POT2, POT3, POT4):

        # SUBROUTINE  UA0(X,Y,D,POT1,POT2,POT3,POT4,U)
     #     IMPLICIT REAL*8 (A-H,O-Z)
     #     DIMENSION U(12),DU(12)
        """********************************************************************
        C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-A)             *****
        C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
        C********************************************************************
        C***** INPUT
        C*****   X,Y,D : STATION COORDINATES IN FAULT SYSTEM
        C*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
        C***** OUTPUT
        C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES    """
#
# COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
#      COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,
#    *           UY,VY,WY,UZ,VZ,WZ

        F0 = 0.0
        F1 = 1.0
        F3 = 3.0
        PI2 = 6.283185307179586

        ALP1 = self.ALP1
        ALP2 = self.ALP2
        ALP3 = self.ALP3
        ALP4 = self.ALP4
        ALP5 = self.ALP4
        SD = self.SD
        CD = self.CD
        SDSD = self.SDSD
        CDCD = self.CDCD
        SDCD = self.SDCD
        S2D = self.S2D
        C2D = self.C2D
        P = self.P
        Q = self.Q
        S = self.S
        T = self.T
        XY = self.XY
        X2 = self.X2
        Y2 = self.Y2
        D2 = self.D2
        R = self.R
        R2 = self.R2
        R3 = self.R3
        R5 = self.R5
        QR = self.QR
        QRX = self.QRX
        A3 = self.A3
        A5 = self.A5
        B3 = self.B3
        C3 = self.C3
        UY = self.UY
        VY = self.VY
        WY = self.WY
        UZ = self.UZ
        VZ = self.VZ
        WZ = self.WZ

        DUA = num.zeros((12, 1), float)
        DU = num.zeros((12, 1), float)
        U = num.zeros((12, 1), float)
# -----
        for I in range(0, 12):
            U[I] = F0
# ======================================
# =====  STRIKE-SLIP CONTRIBUTION  =====
# ======================================
        if POT1 != F0:
            DU[0] = old_div(ALP1*Q, R3) + ALP2*X2*QR
            DU[1] = old_div(ALP1*X, R3)*SD + ALP2*XY*QR
            DU[2] = old_div(-ALP1*X, R3)*CD + ALP2*X*D*QR
            DU[3] = X*QR*(-ALP1 + ALP2*(F1+A5))
            DU[4] = old_div(ALP1*A3, R3)*SD + ALP2*Y*QR*A5
            DU[5] = old_div(-ALP1*A3, R3)*CD + ALP2*D*QR*A5
            DU[6] = ALP1*(old_div(SD, R3)-Y*QR) + old_div(ALP2*F3*X2, R5)*UY
            DU[7] = old_div(F3*X, R5)*(-ALP1*Y*SD + ALP2*(Y*UY+Q))
            DU[8] = old_div(F3*X, R5)*(ALP1*Y*CD + ALP2*D*UY)
            DU[9] = ALP1*(old_div(CD, R3)+D*QR) + old_div(ALP2*F3*X2, R5)*UZ
            DU[10] = old_div(F3*X, R5)*(ALP1*D*SD + ALP2*Y*UZ)
            DU[11] = old_div(F3*X, R5)*(-ALP1*D*CD + ALP2*(D*UZ-Q))
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT1, PI2)*DU[I]

# ===================================
# =====  DIP-SLIP CONTRIBUTION  =====
# ===================================
        if POT2 != F0:
            DU[0] = ALP2*X*P*QR
            DU[1] = old_div(ALP1*S, R3) + ALP2*Y*P*QR
            DU[2] = old_div(-ALP1*T, R3) + ALP2*D*P*QR
            DU[3] = ALP2*P*QR*A5
            DU[4] = old_div(-ALP1*F3*X*S, R5) - ALP2*Y*P*QRX
            DU[5] = old_div(ALP1*F3*X*T, R5) - ALP2*D*P*QRX
            DU[6] = old_div(ALP2*F3*X, R5)*VY
            DU[7] = ALP1*(old_div(S2D, R3)-old_div(F3*Y*S, R5)) + \
                ALP2*(old_div(F3*Y, R5)*VY+P*QR)
            DU[8] = -ALP1*(old_div(C2D, R3)-old_div(F3*Y*T, R5)
                           ) + old_div(ALP2*F3*D, R5)*VY
            DU[9] = old_div(ALP2*F3*X, R5)*VZ
            DU[10] = ALP1*(old_div(C2D, R3)+old_div(F3*D*S, R5)
                           ) + old_div(ALP2*F3*Y, R5)*VZ
            DU[11] = ALP1*(old_div(S2D, R3)-old_div(F3*D*T, R5)
                           ) + ALP2*(old_div(F3*D, R5)*VZ-P*QR)
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT2, PI2)*DU[I]

# ========================================
# =====  TENSILE-FAULT CONTRIBUTION  =====
# ========================================
        if POT3 != F0:
            DU[0] = old_div(ALP1*X, R3) - ALP2*X*Q*QR
            DU[1] = old_div(ALP1*T, R3) - ALP2*Y*Q*QR
            DU[2] = old_div(ALP1*S, R3) - ALP2*D*Q*QR
            DU[3] = old_div(ALP1*A3, R3) - ALP2*Q*QR*A5
            DU[4] = old_div(-ALP1*F3*X*T, R5) + ALP2*Y*Q*QRX
            DU[5] = old_div(-ALP1*F3*X*S, R5) + ALP2*D*Q*QRX
            DU[6] = old_div(-ALP1*F3*XY, R5) - ALP2*X*QR*WY
            DU[7] = ALP1*(old_div(C2D, R3)-old_div(F3*Y*T, R5)) - \
                ALP2*(Y*WY+Q)*QR
            DU[8] = ALP1*(old_div(S2D, R3)-old_div(F3*Y*S, R5)) - ALP2*D*QR*WY
            DU[9] = old_div(ALP1*F3*X*D, R5) - ALP2*X*QR*WZ
            DU[10] = -ALP1*(old_div(S2D, R3) -
                            old_div(F3*D*T, R5)) - ALP2*Y*QR*WZ
            DU[11] = ALP1*(old_div(C2D, R3)+old_div(F3*D*S, R5)
                           ) - ALP2*(D*WZ-Q)*QR
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT3, PI2)*DU[I]

# =========================================
# =====  INFLATE SOURCE CONTRIBUTION  =====
# =========================================
        if POT4 != F0:
            DU[0] = old_div(-ALP1*X, R3)
            DU[1] = old_div(-ALP1*Y, R3)
            DU[2] = old_div(-ALP1*D, R3)
            DU[3] = old_div(-ALP1*A3, R3)
            DU[4] = old_div(ALP1*F3*XY, R5)
            DU[5] = old_div(ALP1*F3*X*D, R5)
            DU[6] = DU[4]
            DU[7] = old_div(-ALP1*B3, R3)
            DU[8] = old_div(ALP1*F3*Y*D, R5)
            DU[9] = -DU[5]
            DU[10] = -DU[8]
            DU[11] = old_div(ALP1*C3, R3)
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT4, PI2)*DU[I]

        # for I in range(0,12):
            # DUA[I]=U[I]

        self.DUA = U

    def UB0(self, X, Y, D, Z, POT1, POT2, POT3, POT4):
     #     SUBROUTINE  UB0(X,Y,D,Z,POT1,POT2,POT3,POT4,U)
      #    IMPLICIT REAL*8 (A-H,O-Z)
       #   DIMENSION U(12),DU(12)
        #
        """********************************************************************
        C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****
        C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
        C********************************************************************
        C
        C***** INPUT
        C*****   X,Y,D,Z : STATION COORDINATES IN FAULT SYSTEM
        C*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
        C***** OUTPUT
        C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES    """
#
 #     COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
 #     COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,
 #    *           UY,VY,WY,UZ,VZ,WZ
# F0,F1,F2,F3,F4,F5,F8,F9
#     *        /0.D0,1.D0,2.D0,3.D0,4.D0,5.D0,8.D0,9.D0/
#      DATA PI2/6.283185307179586D0/

        DUB = num.zeros((12, 1), float)
        DU = num.zeros((12, 1), float)
        U = num.zeros((12, 1), float)

        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

        ALP1 = self.ALP1
        ALP2 = self.ALP2
        ALP3 = self.ALP3
        ALP4 = self.ALP4
        ALP5 = self.ALP4
        SD = self.SD
        CD = self.CD
        SDSD = self.SDSD
        CDCD = self.CDCD
        SDCD = self.SDCD
        S2D = self.S2D
        C2D = self.C2D
        P = self.P
        Q = self.Q
        S = self.S
        T = self.T
        XY = self.XY
        X2 = self.X2
        Y2 = self.Y2
        D2 = self.D2
        R = self.R
        R2 = self.R2
        R3 = self.R3
        R5 = self.R5
        QR = self.QR
        QRX = self.QRX
        A3 = self.A3
        A5 = self.A5
        B3 = self.B3
        C3 = self.C3
        UY = self.UY
        VY = self.VY
        WY = self.WY
        UZ = self.UZ
        VZ = self.VZ
        WZ = self.WZ

# -----
        C = D+Z
        RD = R+D
        D12 = old_div(F1, (R*RD*RD))
        D32 = old_div(D12*(F2*R+D), R2)
        D33 = old_div(D12*(F3*R+D), (R2*RD))
        D53 = old_div(D12*(F8*R2+F9*R*D+F3*D2), (R2*R2*RD))
        D54 = old_div(D12*(F5*R2+F4*R*D+D2), R3)*D12
# -----
        FI1 = Y*(D12-X2*D33)
        FI2 = X*(D12-Y2*D33)
        FI3 = old_div(X, R3)-FI2
        FI4 = -XY*D32
        FI5 = old_div(F1, (R*RD))-X2*D32
        FJ1 = -F3*XY*(D33-X2*D54)
        FJ2 = old_div(F1, R3)-F3*D12+F3*X2*Y2*D54
        FJ3 = old_div(A3, R3)-FJ2
        FJ4 = old_div(-F3*XY, R5)-FJ1
        FK1 = -Y*(D32-X2*D53)
        FK2 = -X*(D32-Y2*D53)
        FK3 = old_div(-F3*X*D, R5)-FK2
# -----
        for I in range(0, 12):
            U[I] = F0
# ======================================
# =====  STRIKE-SLIP CONTRIBUTION  =====
# ======================================
        if POT1 != F0:
            DU[0] = -X2*QR - ALP3*FI1*SD
            DU[1] = -XY*QR - ALP3*FI2*SD
            DU[2] = -C*X*QR - ALP3*FI4*SD
            DU[3] = -X*QR*(F1+A5) - ALP3*FJ1*SD
            DU[4] = -Y*QR*A5 - ALP3*FJ2*SD
            DU[5] = -C*QR*A5 - ALP3*FK1*SD
            DU[6] = old_div(-F3*X2, R5)*UY - ALP3*FJ2*SD
            DU[7] = old_div(-F3*XY, R5)*UY-X*QR - ALP3*FJ4*SD
            DU[8] = old_div(-F3*C*X, R5)*UY - ALP3*FK2*SD
            DU[9] = old_div(-F3*X2, R5)*UZ + ALP3*FK1*SD
            DU[10] = old_div(-F3*XY, R5)*UZ + ALP3*FK2*SD
            DU[11] = old_div(F3*X, R5)*(-C*UZ + ALP3*Y*SD)
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT1, PI2)*DU[I]
# ===================================
# =====  DIP-SLIP CONTRIBUTION  =====
# ===================================
        if POT2 != F0:
            DU[0] = -X*P*QR + ALP3*FI3*SDCD
            DU[1] = -Y*P*QR + ALP3*FI1*SDCD
            DU[2] = -C*P*QR + ALP3*FI5*SDCD
            DU[3] = -P*QR*A5 + ALP3*FJ3*SDCD
            DU[4] = Y*P*QRX + ALP3*FJ1*SDCD
            DU[5] = C*P*QRX + ALP3*FK3*SDCD
            DU[6] = old_div(-F3*X, R5)*VY + ALP3*FJ1*SDCD
            DU[7] = old_div(-F3*Y, R5)*VY-P*QR + ALP3*FJ2*SDCD
            DU[8] = old_div(-F3*C, R5)*VY + ALP3*FK1*SDCD
            DU[9] = old_div(-F3*X, R5)*VZ - ALP3*FK3*SDCD
            DU[10] = old_div(-F3*Y, R5)*VZ - ALP3*FK1*SDCD
            DU[11] = old_div(-F3*C, R5)*VZ + old_div(ALP3*A3, R3)*SDCD
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT2, PI2)*DU[I]
# ========================================
# =====  TENSILE-FAULT CONTRIBUTION  =====
# ========================================
        if POT3 != F0:
            DU[0] = X*Q*QR - ALP3*FI3*SDSD
            DU[1] = Y*Q*QR - ALP3*FI1*SDSD
            DU[2] = C*Q*QR - ALP3*FI5*SDSD
            DU[3] = Q*QR*A5 - ALP3*FJ3*SDSD
            DU[4] = -Y*Q*QRX - ALP3*FJ1*SDSD
            DU[5] = -C*Q*QRX - ALP3*FK3*SDSD
            DU[6] = X*QR*WY - ALP3*FJ1*SDSD
            DU[7] = QR*(Y*WY+Q) - ALP3*FJ2*SDSD
            DU[8] = C*QR*WY - ALP3*FK1*SDSD
            DU[9] = X*QR*WZ + ALP3*FK3*SDSD
            DU[10] = Y*QR*WZ + ALP3*FK1*SDSD
            DU[11] = C*QR*WZ - old_div(ALP3*A3, R3)*SDSD
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT3, PI2)*DU[I]

# =========================================
# =====  INFLATE SOURCE CONTRIBUTION  =====
# =========================================
        if POT4 != F0:
            DU[0] = old_div(ALP3*X, R3)
            DU[1] = old_div(ALP3*Y, R3)
            DU[2] = old_div(ALP3*D, R3)
            DU[3] = old_div(ALP3*A3, R3)
            DU[4] = old_div(-ALP3*F3*XY, R5)
            DU[5] = old_div(-ALP3*F3*X*D, R5)
            DU[6] = DU[4]
            DU[7] = old_div(ALP3*B3, R3)
            DU[8] = old_div(-ALP3*F3*Y*D, R5)
            DU[9] = -DU[5]
            DU[10] = -DU[8]
            DU[11] = old_div(-ALP3*C3, R3)
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT4, PI2)*DU[I]

        # for I in range(0,12):
            # DUB[I]=U[I]
        self.DUB = U

    def UC0(self, X, Y, D, Z, POT1, POT2, POT3, POT4):
      #    SUBROUTINE  UC0(X,Y,D,Z,POT1,POT2,POT3,POT4,U)
      #    IMPLICIT REAL*8 (A-H,O-Z)
      #    DIMENSION U(12),DU(12)
        """********************************************************************
        C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****
        C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
        C********************************************************************
        C***** INPUT
        C*****   X,Y,D,Z : STATION COORDINATES IN FAULT SYSTEM
        C*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
        C***** OUTPUT
        C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES"""

#      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
 #     COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,um(6)
#      DATA F0,F1,F2,F3,F5,F7,F10,F15
#     *        /0.D0,1.D0,2.D0,3.D0,5.D0,7.D0,10.D0,15.D0/
#      DATA PI2/6.283185307179586D0/

        DUC = num.zeros((12, 1), float)
        DU = num.zeros((12, 1), float)
        U = num.zeros((12, 1), float)

        F0 = 0.0
        F1 = 1.0
        F2 = 2.0
        F3 = 3.0
        F5 = 5.0
        F7 = 7.0
        F10 = 10.0
        F15 = 15.0
        PI2 = 6.283185307179586

        ALP1 = self.ALP1
        ALP2 = self.ALP2
        ALP3 = self.ALP3
        ALP4 = self.ALP4
        ALP5 = self.ALP4
        SD = self.SD
        CD = self.CD
        SDSD = self.SDSD
        CDCD = self.CDCD
        SDCD = self.SDCD
        S2D = self.S2D
        C2D = self.C2D
        P = self.P
        Q = self.Q
        S = self.S
        T = self.T
        XY = self.XY
        X2 = self.X2
        Y2 = self.Y2
        D2 = self.D2
        R = self.R
        R2 = self.R2
        R3 = self.R3
        R5 = self.R5
        QR = self.QR
        QRX = self.QRX
        A3 = self.A3
        A5 = self.A5
        B3 = self.B3
        C3 = self.C3
        UY = self.UY
        VY = self.VY
        WY = self.WY
        UZ = self.UZ
        VZ = self.VZ
        WZ = self.WZ

# -----
        C = D+Z
        Q2 = Q*Q
        R7 = R5*R2
        A7 = F1-old_div(F7*X2, R2)
        B5 = F1-old_div(F5*Y2, R2)
        B7 = F1-old_div(F7*Y2, R2)
        C5 = F1-old_div(F5*D2, R2)
        C7 = F1-old_div(F7*D2, R2)
        D7 = F2-old_div(F7*Q2, R2)
        QR5 = old_div(F5*Q, R2)
        QR7 = old_div(F7*Q, R2)
        DR5 = old_div(F5*D, R2)
# -----
        for I in range(0, 12):
            U[I] = F0
# ======================================
# =====  STRIKE-SLIP CONTRIBUTION  =====
# ======================================
        if POT1 != F0:
            DU[0] = old_div(-ALP4*A3, R3)*CD + ALP5*C*QR*A5
            DU[1] = old_div(F3*X, R5)*(ALP4*Y*CD + ALP5*C*(SD-Y*QR5))
            DU[2] = old_div(F3*X, R5)*(-ALP4*Y*SD + ALP5*C*(CD+D*QR5))
            DU[3] = old_div(ALP4*F3*X, R5)*(F2+A5)*CD - ALP5*C*QRX*(F2+A7)
            DU[4] = old_div(F3, R5)*(ALP4*Y*A5*CD + ALP5*C*(A5*SD-Y*QR5*A7))
            DU[5] = old_div(F3, R5)*(-ALP4*Y*A5*SD + ALP5*C*(A5*CD+D*QR5*A7))
            DU[6] = DU[4]
            DU[7] = old_div(F3*X, R5)*(ALP4*B5*CD -
                                       old_div(ALP5*F5*C, R2)*(F2*Y*SD+Q*B7))
            DU[8] = old_div(F3*X, R5)*(-ALP4*B5*SD +
                                       old_div(ALP5*F5*C, R2)*(D*B7*SD-Y*C7*CD))
            DU[9] = old_div(F3, R5) * (-ALP4*D*A5*CD + ALP5*C*(A5*CD+D*QR5*A7))
            DU[10] = old_div(F15*X, R7)*(ALP4*Y*D*CD +
                                         ALP5*C*(D*B7*SD-Y*C7*CD))
            DU[11] = old_div(F15*X, R7)*(-ALP4*Y*D*SD + ALP5*C*(F2*D*CD-Q*C7))
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT1, PI2)*DU[I]
# ===================================
# =====  DIP-SLIP CONTRIBUTION  =====
# ===================================
        if POT2 != F0:
            DU[0] = old_div(ALP4*F3*X*T, R5) - ALP5*C*P*QRX
            DU[1] = old_div(-ALP4, R3)*(C2D-old_div(F3*Y*T, R2)
                                        ) + old_div(ALP5*F3*C, R5)*(S-Y*P*QR5)
            DU[2] = old_div(-ALP4*A3, R3)*SDCD + \
                old_div(ALP5*F3*C, R5)*(T+D*P*QR5)
            DU[3] = old_div(ALP4*F3*T, R5)*A5 - old_div(ALP5*F5*C*P*QR, R2)*A7
            DU[4] = old_div(F3*X, R5)*(ALP4*(C2D-old_div(F5*Y*T, R2)
                                             )-old_div(ALP5*F5*C, R2)*(S-Y*P*QR7))
            DU[5] = old_div(F3*X, R5)*(ALP4*(F2+A5)*SDCD -
                                       old_div(ALP5*F5*C, R2)*(T+D*P*QR7))
            DU[6] = DU[4]
            DU[7] = (old_div(F3, R5)*(ALP4*(F2*Y*C2D+T*B5)
                                      + ALP5*C*(S2D-old_div(F10*Y*S, R2)-P*QR5*B7)))
            DU[8] = old_div(F3, R5)*(ALP4*Y*A5*SDCD-ALP5 *
                                     C*((F3+A5)*C2D+Y*P*DR5*QR7))
            DU[9] = old_div(F3*X, R5)*(-ALP4*(S2D-T*DR5) -
                                       old_div(ALP5*F5*C, R2)*(T+D*P*QR7))
            DU[10] = (old_div(F3, R5)*(-ALP4*(D*B5*C2D+Y*C5*S2D)
                                       - ALP5*C*((F3+A5)*C2D+Y*P*DR5*QR7)))
            DU[11] = old_div(F3, R5)*(-ALP4*D*A5*SDCD-ALP5 *
                                      C*(S2D-old_div(F10*D*T, R2)+P*QR5*C7))
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT2, PI2)*DU[I]
# ========================================
# =====  TENSILE-FAULT CONTRIBUTION  =====
# ========================================
        if POT3 != F0:
            DU[0] = old_div(F3*X, R5)*(-ALP4*S + ALP5*(C*Q*QR5-Z))
            DU[1] = old_div(ALP4, R3)*(S2D-old_div(F3*Y*S, R2)) + \
                old_div(ALP5*F3, R5)*(C*(T-Y+Y*Q*QR5)-Y*Z)
            DU[2] = old_div(-ALP4, R3)*(F1-A3*SDSD) - \
                old_div(ALP5*F3, R5)*(C*(S-D+D*Q*QR5)-D*Z)
            DU[3] = old_div(-ALP4*F3*S, R5)*A5 + ALP5 * \
                (C*QR*QR5*A7-old_div(F3*Z, R5)*A5)
            DU[4] = (old_div(F3*X, R5)*(-ALP4*(S2D-old_div(F5*Y*S, R2))
                                        - old_div(ALP5*F5, R2)*(C*(T-Y+Y*Q*QR7)-Y*Z)))
            DU[5] = (old_div(F3*X, R5)*(ALP4*(F1-(F2+A5)*SDSD)
                                        + old_div(ALP5*F5, R2)*(C*(S-D+D*Q*QR7)-D*Z)))
            DU[6] = DU[4]
            DU[7] = (old_div(F3, R5)*(-ALP4*(F2*Y*S2D+S*B5)
                                      - ALP5*(C*(F2*SDSD+old_div(F10*Y*(T-Y), R2)-Q*QR5*B7)+Z*B5)))
            DU[8] = (old_div(F3, R5)*(ALP4*Y*(F1-A5*SDSD)
                                      + ALP5*(C*(F3+A5)*S2D-Y*DR5*(C*D7+Z))))
            DU[9] = (old_div(F3*X, R5)*(-ALP4*(C2D+S*DR5)
                                        + ALP5*(old_div(F5*C, R2)*(S-D+D*Q*QR7)-F1-Z*DR5)))
            DU[10] = (old_div(F3, R5)*(ALP4*(D*B5*S2D-Y*C5*C2D)
                                       + ALP5*(C*((F3+A5)*S2D-Y*DR5*D7)-Y*(F1+Z*DR5))))
            DU[11] = (old_div(F3, R5)*(-ALP4*D*(F1-A5*SDSD)
                                       - ALP5*(C*(C2D+old_div(F10*D*(S-D), R2)-Q*QR5*C7)+Z*(F1+C5))))
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT3, PI2)*DU[I]
# =========================================
# =====  INFLATE SOURCE CONTRIBUTION  =====
# =========================================
        if POT4 != F0:
            DU[0] = old_div(ALP4*F3*X*D, R5)
            DU[1] = old_div(ALP4*F3*Y*D, R5)
            DU[2] = old_div(ALP4*C3, R3)
            DU[3] = old_div(ALP4*F3*D, R5)*A5
            DU[4] = old_div(-ALP4*F15*XY*D, R7)
            DU[5] = old_div(-ALP4*F3*X, R5)*C5
            DU[6] = DU[4]
            DU[7] = old_div(ALP4*F3*D, R5)*B5
            DU[8] = old_div(-ALP4*F3*Y, R5)*C5
            DU[9] = DU[5]
            DU[10] = DU[8]
            DU[11] = old_div(ALP4*F3*D, R5)*(F2+C5)
            for I in range(0, 12):
                U[I] = U[I]+old_div(POT4, PI2)*DU[I]

        # for I in range(0,12):
            # DUC[I]=U[I]
        self.DUC = U

    def DC3D(self, ALPHA, X, Y, Z, DEPTH, DIP, AL1, AL2, AW1, AW2, DISL1, DISL2, DISL3):

        #      SUBROUTINE  DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
        #     *              AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
     #    *              UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET)
        #      IMPLICIT REAL*8 (A-H,O-Z)
        #      REAL*4   ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
        #     *         UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,EPS
        """********************************************************************
        C*****                                                          *****
        C*****    DISPLACEMENT AND STRAIN AT DEPTH                      *****
        C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
        C*****              CODED BY  Y.OKADA ... SEP.1991              *****
        C*****              REVISED ... NOV.1991, APR.1992, MAY.1993,   *****
        C*****                          JUL.1993                        *****
        C********************************************************************
        C
        C***** INPUT
        C*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)
        C*****   X,Y,Z : COORDINATE OF OBSERVING POINT
        C*****   DEPTH : DEPTH OF REFERENCE POINT
        C*****   DIP   : DIP-ANGLE (DEGREE)
        C*****   AL1,AL2   : FAULT LENGTH RANGE
        C*****   AW1,AW2   : FAULT WIDTH RANGE
        C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
        C
        C***** OUTPUT
        C*****   UX, UY, UZ  : DISPLACEMENT ( UNIT=(UNIT OF DISL)
        C*****   UXX,UYX,UZX : X-DERIVATIVE ( UNIT=(UNIT OF DISL) /
        C*****   UXY,UYY,UZY : Y-DERIVATIVE        (UNIT OF X,Y,Z,DEPTH,AL,AW) )
        C*****   UXZ,UYZ,UZZ : Z-DERIVATIVE
        C*****   IRET        : RETURN CODE  ( =0....NORMAL,   =1....SINGULAR )  """
#
#      COMMON /C0/DUMMY(5),SD,CD,dumm(5)
 #     DIMENSION  XI(2),ET(2),KXI(2),KET(2)
 #     DIMENSION  U(12),DU(12),DUA(12),DUB(12),DUC(12)
        from math import sqrt
        F0 = 0.0
        EPS = 0.000001

        XI = num.zeros(2, float)
        ET = num.zeros(2, float)
        KXI = num.zeros(2, float)
        KET = num.zeros(2, float)
        U = num.zeros(12, float)
        DU = num.zeros(12, float)
        DUA = num.zeros(12, float)
        DUB = num.zeros(12, float)
        DUC = num.zeros(12, float)

# -----
        if Z > 0:
            log.critical('** POSITIVE Z WAS GIVEN IN SUB-DC3D')

        for I in range(0, 12):
            U[I] = F0
            DUA[I] = F0
            DUB[I] = F0
            DUC[I] = F0

        AALPHA = ALPHA
        DDIP = DIP
        self.DCC0N0(AALPHA, DDIP)
        CD = self.CD
        SD = self.SD

# -----
        ZZ = Z
        DD1 = DISL1
        DD2 = DISL2
        DD3 = DISL3
        XI[0] = X-AL1
        XI[1] = X-AL2
        # if X==-6 and Y==-1:
        # print AL1
        # print AL2
        # print 'hello'

        if abs(XI[0]) < EPS:
            XI[0] = F0
        if abs(XI[1]) < EPS:
            XI[1] = F0
# ======================================
# =====  REAL-SOURCE CONTRIBUTION  =====
# ======================================
        D = DEPTH+Z
        P = Y*CD+D*SD
        Q = Y*SD-D*CD
        ET[0] = P-AW1
        ET[1] = P-AW2
        if abs(Q) < EPS:
            Q = F0
        if abs(ET[0]) < EPS:
            ET[0] = F0
        if abs(ET[1]) < EPS:
            ET[1] = F0

# --------------------------------
# ----- REJECT SINGULAR CASE -----
# --------------------------------
# ----- ON FAULT EDGE
        if (Q == F0 and ((XI[0]*XI[1] < F0 and ET[0]*ET[1] == F0)
                         or (ET[0]*ET[1] < F0 and XI[0]*XI[1] == F0))):
            self.UX = F0
            self.UY = F0
            self.UZ = F0
            self.UXX = F0
            self.UYX = F0
            self.UZX = F0
            self.UXY = F0
            self.UYY = F0
            self.UZY = F0
            self.UXZ = F0
            self.UYZ = F0
            self.UZZ = F0
            self.IRET = 1
            return
# ----- ON NEGATIVE EXTENSION OF FAULT EDGE
        KXI[0] = 0
        KXI[1] = 0
        KET[0] = 0
        KET[1] = 0
        R12 = sqrt(XI[0]*XI[0]+ET[1]*ET[1]+Q*Q)
        R21 = sqrt(XI[1]*XI[1]+ET[0]*ET[0]+Q*Q)
        R22 = sqrt(XI[1]*XI[1]+ET[1]*ET[1]+Q*Q)
        if XI[0] < F0 and R21+XI[1] < EPS:
            KXI[0] = 1
        if XI[0] < F0 and R22+XI[1] < EPS:
            KXI[1] = 1
        if ET[0] < F0 and R12+ET[1] < EPS:
            KET[0] = 1
        if ET[0] < F0 and R22+ET[1] < EPS:
            KET[1] = 1
# =====
        for K in range(0, 2):
            for J in range(0, 2):
                self.DCCON2(XI[J], ET[K], Q, SD, CD, KXI[K], KET[J])
                self.UA(XI[J], ET[K], Q, DD1, DD2, DD3)
                DUA = self.DUA
                for I in range(0, 10, 3):
                    DU[I] = -DUA[I]
                    DU[I+1] = -DUA[I+1]*CD+DUA[I+2]*SD
                    DU[I+2] = -DUA[I+1]*SD-DUA[I+2]*CD
                    if I == 9:
                        DU[I] = -DU[I]
                        DU[I+1] = -DU[I+1]
                        DU[I+2] = -DU[I+2]
                    else:
                        continue
                for I in range(0, 12):
                    if (J+K) != 1:
                        U[I] = U[I]+DU[I]
                    if (J+K) == 1:
                        U[I] = U[I]-DU[I]

# =======================================
# =====  IMAGE-SOURCE CONTRIBUTION  =====
# =======================================
        D = DEPTH-Z
        P = Y*CD+D*SD
        Q = Y*SD-D*CD
        ET[0] = P-AW1
        ET[1] = P-AW2
        if abs(Q) < EPS:
            Q = F0
        if abs(ET[0]) < EPS:
            ET[0] = F0
        if abs(ET[1]) < EPS:
            ET[1] = F0
# --------------------------------
# ----- REJECT SINGULAR CASE -----
# --------------------------------
# ----- ON FAULT EDGE
        if (Q == F0 and ((XI[0]*XI[1] < F0 and ET[0]*ET[1] == F0)
                         or (ET[0]*ET[1] < F0 and XI[0]*XI[1] == F0))):
            self.UX = F0
            self.UY = F0
            self.UZ = F0
            self.UXX = F0
            self.UYX = F0
            self.UZX = F0
            self.UXY = F0
            self.UYY = F0
            self.UZY = F0
            self.UXZ = F0
            self.UYZ = F0
            self.UZZ = F0
            self.IRET = 1
            return
# ----- ON NEGATIVE EXTENSION OF FAULT EDGE
        KXI[0] = 0
        KXI[1] = 0
        KET[0] = 0
        KET[1] = 0
        R12 = sqrt(XI[0]*XI[0]+ET[1]*ET[1]+Q*Q)
        R21 = sqrt(XI[1]*XI[1]+ET[0]*ET[0]+Q*Q)
        R22 = sqrt(XI[1]*XI[1]+ET[1]*ET[1]+Q*Q)
        if XI[0] < F0 and R21+XI[1] < EPS:
            KXI[0] = 1
        if XI[0] < F0 and R22+XI[1] < EPS:
            KXI[1] = 1
        if ET[0] < F0 and R12+ET[1] < EPS:
            KET[0] = 1
        if ET[0] < F0 and R22+ET[1] < EPS:
            KET[1] = 1
# =====
        for K in range(0, 2):
            for J in range(0, 2):
                self.DCCON2(XI[J], ET[K], Q, SD, CD, KXI[K], KET[J])
                self.UA(XI[J], ET[K], Q, DD1, DD2, DD3)
                DUA = self.DUA
                self.UB(XI[J], ET[K], Q, DD1, DD2, DD3)
                DUB = self.DUB
                self.UC(XI[J], ET[K], Q, ZZ, DD1, DD2, DD3)
                DUC = self.DUC
# -----
                for I in range(0, 10, 3):
                    DU[I] = DUA[I]+DUB[I]+Z*DUC[I]
                    DU[I+1] = ((DUA[I+1]+DUB[I+1]+Z*DUC[I+1])*CD
                               - (DUA[I+2]+DUB[I+2]+Z*DUC[I+2])*SD)
                    DU[I+2] = ((DUA[I+1]+DUB[I+1]-Z*DUC[I+1])*SD
                               + (DUA[I+2]+DUB[I+2]-Z*DUC[I+2])*CD)
                    if I == 9:
                        DU[9] = DU[9]+DUC[0]
                        DU[10] = DU[10]+DUC[1]*CD-DUC[2]*SD
                        DU[11] = DU[11]-DUC[1]*SD-DUC[2]*CD
                for I in range(0, 12):
                    if (J+K) != 1:
                        U[I] = U[I]+DU[I]
                    if (J+K) == 1:
                        U[I] = U[I]-DU[I]
# =====
        self.UX = U[0]
        self.UY = U[1]
        self.UZ = U[2]
        self.UXX = U[3]
        self.UYX = U[4]
        self.UZX = U[5]
        self.UXY = U[6]
        self.UYY = U[7]
        self.UZY = U[8]
        self.UXZ = U[9]
        self.UYZ = U[10]
        self.UZZ = U[11]
        self.IRET = 0

    def UA(self, XI, ET, Q, DISL1, DISL2, DISL3):
        #      IMPLICIT REAL*8 (A-H,O-Z)
        #      DIMENSION U(12),DU(12)
        #
        # ********************************************************************
        # *****    DISPLACEMENT AND STRAIN AT DEPTH (PART-A)             *****
        # *****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
        # ********************************************************************
        #
        # ***** INPUT
        # *****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM
        # *****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
        # ***** OUTPUT
        # *****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES
        #
     #     COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
     #     COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,
        #     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ
        #      DATA F0,F2,PI2/0.D0,2.D0,6.283185307179586D0/

        U = num.zeros(12, float)
        DU = num.zeros(12, float)
        DUA = num.zeros(12, float)
        F0 = 0.0
        F2 = 2.0
        PI2 = 6.283185307179586

        ALP1 = self.ALP1
        ALP2 = self.ALP2
        ALP3 = self.ALP3
        ALP4 = self.ALP4
        ALP5 = self.ALP4
        SD = self.SD
        CD = self.CD
        SDSD = self.SDSD
        CDCD = self.CDCD
        SDCD = self.SDCD
        S2D = self.S2D
        C2D = self.C2D
        XI2 = self.XI2
        ET2 = self.ET2
        Q2 = self.Q2
        R = self.R
        R2 = self.R2
        R3 = self.R3
        R5 = self.R5
        Y = self.Y
        D = self.D
        TT = self.TT
        ALX = self.ALX
        ALE = self.ALE
        X11 = self.X11
        Y11 = self.Y11
        X32 = self.X32
        Y32 = self.Y32
        EY = self.EY
        EZ = self.EZ
        FY = self.FY
        FZ = self.FZ
        GY = self.GY
        GZ = self.GZ
        HY = self.HY
        HZ = self.HZ

# -----
        for I in range(0, 12):
            U[I] = F0
        XY = XI*Y11
        QX = Q * X11
        QY = Q * Y11
# ======================================
# =====  STRIKE-SLIP CONTRIBUTION  =====
# ======================================
        if DISL1 != F0:
            DU[0] = old_div(TT, F2) + ALP2*XI*QY
            DU[1] = old_div(ALP2*Q, R)
            DU[2] = ALP1*ALE - ALP2*Q*QY
            DU[3] = -ALP1*QY - ALP2*XI2*Q*Y32
            DU[4] = old_div(-ALP2*XI*Q, R3)
            DU[5] = ALP1*XY + ALP2*XI*Q2*Y32
            DU[6] = ALP1*XY*SD + ALP2*XI*FY+old_div(D, F2)*X11
            DU[7] = ALP2*EY
            DU[8] = ALP1*(old_div(CD, R)+QY*SD) - ALP2*Q*FY
            DU[9] = ALP1*XY*CD + ALP2*XI*FZ+old_div(Y, F2)*X11
            DU[10] = ALP2*EZ
            DU[11] = -ALP1*(old_div(SD, R)-QY*CD) - ALP2*Q*FZ
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL1, PI2)*DU[I]
# ======================================
# =====    DIP-SLIP CONTRIBUTION   =====
# ======================================
        if DISL2 != F0:
            DU[0] = old_div(ALP2*Q, R)
            DU[1] = old_div(TT, F2) + ALP2*ET*QX
            DU[2] = ALP1*ALX - ALP2*Q*QX
            DU[3] = old_div(-ALP2*XI*Q, R3)
            DU[4] = old_div(-QY, F2) - old_div(ALP2*ET*Q, R3)
            DU[5] = old_div(ALP1, R) + old_div(ALP2*Q2, R3)
            DU[6] = ALP2*EY
            DU[7] = ALP1*D*X11+old_div(XY, F2)*SD + ALP2*ET*GY
            DU[8] = ALP1*Y*X11 - ALP2*Q*GY
            DU[9] = ALP2*EZ
            DU[10] = ALP1*Y*X11+old_div(XY, F2)*CD + ALP2*ET*GZ
            DU[11] = -ALP1*D*X11 - ALP2*Q*GZ
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL2, PI2)*DU[I]
# ========================================
# =====  TENSILE-FAULT CONTRIBUTION  =====
# ========================================
        if DISL3 != F0:
            DU[0] = -ALP1*ALE - ALP2*Q*QY
            DU[1] = -ALP1*ALX - ALP2*Q*QX
            DU[2] = old_div(TT, F2) - ALP2*(ET*QX+XI*QY)
            DU[3] = -ALP1*XY + ALP2*XI*Q2*Y32
            DU[4] = old_div(-ALP1, R) + old_div(ALP2*Q2, R3)
            DU[5] = -ALP1*QY - ALP2*Q*Q2*Y32
            DU[6] = -ALP1*(old_div(CD, R)+QY*SD) - ALP2*Q*FY
            DU[7] = -ALP1*Y*X11 - ALP2*Q*GY
            DU[8] = ALP1*(D*X11+XY*SD) + ALP2*Q*HY
            DU[9] = ALP1*(old_div(SD, R)-QY*CD) - ALP2*Q*FZ
            DU[10] = ALP1*D*X11 - ALP2*Q*GZ
            DU[11] = ALP1*(Y*X11+XY*CD) + ALP2*Q*HZ
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL3, PI2)*DU[I]
        # for I in range (0,12):
            # DUA[I]=U[I]
        self.DUA = U

    def UB(self, XI, ET, Q, DISL1, DISL2, DISL3):
        from math import sqrt, atan, log
 # SUBROUTINE  UB(XI,ET,Q,DISL1,DISL2,DISL3,U)
 #     IMPLICIT REAL*8 (A-H,O-Z)
 #     DIMENSION U(12),DU(12)

        """********************************************************************
      C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****
      C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
      C********************************************************************
      C
      C***** INPUT
      C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM
      C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
      C***** OUTPUT
      C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES """

#      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
#      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,
#     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ
#      DATA  F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/

        DUB = num.zeros(12, float)
        DU = num.zeros(12, float)
        U = num.zeros(12, float)

        F0 = 0.0
        F1 = 1.0
        F2 = 2.0
        PI2 = 6.283185307179586

        ALP1 = self.ALP1
        ALP2 = self.ALP2
        ALP3 = self.ALP3
        ALP4 = self.ALP4
        ALP5 = self.ALP4
        SD = self.SD
        CD = self.CD
        SDSD = self.SDSD
        CDCD = self.CDCD
        SDCD = self.SDCD
        S2D = self.S2D
        C2D = self.C2D
        XI2 = self.XI2
        ET2 = self.ET2
        Q2 = self.Q2
        R = self.R
        R2 = self.R2
        R3 = self.R3
        R5 = self.R5
        Y = self.Y
        D = self.D
        TT = self.TT
        ALX = self.ALX
        ALE = self.ALE
        X11 = self.X11
        Y11 = self.Y11
        X32 = self.X32
        Y32 = self.Y32
        EY = self.EY
        EZ = self.EZ
        FY = self.FY
        FZ = self.FZ
        GY = self.GY
        GZ = self.GZ
        HY = self.HY
        HZ = self.HZ

        RD = R+D
        D11 = old_div(F1, (R*RD))
        AJ2 = old_div(XI*Y, RD)*D11
        AJ5 = -(D+old_div(Y*Y, RD))*D11
        if CD != F0:
            if XI == F0:
                AI4 = F0
            else:
                X = sqrt(XI2+Q2)
                AI4 = (old_div(F1, CDCD)*(old_div(XI, RD)*SDCD
                                          + F2*atan(old_div((ET*(X+Q*CD)+X*(R+X)*SD), (XI*(R+X)*CD)))))
            AI3 = old_div((old_div(Y*CD, RD)-ALE+SD*log(RD)), CDCD)
            AK1 = old_div(XI*(D11-Y11*SD), CD)
            AK3 = old_div((Q*Y11-Y*D11), CD)
            AJ3 = old_div((AK1-AJ2*SD), CD)
            AJ6 = old_div((AK3-AJ5*SD), CD)
        else:
            RD2 = RD*RD
            AI3 = old_div((old_div(ET, RD)+old_div(Y*Q, RD2)-ALE), F2)
            AI4 = old_div(old_div(XI*Y, RD2), F2)
            AK1 = old_div(XI*Q, RD)*D11
            AK3 = old_div(SD, RD)*(XI2*D11-F1)
            AJ3 = old_div(-XI, RD2)*(Q2*D11-old_div(F1, F2))
            AJ6 = old_div(-Y, RD2)*(XI2*D11-old_div(F1, F2))
# -----
        XY = XI*Y11
        AI1 = old_div(-XI, RD)*CD-AI4*SD
        AI2 = log(RD)+AI3*SD
        AK2 = old_div(F1, R)+AK3*SD
        AK4 = XY*CD-AK1*SD
        AJ1 = AJ5*CD-AJ6*SD
        AJ4 = -XY-AJ2*CD+AJ3*SD
# =====
        for I in range(0, 12):
            U[I] = F0
            QX = Q*X11
            QY = Q*Y11
# ======================================
# =====  STRIKE-SLIP CONTRIBUTION  =====
# ======================================
        if DISL1 != F0:
            DU[0] = -XI*QY-TT - ALP3*AI1*SD
            DU[1] = old_div(-Q, R) + old_div(ALP3*Y, RD)*SD
            DU[2] = Q*QY - ALP3*AI2*SD
            DU[3] = XI2*Q*Y32 - ALP3*AJ1*SD
            DU[4] = old_div(XI*Q, R3) - ALP3*AJ2*SD
            DU[5] = -XI*Q2*Y32 - ALP3*AJ3*SD
            DU[6] = -XI*FY-D*X11 + ALP3*(XY+AJ4)*SD
            DU[7] = -EY + ALP3*(old_div(F1, R)+AJ5)*SD
            DU[8] = Q*FY - ALP3*(QY-AJ6)*SD
            DU[9] = -XI*FZ-Y*X11 + ALP3*AK1*SD
            DU[10] = -EZ + ALP3*Y*D11*SD
            DU[11] = Q*FZ + ALP3*AK2*SD
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL1, PI2)*DU[I]
# ======================================
# =====    DIP-SLIP CONTRIBUTION   =====
# ======================================
        if DISL2 != F0:
            DU[0] = old_div(-Q, R) + ALP3*AI3*SDCD
            DU[1] = -ET*QX-TT - old_div(ALP3*XI, RD)*SDCD
            DU[2] = Q*QX + ALP3*AI4*SDCD
            DU[3] = old_div(XI*Q, R3) + ALP3*AJ4*SDCD
            DU[4] = old_div(ET*Q, R3)+QY + ALP3*AJ5*SDCD
            DU[5] = old_div(-Q2, R3) + ALP3*AJ6*SDCD
            DU[6] = -EY + ALP3*AJ1*SDCD
            DU[7] = -ET*GY-XY*SD + ALP3*AJ2*SDCD
            DU[8] = Q*GY + ALP3*AJ3*SDCD
            DU[9] = -EZ - ALP3*AK3*SDCD
            DU[10] = -ET*GZ-XY*CD - ALP3*XI*D11*SDCD
            DU[11] = Q*GZ - ALP3*AK4*SDCD
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL2, PI2)*DU[I]
# ========================================
# =====  TENSILE-FAULT CONTRIBUTION  =====
# ========================================
        if DISL3 != F0:
            DU[0] = Q*QY - ALP3*AI3*SDSD
            DU[1] = Q*QX + old_div(ALP3*XI, RD)*SDSD
            DU[2] = ET*QX+XI*QY-TT - ALP3*AI4*SDSD
            DU[3] = -XI*Q2*Y32 - ALP3*AJ4*SDSD
            DU[4] = old_div(-Q2, R3) - ALP3*AJ5*SDSD
            DU[5] = Q*Q2*Y32 - ALP3*AJ6*SDSD
            DU[6] = Q*FY - ALP3*AJ1*SDSD
            DU[7] = Q*GY - ALP3*AJ2*SDSD
            DU[8] = -Q*HY - ALP3*AJ3*SDSD
            DU[9] = Q*FZ + ALP3*AK3*SDSD
            DU[10] = Q*GZ + ALP3*XI*D11*SDSD
            DU[11] = -Q*HZ + ALP3*AK4*SDSD
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL3, PI2)*DU[I]

        # for I in range(0,12):
                # DUB[I]=U[I]
        self.DUB = U

    def UC(self, XI, ET, Q, Z, DISL1, DISL2, DISL3):
     #     SUBROUTINE  UC(XI,ET,Q,Z,DISL1,DISL2,DISL3,U)
     #     IMPLICIT REAL*8 (A-H,O-Z)
      #    DIMENSION U(12),DU(12)
        """********************************************************************
        C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-C)             *****
        C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
        C********************************************************************
        C
        C***** INPUT
        C*****   XI,ET,Q,Z   : STATION COORDINATES IN FAULT SYSTEM
        C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
        C***** OUTPUT
        C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES  """
#     COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
#      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,
 #    *           EY,EZ,FY,FZ,GY,GZ,HY,HZ
 #     DATA F0,F1,F2,F3,PI2/0.D0,1.D0,2.D0,3.D0,6.283185307179586D0/

        DUC = num.zeros(12, float)
        DU = num.zeros(12, float)
        U = num.zeros(12, float)

        F0 = 0.0
        F1 = 1.0
        F2 = 2.0
        F3 = 3.0
        PI2 = 6.283185307179586

        ALP1 = self.ALP1
        ALP2 = self.ALP2
        ALP3 = self.ALP3
        ALP4 = self.ALP4
        ALP5 = self.ALP4
        SD = self.SD
        CD = self.CD
        SDSD = self.SDSD
        CDCD = self.CDCD
        SDCD = self.SDCD
        S2D = self.S2D
        C2D = self.C2D
        XI2 = self.XI2
        ET2 = self.ET2
        Q2 = self.Q2
        R = self.R
        R2 = self.R2
        R3 = self.R3
        R5 = self.R5
        Y = self.Y
        D = self.D
        TT = self.TT
        ALX = self.ALX
        ALE = self.ALE
        X11 = self.X11
        Y11 = self.Y11
        X32 = self.X32
        Y32 = self.Y32
        EY = self.EY
        EZ = self.EZ
        FY = self.FY
        FZ = self.FZ
        GY = self.GY
        GZ = self.GZ
        HY = self.HY
        HZ = self.HZ
# -----
        C = D+Z
        X53 = old_div((8.0*R2+9.0*R*XI+F3*XI2)*X11*X11*X11, R2)
        Y53 = old_div((8.0*R2+9.0*R*ET+F3*ET2)*Y11*Y11*Y11, R2)
        H = Q*CD-Z
        Z32 = old_div(SD, R3)-H*Y32
        Z53 = old_div(F3*SD, R5)-H*Y53
        Y0 = Y11-XI2*Y32
        Z0 = Z32-XI2*Z53
        PPY = old_div(CD, R3)+Q*Y32*SD
        PPZ = old_div(SD, R3)-Q*Y32*CD
        QQ = Z*Y32+Z32+Z0
        QQY = old_div(F3*C*D, R5)-QQ*SD
        QQZ = old_div(F3*C*Y, R5)-QQ*CD+Q*Y32
        XY = XI*Y11
        QX = Q*X11
        QY = Q*Y11
        QR = old_div(F3*Q, R5)
        CQX = C*Q*X53
        CDR = old_div((C+D), R3)
        YY0 = old_div(Y, R3)-Y0*CD
# =====
        for I in range(1, 12):
            U[I] = F0
# ======================================
# =====  STRIKE-SLIP CONTRIBUTION  =====
# ======================================
        if DISL1 != F0:
            DU[0] = ALP4*XY*CD - ALP5*XI*Q*Z32
            DU[1] = ALP4*(old_div(CD, R)+F2*QY*SD) - old_div(ALP5*C*Q, R3)
            DU[2] = ALP4*QY*CD - ALP5*(old_div(C*ET, R3)-Z*Y11+XI2*Z32)
            DU[3] = ALP4*Y0*CD - ALP5*Q*Z0
            DU[4] = -ALP4*XI*(old_div(CD, R3)+F2*Q*Y32*SD) + ALP5*C*XI*QR
            DU[5] = -ALP4*XI*Q*Y32*CD + ALP5*XI*(old_div(F3*C*ET, R5)-QQ)
            DU[6] = -ALP4*XI*PPY*CD - ALP5*XI*QQY
            DU[7] = (ALP4*F2*(old_div(D, R3)-Y0*SD)*SD-old_div(Y, R3)*CD
                     - ALP5*(CDR*SD-old_div(ET, R3)-C*Y*QR))
            DU[8] = old_div(-ALP4*Q, R3)+YY0*SD + ALP5 * \
                (CDR*CD+C*D*QR-(Y0*CD+Q*Z0)*SD)
            DU[9] = ALP4*XI*PPZ*CD - ALP5*XI*QQZ
            DU[10] = ALP4*F2*(old_div(Y, R3)-Y0*CD)*SD + \
                old_div(D, R3)*CD - ALP5*(CDR*CD+C*D*QR)
            DU[11] = YY0*CD - ALP5*(CDR*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD)
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL1, PI2)*DU[I]

# ======================================
# =====    DIP-SLIP CONTRIBUTION   =====
# ======================================
        if DISL2 != F0:
            DU[0] = old_div(ALP4*CD, R) - QY*SD - old_div(ALP5*C*Q, R3)
            DU[1] = ALP4*Y*X11 - ALP5*C*ET*Q*X32
            DU[2] = -D*X11-XY*SD - ALP5*C*(X11-Q2*X32)
            DU[3] = old_div(-ALP4*XI, R3)*CD + ALP5*C*XI*QR + XI*Q*Y32*SD
            DU[4] = old_div(-ALP4*Y, R3) + ALP5*C*ET*QR
            DU[5] = old_div(D, R3)-Y0*SD + old_div(ALP5 *
                                                   C, R3)*(F1-old_div(F3*Q2, R2))
            DU[6] = old_div(-ALP4*ET, R3)+Y0*SDSD - ALP5*(CDR*SD-C*Y*QR)
            DU[7] = ALP4*(X11-Y*Y*X32) - ALP5*C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)
            DU[8] = XI*PPY*SD+Y*D*X32 + ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53)
            DU[9] = old_div(-Q, R3)+Y0*SDCD - ALP5*(CDR*CD+C*D*QR)
            DU[10] = ALP4*Y*D*X32 - ALP5*C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)
            DU[11] = -XI*PPZ*SD+X11-D*D*X32-ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53)
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL2, PI2)*DU[I]
# ========================================
# =====  TENSILE-FAULT CONTRIBUTION  =====
# ========================================
        if DISL3 != F0:
            DU[0] = -ALP4*(old_div(SD, R)+QY*CD) - ALP5*(Z*Y11-Q2*Z32)
            DU[1] = ALP4*F2*XY*SD+D*X11 - ALP5*C*(X11-Q2*X32)
            DU[2] = ALP4*(Y*X11+XY*CD) + ALP5*Q*(C*ET*X32+XI*Z32)
            DU[3] = old_div(ALP4*XI, R3)*SD+XI*Q*Y32*CD+ALP5 * \
                XI*(old_div(F3*C*ET, R5)-F2*Z32-Z0)
            DU[4] = ALP4*F2*Y0*SD - \
                old_div(D, R3) + old_div(ALP5*C, R3)*(F1-old_div(F3*Q2, R2))
            DU[5] = -ALP4*YY0 - ALP5*(C*ET*QR-Q*Z0)
            DU[6] = ALP4*(old_div(Q, R3)+Y0*SDCD) + ALP5 * \
                (old_div(Z, R3)*CD+C*D*QR-Q*Z0*SD)
            DU[7] = (-ALP4*F2*XI*PPY*SD-Y*D*X32
                     + ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53))
            DU[8] = (-ALP4*(XI*PPY*CD-X11+Y*Y*X32)
                     + ALP5*(C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)+XI*QQY))
            DU[9] = old_div(-ET, R3)+Y0*CDCD - ALP5 * \
                (old_div(Z, R3)*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD)
            DU[10] = (ALP4*F2*XI*PPZ*SD-X11+D*D*X32
                      - ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53))
            DU[11] = (ALP4*(XI*PPZ*CD+Y*D*X32)
                      + ALP5*(C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)+XI*QQZ))
            for I in range(0, 12):
                U[I] = U[I]+old_div(DISL3, PI2)*DU[I]

        # for I in range(0,12):
            # DUC[I]=U[I]
        self.DUC = U

    def DCC0N0(self, ALPHA, DIP):
        from math import sin, cos
#      SUBROUTINE  DCCON0(ALPHA,DIP)
#      IMPLICIT REAL*8 (A-H,O-Z)
#
        """*******************************************************************
      C*****   CALCULATE MEDIUM CONSTANTS AND FAULT-DIP CONSTANTS    *****
      C*******************************************************************
      C
      C***** INPUT
      C*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)
      C*****   DIP   : DIP-ANGLE (DEGREE)
      C### CAUTION ### IF COS(DIP) IS SUFFICIENTLY SMALL, IT IS SET TO ZERO """

#      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
 #     DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/
#      DATA EPS/1.D-6/
        F0 = 0.0
        F1 = 1.0
        F2 = 2.0
        PI2 = 6.283185307179586
        EPS = 1.0e-6

# -----
        ALP1 = old_div((F1-ALPHA), F2)
        ALP2 = old_div(ALPHA, F2)
        ALP3 = old_div((F1-ALPHA), ALPHA)
        ALP4 = F1-ALPHA
        ALP5 = ALPHA
        # print ALP1
# -----
        P18 = PI2/360.0
        SD = sin(DIP*P18)
        CD = cos(DIP*P18)
        if abs(CD) < EPS:
            CD = F0
            if SD > F0:
                SD = F1
            if SD < F0:
                SD = -F1
        SDSD = SD*SD
        CDCD = CD*CD
        SDCD = SD*CD
        S2D = F2*SDCD
        C2D = CDCD-SDSD

        self.ALP1 = ALP1
        self.ALP2 = ALP2
        self.ALP3 = ALP3
        self.ALP4 = ALP4
        self.ALP5 = ALP5
        self.SD = SD
        self.CD = CD
        self.SDSD = SDSD
        self.CDCD = CDCD
        self.SDCD = SDCD
        self.S2D = S2D
        self.C2D = C2D

    def DCCON1(self, X, Y, D):
        from math import sqrt
#     SUBROUTINE  DCCON1(X,Y,D)
 #     IMPLICIT REAL*8 (A-H,O-Z)

        """**********************************************************************
      C*****   CALCULATE STATION GEOMETRY CONSTANTS FOR POINT SOURCE    *****
      C**********************************************************************
      C
      C***** INPUT
      C*****   X,Y,D : STATION COORDINATES IN FAULT SYSTEM
      C### CAUTION ### IF X,Y,D ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZERO  """

#      COMMON /C0/DUMMY(5),SD,CD,dumm(5)
#      COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,
#     *           UY,VY,WY,UZ,VZ,WZ
        F0 = 0.0
        F1 = 1.0
        F3 = 3.0
        F5 = 5.0
        EPS = 1.0e-6

        SD = self.SD
        CD = self.CD
# -----
        if abs(X) < EPS:
            X = F0
        if abs(Y) < EPS:
            Y = F0
        if abs(D) < EPS:
            D = F0
        P = Y*CD+D*SD
        Q = Y*SD-D*CD
        S = P*SD+Q*CD
        T = P*CD-Q*SD
        XY = X*Y
        X2 = X*X
        Y2 = Y*Y
        D2 = D*D
        R2 = X2+Y2+D2
        R = sqrt(R2)
        if R == F0:
            return
        R3 = R * R2
        R5 = R3*R2
        R7 = R5*R2
# -----
        A3 = F1-old_div(F3*X2, R2)
        A5 = F1-old_div(F5*X2, R2)
        B3 = F1-old_div(F3*Y2, R2)
        C3 = F1-old_div(F3*D2, R2)
# -----
        QR = old_div(F3*Q, R5)
        QRX = old_div(F5*QR*X, R2)
# -----
        UY = SD-old_div(F5*Y*Q, R2)
        UZ = CD+old_div(F5*D*Q, R2)
        VY = S - old_div(F5*Y*P*Q, R2)
        VZ = T + old_div(F5*D*P*Q, R2)
        WY = UY+SD
        WZ = UZ+CD

        self.P = P
        self.Q = Q
        self.S = S
        self.T = T
        self.XY = XY
        self.X2 = X2
        self.Y2 = Y2
        self.D2 = D2
        self.R2 = R2
        self.R = R
        self.R3 = R3
        self.R5 = R5
        self.R7 = R7
        self.A3 = A3
        self.A5 = A5
        self.B3 = B3
        self.C3 = C3
        self.QR = QR
        self.QRX = QRX
        self.UY = UY
        self.UZ = UZ
        self.VZ = VZ
        self.VY = VY
        self.WY = WY
        self.WZ = WZ

    def DCCON2(self, XI, ET, Q, SD, CD, KXI, KET):
        from math import sqrt, atan, log
#      SUBROUTINE  DCCON2(XI,ET,Q,SD,CD,KXI,KET)
#      IMPLICIT REAL*8 (A-H,O-Z)

        """**********************************************************************
      C*****   CALCULATE STATION GEOMETRY CONSTANTS FOR FINITE SOURCE   *****
      C**********************************************************************
      C
      C***** INPUT
      C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM
      C*****   SD,CD   : SIN, COS OF DIP-ANGLE
      C*****   KXI,KET : KXI=1, KET=1 MEANS R+XI<EPS, R+ET<EPS, RESPECTIVELY
      C
      C### CAUTION ### IF XI,ET,Q ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZER0"""

        F0 = 0.0
        F1 = 1.0
        F2 = 2.0
        EPS = 1.0e-6
# -----
        if abs(XI) < EPS:
            XI = F0
        if abs(ET) < EPS:
            ET = F0
        if abs(Q) < EPS:
            Q = F0
        XI2 = XI*XI
        ET2 = ET*ET
        Q2 = Q*Q
        R2 = XI2+ET2+Q2
        R = sqrt(R2)
        if R == F0:
            return
        R3 = R * R2
        R5 = R3*R2
        Y = ET*CD+Q*SD
        D = ET*SD-Q*CD
# -----
        if Q == F0:
            TT = F0
        else:
            TT = atan(old_div(XI*ET, (Q*R)))
# -----
        if KXI == 1:
            ALX = -log(R-XI)
            X11 = F0
            X32 = F0
        else:
            RXI = R+XI
            ALX = log(RXI)
            X11 = old_div(F1, (R*RXI))
            X32 = old_div((R+RXI)*X11*X11, R)
# -----
        if KET == 1:
            ALE = -log(R-ET)
            Y11 = F0
            Y32 = F0
        else:
            RET = R+ET
            ALE = log(RET)
            Y11 = old_div(F1, (R*RET))
            Y32 = old_div((R+RET)*Y11*Y11, R)
# -----
        EY = old_div(SD, R)-old_div(Y*Q, R3)
        EZ = old_div(CD, R)+old_div(D*Q, R3)
        FY = old_div(D, R3)+XI2*Y32*SD
        FZ = old_div(Y, R3)+XI2*Y32*CD
        GY = F2*X11*SD-Y*Q*X32
        GZ = F2*X11*CD+D*Q*X32
        HY = D*Q*X32+XI*Q*Y32*SD
        HZ = Y*Q*X32+XI*Q*Y32*CD

        self.XI2 = XI2
        self.Q2 = Q2
        self.R = R
        self.R2 = R2
        self.R3 = R3
        self.R5 = R5
        self.Y = Y
        self.D = D
        self.TT = TT
        self.ALX = ALX
        self.ALE = ALE
        self.X11 = X11
        self.Y11 = Y11
        self.X32 = X32
        self.Y32 = Y32
        self.EY = EY
        self.EZ = EZ
        self.FY = FY
        self.FZ = FZ
        self.GY = GY
        self.GZ = GZ
        self.HY = HY
        self.HZ = HZ
        self.ET2 = ET2
