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
        !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
        if (lb1 < 0) then
        end if
        ! Calculation of Strike, sip and rake of the auxiliary fault plane
        del2=acos(sin(lb1)*sin(del1))   !dip of auxialiary fault plane
        lb2=acos(coslb2)        !rake of auxilairy fault plane

        !Checking for the quadrant of the strike angle
        if (sinph1_ph2 >= 0 .and. cosph1_ph2 >=0) then
        else if (sinph1_ph2 > 0 .and. cosph1_ph2 < 0) then
        else if (sinph1_ph2 < 0 .and. cosph1_ph2 < 0) then
        else if (sinph1_ph2 < 0 .and. cosph1_ph2 > 0) then
        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
    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:

gfortran -ffree-form auxiliary_fault_plane.f -o 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).

Screenshot from 2016-12-05 14-34-26.png


-Utpal Kumar (IES, Academia Sinica)