We have calculated the auxiliary fault plane solution using the input of main fault solution in perl. Here, we do the same in Fortran. And we will also plot to solution to visualize the results.
Fortran Code to get the auxiliary fault plane solutions:
program auxiliary_fault_plane
! Program to calculate the strike, dip and rake of the auxiliary fault plane solutions
! given the strike, dip and rake of the fault plane.
! Authors: Utpal Kumar, Li Zhao
implicit none
real (kind=8):: ph1,del1,lb1,ph2,del2,lb2,sinlb2,coslb2,sinph1_ph2,cosph1_ph2,ph1_ph2,phI,delI,lbI
real, parameter :: pi = 2*asin(1.0)
integer :: n,lbsgn
print *, 'Enter strike, dip and rake respectively of the fault plane (in degrees) (e.g., 40,50,60):'
read *, ph1,del1,lb1
!ph1=40
!del1=80
!lb1=200
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
!Adaptation to include the negative values of rake
lbsgn=1
if (lb1 < 0) then
lbsgn=-1
end if
! Calculation of Strike, sip and rake of the auxiliary fault plane
del2=acos(sin(lb1)*sin(del1)) !dip of auxialiary fault plane
sinlb2=cos(del1)/sin(del2)
coslb2=-(sin(del1)*cos(lb1)/sin(del2))
lb2=acos(coslb2) !rake of auxilairy fault plane
sinph1_ph2=cos(lb1)/sin(del2)
cosph1_ph2=-1/(tan(del1)*tan(del2))
ph1_ph2=acos(cosph1_ph2)
!Checking for the quadrant of the strike angle
if (sinph1_ph2 >= 0 .and. cosph1_ph2 >=0) then
ph1_ph2=ph1_ph2
else if (sinph1_ph2 > 0 .and. cosph1_ph2 < 0) then
ph1_ph2=ph1_ph2
else if (sinph1_ph2 < 0 .and. cosph1_ph2 < 0) then
ph1_ph2=-ph1_ph2
else if (sinph1_ph2 < 0 .and. cosph1_ph2 > 0) then
ph1_ph2=-ph1_ph2
end if
ph2=ph1-ph1_ph2 !strike of auxialiary fault plane
! For dip > 90 degrees and less than 180 degrees
if (del2 > pi/2 .and. del2 < pi) then
ph2= pi + ph2
del2= pi - del2
lb2= 2*pi - lb2
end if
if (lbsgn < 0) then
lb2 = -(2*pi - lb2)
end if
!Adaptation to give the strike value in the range of 0 to 360 degrees
if (ph2 > 2*pi) then
ph2 = ph2 - 2*pi
end if
101 format("The strike of the auxiliary plane is: ",f9.4, " degrees")
102 format("The dip of the auxiliary plane is: ",f9.4, " degrees")
103 format("The rake of the auxiliary plane is: ",f9.4, " degrees")
print 101, ph2*180/pi
print 102, del2*180/pi
print 103, lb2*180/pi
open(unit=10,file='plt.dat')
10 format("25 25 0 ", f8.2, f8.2, f8.2, f8.2, f8.2, f8.2," MainFault: ",f7.2,"/",f7.2,"/",f7.2, & !continuation of line
" AuxFault: ",f7.2,"/",f7.2,"/",f7.2)
write(10,10) phI,delI,lbI,ph2*180/pi, del2*180/pi, lb2*180/pi, phI,delI,lbI,ph2*180/pi, del2*180/pi, lb2*180/pi !writing in the file
end program auxiliary_fault_plane
Then we save this in the file called “auxiliary_fault_plane.f”. This script basically aim for calculation of the fault plane solution. For visualizing the results, we can execute the program and plot the results using the psmeca command of GMT.
The bash script for compiling and executing the above program and then plotting it is following:
#!/bin/bash
gfortran -ffree-form auxiliary_fault_plane.f -o auxiliary_fault_plane
./auxiliary_fault_plane
psbasemap -Bwsne -R0/50/0/50 -Jm1.0 -Xc -Yc -P -K >output.ps
echo "24.6 30 20 0 2 MC Focal Mechanism Plot"| pstext -Jm -R -K -O -P -N >>output.ps
psmeca plt.dat -Jm -R -Sa0.40/14/-0.2i -Gbrown -Fa0.2c/cd -Egray -P -O -N -L -V >>output.ps
gs output.ps
ps2pdf output.ps
The P axis is plotted using circle and T axis using the diamond symbol (-Fa0.2c/cd).
-Utpal Kumar (IES, Academia Sinica)
Institute of Earth Sciences Coders 