!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Esta subrotina calcula a razao entre o polinomio e a derivada do polinomio de Jacobi
!NAO E RECOMENDADO ALTERACOES NESTA ROTINA POR PARTE DO USUARIO
REAL*8 FUNCTION pjac (n, al, be, x)

	IMPLICIT NONE
		
	!Variaveis definidas na chamada da funcao
	INTEGER :: n		!Numero de pontos de colocacao
	REAL*8 :: al		!Valor de alfa
	REAL*8 :: be		!Valor de beta
	REAL*8 :: x			!Valor de x	

	!Variaveis utilizadas no modulo
	INTEGER :: i			!Contador
	REAL*8 :: p, q, pol		!Polinomio P(i-1)=p ; P(i-2)=q e p(i) = pol
	REAL*8 :: dp, dq, dpol	!Derivada do polinomio dP(i-1)=dp ; dP(i-2)=dq e dp(i) = dpol 
	REAL*8 :: g, h			!Declaracao de g e h

	!Inicializacao das variaveis
	q = 1.0d0
	dq = 0.0d0
	p = (x - ( (1.0d0+be) / (2.0d0+al+be) ))		
	dp = 1.0d0
	pol = p
	dpol = dp
	i=2

	loop1 : DO i =2, n
	
		g = 0.5d0 * (1.0d0 - ( (al**2 - be**2) / ( (2.0d0*i+al+be-1)**2 - 1.0d0) ) )

		h =  ((i-1.0d0)*(i+al-1.0d0)*(i+be-1.0d0)*(i+al+be-1.0d0)) / ((2.0d0*i+al+be-1.0d0)*(2.0d0*i+al+be-2.0d0)**2*(2.0d0*i+al+be-3.0d0))

		!Calculo do polinomio 
		pol = ((x-g)*p)-(h*q)

		!Calculo da derivada do polinomio 
		dpol = p + ((x-g)*dp) - (h*dq)

		!Reatribuicao dos valores para o calculo iterativo
												
		q = p	!P(i-2) = P(i-1)

		dq = dp	!dP(i-2) = dP(i-1)

		p = pol	!P(i-1) = P(i)
		!WRITE (*,*) p
		!dP(i-1) = dP(i)

		dp = dpol
		!WRITE (*,*) dp
							
	ENDDO loop1

	!Retorno da funcao 
	pjac = pol / dpol
				
ENDFUNCTION pjac

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Esta subrotina calcula o valor das raizes do polinomio de jacobi
!NAO E RECOMENDADO ALTERACOES NESTA ROTINA POR PARTE DO USUARIO
SUBROUTINE raizes(n,nt,al,be,Raiz)

	IMPLICIT NONE
		
	!Variaveis declaradas na chamada da rotina
	INTEGER,INTENT(IN) :: n			!Numero de pontos de colocacao interno
	INTEGER,INTENT(IN) :: nt		!Numero de pontos de colocacao totais (internos + externos)
	REAL*8, INTENT(IN) :: al		!Valor de alfa
	REAL*8, INTENT(IN) :: be		!Valor de beta
	REAL*8, INTENT(OUT) :: Raiz(nt)	!Recebe as raizes do polinomio de Jacobi

	!Variaveis utilizadas no modulo
	REAL*8 :: pjac		!Recebe a razao entre o polinomio e a derivada do polinomio de Jacobi
	REAL*8 :: x			!Valor de Xi(i)
	REAL*8 :: y			!Valor de Xi(i+1)
	REAL*8 :: dx		!Diferenca entre Xi(i+1) e Xi(i)
	REAL*8 :: Razao		!Recebe a razão p/p'
	REAL*8 :: soma		!Variavel auxiliar para soma de 1/x-Raiz(i)
	INTEGER :: i, j		!Contadores 

	!Inicializacao das variaveis
	Soma=0.0d0
	dx=1.0d0
	x = 0.0d0

	!Atribuindo o valor para primeira raiz
	Raiz(1)=0.0d0

	!Atribuindo o valor um a n'essima raiz
	Raiz(nt)=1.0d0

	!Calculo da segunda raiz do polinomio
	loop1 : DO 
		
		!x(k+1)= x(k) - f(x)/f'(x)			
		y = x - (pjac(n, al, be, x))
		
		!Modulo de x(k+1)-x(k)
		dx = DABS(y - x)

		!x(k+1)= x(k) para nova iteracao
		x = y
			
		if1 : IF (dx < 10.0d-9) THEN
						
			EXIT

		END IF if1
				
	ENDDO loop1

	!Atribuicao do valor da Raiz ao vetor
	Raiz(2) = x

	!Cálculo das demais raízes do polinomio
	loop2 : DO i=2, n

		!Inicializacao das variaveis
		x = (Raiz(i) + (10.0d-3))
		dx = 1.0d0
		soma= 0.0d0

		loop3 : DO 

			!Reducao do grau do polinomio para n-i
			loop4 : DO j=1, i

				soma = soma + (1/ (x-Raiz(j)))

			ENDDO loop4 
								
			!f(x)/f'(x)				
			Razao = pjac(n, al, be, x)
			
			!Calculo de x(k+1)
			y = x - ( Razao / ( 1.0d0 - (Razao * soma)))
			
			!Modulo de x(k+1)-x(k)
			dx = DABS(y - x)
			
			!x(k+1)= x(k) para nova iteracao
			x = y
			
			!Verifica se a tolerancia foi atingida
			if2 : IF (dx < 10.0d-9) THEN
				
				EXIT

			ENDIF if2
													
		ENDDO loop3
					
		!Alocação dos valores das raízes
		Raiz(i+1) = x
					
	ENDDO loop2 
	
END SUBROUTINE raizes

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Esta subrotina calcula a matriz A(aproximacao para dy/dx) e a Matriz B(aproximacao para d2y/dx2)
!NAO E RECOMENDADO ALTERACOES NESTA ROTINA POR PARTE DO USUARIO
SUBROUTINE Matriz_AB(n,nt,Raiz,A,B)

	IMPLICIT NONE

	!Variaveis declaradas na chamada da rotina
	INTEGER,INTENT(IN) :: n				!Numero de pontos de colocacao interno
	INTEGER,INTENT(IN) :: nt			!Numero de pontos de colocacao totais (internos + externos)
	REAL*8,INTENT(IN) :: Raiz(nt)		!Recebe as raizes do polinomio de Jacobi
	REAL*8,INTENT(OUT) :: A(nt,nt)		!Recebe a matriz A - Aproximacao para primeira derivada
	REAL*8,INTENT(OUT) :: B(nt,nt)		!Recebe a matriz B - Aproximacao para segunda derivada

	!Variaveis utilizadas no modulo
	REAL*8 :: dxi			!Diferenca entre (x(i)-x(j))
	REAL*8 :: pi			!Polinomio nodal. Procedimento recursivo
	REAL*8 :: dpi			!Primeira derivada do polinomio nodal. Procedimento recursivo
	REAL*8 :: d2pi			!Segunda derivada do polinomio nodal. Procedimento recursivo
	REAL*8 :: d3pi			!Terceira derivada do polinomio nodal. Procedimento recursivo
	REAL*8 :: dpnodal(nt)	!Valor da primeira derivada do polinomio nodal em x(i)
	REAL*8 :: d2pnodal(nt)	!Valor da segunda derivada do polinomio nodal em x(i)
	REAL*8 :: d3pnodal(nt)	!Valor da terceira derivada do polinomio nodal em x(i)
	INTEGER :: i,j			!Contadores

	!Inicio do procedimento
	loop1 : DO i=1, nt
			
		!Inicializacao das variaveis
		pi=1.0d0 
		dpi=0.0d0                          
		d2pi=0.0d0                         
		d3pi=0.0d0                         
			
		!Calculo iterativo para determinacao p,dp, d2p e d3p no ponto i
		loop2 : DO j=1, nt			   
		
			!Primeira derivada do polinomio nodal	
			dxi=(Raiz(i)-Raiz(j))
			!Terceira derivada do polinomio nodal
			d3pi = (dxi*d3pi) + (3.0d0*d2pi)
			!Segunda derivada do polinomio nodal
			d2pi = (dxi*d2pi) + (2.0d0*dpi)
			!Primeira derivada do polinomio nodal
			dpi = (dxi*dpi) + pi
			!Polinomio nodal
			pi = (dxi*pi)
				
		ENDDO loop2

		!Atribuicao do valor da primeira derivada do polinomio nodal em x(i)
		dpnodal(i)=dpi						
		!Atribuicao do valor da segunda derivada do polinomio nodal em x(i)
		d2pnodal(i)=d2pi                    
		!Atribuicao do valor da terceira derivada do polinomio nodal em x(i)	
		d3pnodal(i)=d3pi                    

	ENDDO loop1

	!Procedimento para determinacao da matriz A e B
	loop3 : DO i=1, nt						

		!Valor de A quando i=j
		A(i,i) = (d2pnodal(i)/(2.0d0*dpnodal(i)))			
		!Valor de B quando i=j
		B(i,i) = (d3pnodal(i)/(3.0d0*dpnodal(i)))			
	
		loop4 : DO j=1, nt
				
			!Valor de A e B para i diferente de j
			if1 : IF (i.NE.j) THEN									
				
				!Valor de A quando i/=j
				A(i,j) = dpnodal(i)/ ((Raiz(i)-Raiz(j))*dpnodal(j))	
				!Valor de B quando i/=j
				B(i,j) = 2.0d0*A(i,j)*(A(i,i)-(1.0d0/(Raiz(i)-Raiz(j))))		

			ENDIF if1

		ENDDO loop4

	ENDDO loop3

ENDSUBROUTINE Matriz_AB

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Esta subrotina calcula os pesos das quadraturas de Gauss-Lobatto
!NAO E RECOMENDADO ALTERACOES NESTA ROTINA POR PARTE DO USUARIO
SUBROUTINE Pesos(nt,Raiz,W)
			
	IMPLICIT NONE

	!Variaveis declaradas na chamada da rotina
	INTEGER,INTENT(IN) :: nt		!Numero de pontos de colocacao totais (internos + externos)
	REAL*8,INTENT(IN) :: Raiz(nt)	!Recebe as raizes do polinomio de Jacobi
	REAL*8,INTENT(OUT) :: W(nt)		!Recebe os pesos da quadratura de Lobatto 

	!Variaveis utilizadas no modulo
	REAL*8 :: dxi		!Diferença entre (x(i)-x(j))
	REAL*8 :: pi		!Polinomio nodal. Procedimento recursivo
	REAL*8 :: dpi		!Primeira derivada do polinomio nodal. Procedimento recursivo
	REAL*8 :: SomaP		!Somatorio de todos os pesos para normalizacao
	INTEGER :: i,j		!Contadores

	!Inicializacao das variaveis
	SomaP=0.0d0	

	!Inicio do procedimento
	loop1 : DO i=1, nt
	
		!Inicializacao das variaveis
		pi=1.0d0 
		dpi=0.0d0                          

		!Calculo iterativo para determinacao p,dp no ponto i
		loop2 : DO j=1, nt			   
				
			!Primeira derivada do polinomio nodal	
			dxi=(Raiz(i)-Raiz(j))
			!Terceira derivada do polinomio nodal
			dpi = (dxi*dpi) + pi
			!Polinomio nodal
			pi = (dxi*pi)
				
		ENDDO loop2

		!Atribuicao do valor do peso(i)
		W(i) = 1.0d0 / (dpi**2.0d0)
			
		!Somatorio dos valores dos pesos
		SomaP=W(i)+SomaP
		!WRITE(*,*) W(i)
		!WRITE(*,*) SomaP
		!WRITE(*,*) ''
		!PAUSE
		
	ENDDO loop1

	!Normalizacao dos valores dos pesos
	loop3 : DO i=1, nt

		W(i)=W(i)/SomaP

	ENDDO loop3

ENDSUBROUTINE Pesos

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Esta subrotina calcula a matriz V para aplicacao do metodo dos momentos
!NAO E RECOMENDADO ALTERACOES NESTA ROTINA POR PARTE DO USUARIO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Esta subrotina calcula o polinomio interpolador de Lagrange
!NAO E RECOMENDADO ALTERACOES NESTA ROTINA POR PARTE DO USUARIO

SUBROUTINE lagrange(nt, Raiz, z, lag);

	IMPLICIT NONE

	!Variaveis declaradas na chamada da rotina
	INTEGER,INTENT(IN) :: nt	    	!Numero de pontos de colocacao totais (internos + externos)
	REAL*8,INTENT(IN)  :: Raiz(nt)   	!Recebe as raizes do polinomio de Jacobi
	REAL*8,INTENT(IN)  :: z		     	!Ponto para o qual se deseja o valor de Y interpolado
	REAL*8,INTENT(OUT) :: lag(nt)      !y interpolado

	!Variaveis auxiliares do modulo
	INTEGER :: j,k	!Contadores
	REAL*8 :: aux !auxiliar
	!Inicializacao das variaveis
	

	!Inicio do calculo dos polinomios interpoladores de Lagrange	
	loop1 : DO j=1, nt
			 
		!Variavel auxiliar utilizada para fazer o produtorio
		
		aux = 1.0d0				  
		
		loop2 : DO k=1, nt

		
	
			!Verifica se i e diferente de j
			if1 : IF (j.NE.k) THEN
			
				!Produtorio(j=1 k./=.j,n) de (x-xk)/(xj-xk)
				aux = ((z-Raiz(k))/(Raiz(j)-Raiz(k)))*aux
	
			ENDIF if1
			  
			lag(j) = aux
		ENDDO loop2
		

	ENDDO loop1


ENDSUBROUTINE lagrange
