Calculating Pressure and Tension axis given the fault plane (strike, dip and rake): Analytical solution

The fault plane is characterized by its normal vector and the direction of its motion is given by the slip vector. Here, we calculate the normal vector and the slip vector for a given fault geometry.

The orientation of the normal vector and the slip vector is given in a geographical coordinate system with x pointing in north direction, y pointing in west, and z pointing up (see Stein and Wysession (2003) pg. 218 for details).

program p_t_axis_calc
! Program to calculate the pressure and tension axis 
! given the strike, dip and rake of the fault plane.
! The coordinate system used is the geographical coordinate system
! x/1 - north, y/2 - west, z/3 - upwards
! Author: Utpal Kumar
        implicit none
        real (kind=8):: ph1,del1,lb1,phI,delI,lbI,tVecMag,pVecMag
    real, parameter :: pi = 2*asin(1.0)
    integer :: i,lbsgn
    real(kind=16), dimension(1:3) :: normVec, slipVec, nullVec, tVec, pVec, tcos, pcos
        print *, 'Enter strike, dip and rake respectively of the fault plane (in degrees) (e.g., 40,50,60):'
        read *, ph1,del1,lb1

        phI=ph1
    delI=del1
    lbI=lb1
        !Converting the strike, dip and rake in radians
        ph1=ph1*pi/180  !strike
        del1=del1*pi/180        !dip    
        lb1=lb1*pi/180          !rake   
        
        !Normal Vector Components
        normVec(1)=-sin(del1)*sin(ph1)
    normVec(2)=-sin(del1)*cos(ph1)
    normVec(3)=cos(del1)

    !Slip Vector Components
    slipVec(1)=cos(lb1)*cos(ph1) + sin(lb1)*cos(del1)*sin(ph1)
    slipVec(2)=-cos(lb1)*sin(ph1) + sin(lb1)*cos(del1)*cos(ph1)
    slipVec(3)=sin(lb1)*sin(del1)

    !Null Vector components
    nullVec(1)=normVec(2)*slipVec(3) - slipVec(2)*normVec(3)
    nullVec(2)=-(normvec(1)*slipVec(3) - slipvec(1)*normVec(3))
    nullVec(3)=normVec(1)*slipVec(2) - slipVec(1)*normVec(2)
    
    !Tension Vector and Pressure Vector
    do i=1,3
        tVec(i)=normVec(i)+slipVec(i)
        pVec(i)=normVec(i)+slipVec(i)
    end do

    !Magnitude of t and p Vectors
    tVecMag=sqrt(tVec(1)**2 + tVec(2)**2 + tVec(3)**2)
    pVecMag=sqrt(pVec(1)**2 + pVec(2)**2 + pVec(3)**2)

    !Direction cosines of the pressure and tension vectors w.r.t x,y,z axis
    do i=1,3 
    tcos(i) = tVec(i)/tVecMag
    pcos(i) = pVec(i)/pVecMag
    end do

    write(*,20) "Normal Vec", "  Slip Vec", " Null Vec", " T- Vec", "  T-cosines", " P-Vec", " P -cosines"
    do i=1,3
    write(*,10) normVec(i), slipVec(i), nullVec(i), tVec(i), tcos(i), pVec(i), pcos(i)
    end do
    10 format(f12.4, 5X,f12.4, 5X,f12.4, 5X,f12.4, 5X,f12.4, 5X,f12.4, 5X,f12.4)
    20 format(A12,5X,A12,5X,A12,5X,A12,5X,A12,5X,A12,5X,A12)
        
end program p_t_axis_calc

For compiling and running the above program

gfortran -ffree-form p_t_axis_calc.f -o p-t_axis_calc
./p-t_axis_calc

Screenshot from 2016-12-12 18-32-08.png

 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Powered by WordPress.com.

Up ↑

%d bloggers like this: