program two_layers

c* A program to compute theoretical apparent splitting parameters
c*	for two pairs of spitting parameters.

c* Based on the method proposed by Silver and Savage, 1994.
c* The output file is two_layers.fitted. 

c* I also made a GMT program (plot_2layers.gmt) to plot the results
c
c* For questions, please contact:
c*  Stephen S. Gao, Missouri University of Science and Technology.
c* (EmailL sgao@mst.edu;  http://www.mst.edu/~sgao) 

c* See Figure 9 of the following paper (Liu and Gao, 2013, BSSA) 
c*	for more info about two-layer modeling:
c* http://web.mst.edu/~sgao/publications/BSSA2013.pdf


	write(6,1)
1	format('Enter phi and dt of the lower layer: ',$)
	read(5,*) phi1, dt1

	write(6,2)
2	format('Enter phi and dt of the upper layer: ',$)
	read(5,*) phi2, dt2

	write(6,3)
3	format('Enter 1 if the range if 0 to 180,',
     &		' and enter 2 if the range is -90 to 90:',$)
	read(5,*) irange

	write(6,4)
4	format('Enter peak frequency (e.g., 0.1): ',$)
	read(5,*) freq

	call layer2_allaz(phi1, dt1, phi2, dt2, irange, freq)

	print*,'The output file is two_layers.fitted'	
	print*,'Column 1: baz, 2: fitted phi, 3: fitted dt'

	stop
	end

*-------------------------------

      subroutine layer2_allaz(phi1, dt1, phi2, dt2, irange,freq)

      character range*8
      complex ak

      pi=4.0*atan(1.0)

      rad= 180./pi

	if(irange.eq.1) range='000to180'
	if(irange.eq.2) range='m90top90'

	open(4,file='two_layers.fitted')

      t1=pi*freq*dt1
      t2=pi*freq*dt2

      do 10 i = 1,361
        phat=float(i-1)
        al1=2.*(phi1-phat)/rad
        al2=2.*(phi2-phat)/rad
        cc  =cos(t1)*sin(t2)*cos(al2) + cos(t2)*sin(t1)*cos(al1)
        cs  =cos(t1)*sin(t2)*sin(al2) + cos(t2)*sin(t1)*sin(al1)
        ap  =cos(t1)*cos(t2) - sin(t1)*sin(t2)*cos(al2-al1)
        app=                 - sin(t1)*sin(t2)*sin(al2-al1)
        al0=atan((app**2 + cs**2)/(app*ap + cs*cc))
        t0= atan(app/(cs*cos(al0) - cc*sin(al0)))
        ak=cmplx(ap,cc)/cmplx(cos(t0),sin(t0)*cos(al0))
        amp=cabs(ak)
        phase=atan2(aimag(ak),real(ak))
        phi0 = .5*al0*rad+phat
        dt0=t0/(pi*freq)
        phase=phase*rad
c     if dt0 is negative, make pos and add 90 to phi0
        if(dt0.lt.0.)then
          dt0=-dt0
          phi0=phi0+90.
        endif

        if(phi0.lt.-90.) phi0=phi0+180.  
        if(phi0.gt.360.) phi0=phi0-360. 
        if(phi0.gt.180.) phi0=phi0-180. 

	if(range.eq.'m90top90') then
        if(phi0.gt.90.) phi0=phi0-180. 
	endif

	if(range.eq.'000to180') then
        if(phi0.lt.00.) phi0=phi0+180. 
	endif

        write(4,110) phat, phi0,dt0,amp, phase, 
     &               phi1,dt1,phi2, dt2
110     format(2(1x,f7.2),f10.2,1x,6(f10.3,1x))
10      continue

	close(4)

	call writegmt(irange)

	print*,' '
	print*,'***Note: freq=',freq,' Hz. Results are freq dependent!'
	print*,' '

	end

*-------------------------------

* Create a GMT program to plot the results

        subroutine writegmt(irange)
        character cmd*100

	open(92,file='plot_2layers.gmt')
 
        write(92,921)
921     format('/bin/rm tmp.ps'/)


        if(irange.eq.2) then

	write(92,922) 
922    format('awk ',
     &  '''{print $1, $2}'' two_layers.fitted | \\'/,
     & ' psxy -R0/90/-90/90 -JX5/3 -K -W15/255/0/0 \\'/,
     & '-Ba10g10f5  -P \\'/,
     & '-X2 -Y6 >tmp.ps'/)
 

        write(92,2)
2       format('awk ',
     &  '''{print "45 100 15 0 0 6", ',
     &  '" Lower:", $6, $7, "; Upper: ",$8, $9}'' \\'/,
     & 	'two_layers.fitted | pstext -R -JX -O -K -N >>tmp.ps'/)


	write(92, 3)
3	format('pstext -R -J -K -O -N <<END>>tmp.ps',/
     &	       '45 -120 13 0 5 6 Back azimuth (modulo 90 deg.)',/
     &	       'END',/)


	write(92,923)
923     format('awk ',
     &  '''{print $1, $3}'' two_layers.fitted |  \\'/,
     &  ' psxy -R0/90/0/5 -JX5/3 -O -W15/255/0/0 \\'/,
     & '-Ba10g10f5/a1f0.5g1  -P \\'/,
     & ' -Y-4 >>tmp.ps'/)


	endif


        if(irange.eq.1) then


	write(92,9221) 
9221    format('awk ',
     &  '''{print $1, $2}'' two_layers.fitted | \\'/,
     & 'psxy -R0/90/0/180 -JX5/3 -K -W15/255/0/0 \\'/,
     & '-Ba10g10f5  -P \\'/,
     & '-X2 -Y6 >tmp.ps'/)
 

        write(92,222)
222       format('awk ',
     &  '''{print "45 190 15 0 0 6", ',
     &  '" Lower:", $6, $7, "; Upper: ",$8, $9}'' \\'/,
     & 	'two_layers.fitted | pstext -R -JX -O -K -N >>tmp.ps'/)

	write(92, 31)
31	format('pstext -R -J -K -O -N <<END>>tmp.ps',/
     &	       '45 -30 13 0 5 6 Back azimuth (modulo 90 deg.)',/
     &	       'END',/)



	write(92,9231)
9231     format('awk ',
     &  '''{print $1, $3}'' two_layers.fitted |  \\'/,
     &  ' psxy -R0/90/0/5 -JX5/3 -O -W15/255/0/0 \\'/,
     & '-Ba10g10f5/a1f0.5g1  -P \\'/,
     & ' -Y-4 >>tmp.ps'/)

	endif


	close(92)

	cmd='chmod a+x plot*gmt'
	call system(cmd)

	cmd='plot_2layers.gmt'
	call system(cmd)

	cmd='gv tmp.ps&'
	call system(cmd)

        end