module processamento

use variaveis

contains

!------------------------------------------------

subroutine funcoes(x,F)

	real*8,dimension(2),intent(inout) :: x, F

	F(1) = K1 - ( (2.0d0*x(1) + x(2) )/ &
	              (1.0d0 - 2.0d0*x(1) - 2.0d0*x(2) ) )**2 *	&
			    x(1) / (1.0d0 + x(1) + x(2))
	F(2) = K2 - ( 2.0d0*x(1) + x(2) )/(1.0d0 + x(1) + x(2)) *&
	            (2.0d0*x(2)/(1.0d0-2.0d0*x(1)-2.0d0*x(2)))**2

end subroutine funcoes

!------------------------------------------------

subroutine jacobiano

	real*8,dimension(:),allocatable :: xe, xw, Fe, Fw

 	integer :: ii, jj

	allocate ( xe(2), xw(2), Fe(2), Fw(2) )

	do jj = 1, 2

		xw = x; xe = x

		xw(jj) = x(jj) - dist*x(jj)
		xe(jj) = x(jj) + dist*x(jj)

		do ii = 1, 2

			call funcoes(xw,Fw)
			call funcoes(xe,Fe)

			J(ii,jj) = ( Fe(ii)-Fw(ii) )/( xe(jj) - xw(jj) )

		end do

	end do

	deallocate ( xe, xw, Fe, Fw ) 

end subroutine jacobiano

!------------------------------------------------

subroutine gauss_seidel

	integer :: k

	do k = 1, 50

		y(1) = ( -F(1)-J(1,2)*y(2) )/J(1,1)
		y(2) = ( -F(2)-J(2,1)*y(1) )/J(2,2)

	end do

end subroutine gauss_seidel

!------------------------------------------------

subroutine newton
	
	integer :: it

	open(1,file = 'saida.txt')

	write(1,1)

	do it = 1, itmax

		call funcoes(x,F)
		call jacobiano
		call gauss_seidel
		x = x+y

		write(1,2) it, x(1), x(2)

	end do

	close(1)

	1 format ('it.', t6, 'X(1)', t31, 'X(2)')
	2 format ( i5, 2(1pe25.15) )

end subroutine newton

!------------------------------------------------

end module processamento