subroutine INTERPOL_2D_TCHE(MAILLX,MAILLY,CHAMP,POSIT,INTERPOL)
!=========INTERPOL_2D_TCHE=================
!Realise l'interpolation 2D a base de polynomes orthogonaux de Tchebychev.
!Entrees : MAILLX, MAILLY : les points de collocations, VALEURS : du champ
!aux points de collocation, POSIT : l'abscice et ordonnee du point ou l'on souhaite realiser
!l'interpolation.
!Sortie : INTERPOL : la valeur interpolee de la fonction au point POSIT
!La resolution du systeme lineaire afin d'obtenir les different coefficients a_k se fait
!via une methode LU montante/descendante.
!MARLAN 2011
implicit none
	real(kind=8), dimension(:),intent(in) :: MAILLX,MAILLY, POSIT
	real(kind=8), dimension(:,:), intent(in) :: CHAMP
	real(kind=8), intent(out) :: INTERPOL
	
	interface
		subroutine INTERPOL_TCHE(MAILL,VALEURS,POSIT,INTERPOL)
		implicit none
		real(kind=8),dimension(:),intent(in) :: MAILL
		real(kind=8),dimension(:),intent(in) :: VALEURS
		real(kind=8),intent(in) :: POSIT
		real(kind=8),intent(out) :: INTERPOL
		end subroutine INTERPOL_TCHE
	end interface
	
	real(kind=8), dimension(size(MAILLX)) :: INTERPOLX
	integer :: i

	!Creation du vecteur d'interpolation sur toutes les colonnes du champ entrant
	do I=1,size(MAILLY)
		call INTERPOL_TCHE(MAILLX,CHAMP(1:size(MAILLX),I),POSIT(1),INTERPOLX(I))
	end do
	!~ print*,size(INTERPOLX),size(MAILLX),size(MAILLY)
	call INTERPOL_TCHE(MAILLY,INTERPOLX,POSIT(2),INTERPOL)

end subroutine INTERPOL_2D_TCHE



subroutine INTERPOL_TCHE(MAILL,VALEURS,POSIT,INTERPOL)
implicit none
!=========INTERPOL_TCHE=================
!Realise l'interpolation 1D a base de polynomes orthogonaux de Tchebychev.
!Entrees : MAILL : les points de collocations, VALEURS : les valeurs de la fonction
!aux points de collocation, POSIT : l'abscice/ordonnee du point ou l'on souhaite realiser
!l'interpolation.
!Sortie : INTERPOL : la valeur interpolee de la fonction au point POSIT
!La resolution du systeme lineaire afin d'obtenir les different coefficients a_k se fait
!via une methode LU montante/descendante.
!MARLAN 2011
	real(kind=8),dimension(:),intent(in) :: MAILL
	real(kind=8),dimension(:),intent(in) :: VALEURS
	real(kind=8),intent(in) :: POSIT
	real(kind=8),intent(out) :: INTERPOL
	
	interface
		subroutine MAILL_MOD_TCHE(MAILL,MAILL_MOD,X,X_MOD)
		implicit none
		real(kind=8),dimension(:),intent(in) :: MAILL
		real(kind=8),dimension(:),intent(OUT) :: MAILL_MOD
		real(kind=8), intent(in) :: X
		real(kind=8), intent(out) :: X_MOD
		end subroutine MAILL_MOD_TCHE
		
		subroutine TCHE_POL(X,NB_COEF,T_N)
		implicit none
		integer,intent(in) :: NB_COEF
		real(kind=8),intent(in) :: X
		real(kind=8),dimension(NB_COEF),intent(out) :: T_N
		end subroutine TCHE_POL
		
		subroutine RESOL_LU(T_NGTAB,a_k,SOL)
		implicit none
		real(kind=8),dimension(:,:),intent(in) :: T_NGTAB
		real(kind=8),dimension(:),intent(in) :: SOL
		real(kind=8),dimension(:),intent(out) :: a_k
		end subroutine RESOL_LU
	end interface
	
	real(kind=8),dimension(size(MAILL)):: MAILL_MOD,T_k,a_k
	real(kind=8),dimension(size(MAILL),size(MAILL)):: T_N
	real(kind=8) :: POSIT_MOD
	integer :: I,J,K
	
	call MAILL_MOD_TCHE(MAILL,MAILL_MOD,POSIT,POSIT_MOD)

	do I=1,size(MAILL)
		call TCHE_POL(MAILL_MOD(I),size(MAILL),T_N(:,i))
	end do

	call TCHE_POL(POSIT_MOD,size(MAILL),T_k)

	call RESOL_LU(T_N(:,:),a_k,VALEURS)
	
	INTERPOL=0.0
	do I=1,size(MAILL)
		INTERPOL=INTERPOL+a_k(I)*T_k(I)
	end do
end subroutine INTERPOL_TCHE

subroutine MAILL_MOD_TCHE(MAILL,MAILL_MOD,X,X_MOD)
implicit none
	real(kind=8),dimension(:),intent(in) :: MAILL
	real(kind=8),dimension(:),intent(OUT) :: MAILL_MOD
	real(kind=8), intent(in) :: X
	real(kind=8), intent(out) :: X_MOD
	
	MAILL_MOD=(2*MAILL-maxval(MAILL,dim=1)-minval(MAILL,dim=1))/(maxval(MAILL,dim=1)-minval(MAILL,dim=1))
	X_MOD=(2*X-maxval(MAILL,dim=1)-minval(MAILL,dim=1))/(maxval(MAILL,dim=1)-minval(MAILL,dim=1))
end subroutine MAILL_MOD_TCHE

subroutine TCHE_POL(X,NB_COEF,T_N)
implicit none
	integer,intent(in) :: NB_COEF
	real(kind=8),intent(in) :: X
	real(kind=8),dimension(NB_COEF),intent(out) :: T_N
	real(kind=8) :: T_NM1,T_NM2,T_0,T_1
	integer :: i,j,k,l
	
	T_N(1)=1.0
	T_N(2)=X
	if(NB_COEF>2) then
		do k=3,NB_COEF
			T_N(k)=2*X*T_N(k-1)-T_N(k-2)
		end do
	end if
end subroutine TCHE_POL

subroutine RESOL_LU(T_NGTAB,a_k,SOL)
implicit none
	real(kind=8),dimension(:,:),intent(in) :: T_NGTAB
	real(kind=8),dimension(:),intent(in) :: SOL
	real(kind=8),dimension(:),intent(out) :: a_k
	real(kind=8) :: TEMPL, TEMPU, TEMP
	real(kind=8),dimension(size(T_NGTAB,dim=1),size(T_NGTAB,dim=2)) :: L, U
	real(kind=8),dimension(size(SOL)) :: Y
	integer :: i,j,k,n
	
	n=size(T_NGTAB,dim=1)
	U=0.0
	L=0.0
	
	forall(i=1:n,j=1:n,i==j)
		U(i,j)=1
	End forall
	
	DO j=1,n
		DO i=1,n
			TEMPU=0.0
			TEMPL=0.0
			IF(i>1) THEN
				DO k=1,i-1
					TEMPU=L(k,i)*U(j,k)+TEMPU
				END DO
			ELSE 
				TEMPU=0.0
			END IF
						
			IF(j>1) THEN
				
				DO k=1,j-1
					TEMPL=L(k,i)*U(j,k)+TEMPL
				END DO
				
			ELSE 
				TEMPL=0.0
			END IF
			
			L(j,i)=(T_NGTAB(j,i)-TEMPL)/U(j,j)
			U(j,i)=(T_NGTAB(j,i)-TEMPU)/L(i,i)
			U(i,i)=1.0
			END DO
	END DO
	
	FORALL(i=1:n,j=1:n,i<j)
		L(j,i)=0.0
		U(i,j)=0.0
	END FORALL

	Y(1)=SOL(1)/L(1,1)
	
	Do i=2,n
		TEMP=0.0
		
		do j=1,i-1
			TEMP=L(j,i)*Y(j)+TEMP
		end do
		
		Y(i)=(SOL(i)-TEMP)/L(i,i)
	end do

	a_k(n)=Y(n)/U(n,n)

	DO i=n-1,1,-1
		TEMP=0.0
		
		DO j=i+1,n
			TEMP=U(j,i)*a_k(j)+TEMP
		END DO
		
		a_k(i)=(Y(i)-TEMP)/U(i,i)
	END DO
end subroutine RESOL_LU