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

        !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

    !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)

    !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
    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

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