!
!  programa de ajuste por minimos quadrados
!
use msflib ! biblioteca que contem o comando systemqq
logical chamada
!
      parameter(nmax=25)   
      dimension x(0:nmax),fx(0:nmax),alfa(0:nmax,0:nmax),b(0:nmax)
      dimension a(nmax+1,nmax+1),aux(nmax+1) 
!
!     arquivo de entrada com Nr pontos dados "n" e os valores "x_i" e "fx(x_i)"
!
      open (9,file='inptab.txt')
!
!     arquivo de saída com a curva ajustada
!
      open (3,file='out-ajust.txt') 
!
!****************************************
!	Leitura do grau do polinomio de ajuste
!	 
!****************************************
!
      write(*,*)'Entre com o grau do polinomio de ajuste:'
      read(*,*)m
!
      write(*,*)'Entre com o valor de n para (n + 1) pontos dados:'
      read(*,*)n
!******************************
!
!
!     ler Nr de pontos dados "n" e os pontos dados a interpolar
!
!	read(9,*)n
      do i=0,n
         read(9,*) x(i),fx(i)
      enddo
!
!	Nr de pontos desejado para a curva ajustada
!
	write(*,*)'Entre com o Nr de pontos desejado para a curva ajustada:'
	read(*,*)naj
!
!	calcular os coeficientes "alfa"
!
	ialfa=0
	do i=0,m
		do j=ialfa,m
			sum=0.
			do k=0,n
				gl=g(i,x(k))
				gr=g(j,x(k))
				sum=sum+gl*gr
			enddo
			if(i.ne.j) then
				alfa(i,j)=sum
				alfa(j,i)=sum
			else
				alfa(i,j)=sum
			endif
		enddo
		ialfa=ialfa+1
	enddo
!
	do i=0,m
		sum=0.
		do k=0,n
			gr=g(i,x(k))
			sum=sum+fx(k)*gr
		enddo
		b(i)=sum
	enddo
!
	do i=0,m
		write(*,*)'b(',i,')=',b(i)
		aux(i+1)=b(i)
		write(*,*)'aux(',i+1,')=',aux(i+1)
		do j=0,m
			write(*,*)'alfa(',i,',',j,')= ',alfa(i,j)
			a(j+1,i+1)=alfa(i,j)
			write(*,*)'a(',i+1,',',j+1,')= ',a(j+1,i+1)
		enddo
	enddo
!
!	obtain the polynomial coefficients
!
	call ludcmp(a,m+1,nmax+1,indx,d)
	call lubksb(a,m+1,nmax+1,indx,aux)
!
	do i=0,m
		write(*,*)'aux(',i+1,')=',aux(i+1)
	enddo
!
	dx=(x(n)-x(0))/naj
	xp=x(0)
	do i=1,naj+1
		fp=0.
		do j=0,m
		fp=fp+aux(j+1)*g(j,xp)
		enddo
		write(3,*)xp,fp
		xp=xp+dx
	enddo

close(3)

chamada = systemqq('notepad out-ajust.txt')       ! listagem dos dados

 chamada = systemqq('wgnuplot dados.gnu')  ! gráfico
      stop
      end

!*************************************************
	function g(i,x)

		aux=x**i
		g=aux

	return
	end



