c===============================================================================
c====================== Programa MinimosCuadrados-PolinomiosOrtogonales - AJUSTE
c===============================================================================
c
c     --------------------------------------------------------------------------
c
c     Dada la nube de puntos {(x ,y )}            con los pesos {w }
c                               k  k  k=0,...,n                   k k=0,...,n
c
c     donde y  = f(x ),
c            k      k
c
c     y el producto escalar <g(x),h(x)> =  SUM   w  g(x ) h(x ),
c                                         k=0,n   k    k     k
c
c     calcula los polinomios P (x) que minimizan los Q = <e (x),e (x)>,
c                             j                       j    j     j
c
c     donde  e (x) = f(x) - P (x), con j=0,...,m.
c             j              j
c
c     Para ello se calcula la base de polinomios ortogonales {B }
c                                                              i i=0,...,m
c     tal que:
c
c        B (x) = 1; B (x) = dl (x-dm ) B (x);
c         0          1        1     1   0
c        B (x) = dl (x-dm ) B   (x) + ds  B   (x),  i=2,...,m
c         i        i     i   i-1        i  i-2
c
c     donde:
c
c         dl  = 1,         dm    =   <xB   ,B   > / <B   ,B   >,   
c           i                i          i-1  i-1      i-1  i-1
c
c         ds = -(dl / dl   ) < B   ,B   > / <B   ,B   >.   
c                  i    i-1     i-1  i-1      i-2  i-2  
c     
c     Finalmente, el polinomio de mínimos cuadrados se escribe en la forma:
c
c        P (x) =  SUM   al  B (x),  con al = <B ,f> / <B ,B >.
c         j      i=0,j    i  i            i    i        i  i
c
c     Se calculan también los valores E = sqrt(Q  / <1,1>),   j=0,...,m.
c                                      j        j
c
c     --------------------------------------------------------------------------

      implicit real*8 (a-h,o-z)

      allocatable x(:),y(:),w(:),p(:),b(:,:)   !b(di,k)->B_{i-di}(x_k), di=0,1,2
      allocatable dl(:),dm(:),ds(:),
     .            bb(:),xbb(:),bf(:),al(:),e(:)  !e(j)=sqrt(<f-P_j,f-P_j>/<1,1>)          

      character formato*40

c.....LECTURA DE DATOS
      open(unit=1,file='datos.dat',IOSTAT=ist)
      if (ist.ne.0) call FError(0,'imposible abrir datos.dat',1)

      read(1,'(i10,1x,a)') n,formato
      allocate(x(0:n),y(0:n),w(0:n),p(0:n),b(0:2,0:n),STAT=ist)
      if (ist.ne.0) call FError(0,'no hay memoria para datos',2)

      do i=0,n
        read(1,formato) x(i),y(i),w(i)
      enddo
      close(unit=1)

C.....GRADO MAXIMO
    1 write(6,100)
  100 format(/' Introduzca GRADO del polinomio [m>0] = ',$)
      read(5,*) m
      if (m.lt.1) goto 1
      allocate(dl(1:m),dm(1:m),ds(1:m),
     .         bb(0:m),xbb(0:m),bf(0:m),al(0:m),e(0:m),
     .         STAT=ist)
      if (ist.ne.0) call FError(0,'no hay memoria para calculos',3)


c.....APROXIMACION POR MINIMOS CUADRADOS
      do i=0,m

c.......Obtencion de polinomios ortogonales
        if (i.eq.1) then
          dl(i)=1.d+00
          dm(i)=xbb(i-1)/bb(i-1)
          ds(i)=0.d+00
        else if (i.gt.1) then
          dl(i)=1.d+00
          dm(i)=xbb(i-1)/bb(i-1)
          ds(i)=-(dl(i)/dl(i-1))*(bb(i-1)/bb(i-2))
        endif

        bb (i)=0.d+00
        xbb(i)=0.d+00
        bf (i)=0.d+00
        do k=0,n
          if (i.eq.0) then
            b(0,k)=1.d+00
          else if (i.eq.1) then
            b(1,k)=b(0,k)
            b(0,k)=dl(i)*(x(k)-dm(i))*b(1,k)
          else
            b(2,k)=b(1,k)
            b(1,k)=b(0,k)
            b(0,k)=dl(i)*(x(k)-dm(i))*b(1,k)+ds(i)*b(2,k)
          endif
          bb (i)=bb (i)+w(k)*b(0,k)*b(0,k)
          xbb(i)=xbb(i)+w(k)*x(k)*b(0,k)*b(0,k) !aunque xbb(m) no se utiliza
          bf (i)=bf (i)+w(k)*b(0,k)*y(k)
        enddo

c.......Ajuste por minimos cuadrados
        j=i
        al(j)=bf(j)/bb(j) !si los datos son incorrectos bb(j) podría ser 0
        e(j)=0.d+00
        do k=0,n
          if (j.eq.0) then
            p(k)=     al(j)*b(0,k)
          else
            p(k)=p(k)+al(j)*b(0,k)
          endif
          e(j)=e(j)+w(k)*(y(k)-p(k))**2
        enddo
        e(j)=e(j)/bb(0)
        e(j)=sqrt(e(j))
      enddo

c.....ESCRITURA DE RESULTADOS
      open(unit=1,file='polcoef.dat')
      write(1,'(i10)') m
        write(1,'(i10,75x,2e25.16)') 0,              al(0),e(0)
      do i=1,m
        write(1,'(i10,5e25.16)') i,dl(i),dm(i),ds(i),al(i),e(i)
      enddo
      close(unit=1)

      call exit(0)

      end


c_________________________________________________________________________FError
      Subroutine FError(logerr,error,ierr)
      character*(*) error

      write(logerr,100) error
  100 format(/,' MCPO:: Error: ',a)

      write(logerr,'(/a,$)')
     .      ' Fin de Ejecucion. Oprima RETURN para salir.'
      read(5,'()')
      call exit(ierr)

      end



