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):'

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