module solve_for

	implicit none 

	contains

	real function solve_for_x1(rho,hh,z1,z2,s)
	
		!----------------------------------------------------------
		! Function to solve a quartic polynomial equation by the 
		! method of bisection. The solution is the distance between
		! the shot and the first bounce of the multiple at the WB.
		! The equation is:
    	! (rho^2-1)x1^4-2(rho^2-1)hh*x1^3+[(rho^2-1)hh^2+rho^2*
		!	z2^2-z1^2]x1^2+2hh*z1^2*x1-z1^2*hh^2=0
		! Ax1^4+Bx1^3+Cx1^2+Dx1+E=0
		!----------------------------------------------------------

		! argument variables
		real, intent(in)	:: rho ! v2/v1
		real, intent(in)	:: hh  ! half-offset
		real, intent(in)	:: z1	! depth of water bottom
		real, intent(in)	:: z2	! depth of sediement reflector
		integer,intent(in)	:: s	! =1 for multiples symmetric wrt the cmp
									! /=1 for other multiples
                                                                                     
		! internal variables
		integer	:: k
		real	:: A, B, C, D, E	! function coefficients
		real	:: x1, x2, x3		! trial values
		real	:: f1, f2, f3		! trial function values

		! compute function coefficients
		A=s**2*(rho**2-1)
		B=-2*s*(rho**2-1)*hh
		C=(rho**2-1)*hh**2+rho**2*z2**2-s**2*z1**2
		D=s**2*hh*z1**2
		E=-z1**2hh**2
                                                                                     
		! set initial guesses
		x1=0;
		x2=hh;
		x3=(x1+x2)/2.;
                                                                                     
		! loop until solution found
		k=1;
		do
			f1=A*x1**4+B*x1**3+C*x1**2+D*x1+E;
			f2=A*x2**4+B*x2**3+C*x2**2+D*x2+E;
			f3=A*x3**4+B*x3**3+C*x3**2+D*x3+E;
			if (abs(f3)<0.01) then
				solve_for_x1=x3
				return
        	end if
        	if (f1*f3<0) then
				x2=x3;
			else
				x1=x3;
			end if
			x3=(x1+x2)/2;

			if (k>100) then
				write(0,*)'ERROR: could not solve for x1 in 1000 iters'
				return;
			end if
			k=k+1;
		end do
		return
	end function solve_for_x1

end module solve_for
