c===============================================================================
c==== Programa MinimosCuadrados-PolinomiosOrtogonales - EVALUACION DEL POLINOMIO
c===============================================================================
c
c     --------------------------------------------------------------------------
c
c     Evalua los polinomios
c
c        P (x) =  SUM   al  B (x),  con al = <B ,f> / <B ,B >,   j=0,...,m
c         j      i=0,j    i  i            i    i        i  i
c
c     en la base ortogonal:
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     --------------------------------------------------------------------------

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

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

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

      read(1,'(i10)') m
      allocate(dl(1:m),dm(1:m),ds(1:m),
     .         al(0:m),p(0:m),e(0:m),
     .         STAT=ist)
      if (ist.ne.0) call FError(0,'no hay memoria para calculos',2)

        read(1,'(10x,75x,2e25.16)')                  al(0),e(0)
      do i=1,m
        read(1,'(10x,5e25.16)')    dl(i),dm(i),ds(i),al(i),e(i)
      enddo

      close(unit=1)

    1 write(6,'(/a,$)') ' Intervalo de muestreo [xmin,xmax] : '
      read(5,*) xmin,xmax
      if (xmin.ge.xmax) goto 1

    2 write(6,'(/a,$)') ' Numero de subintervalos [np > 0]  : '
      read(5,*) np
      if (np.le.0) goto 2

c.....EVALUACION DEL POLINOMIO
      open(unit=1,file='ajuste.dat')
      write(1,'(15x,11e15.6)') (e(j),j=0,m)

      h=(xmax-xmin)/dble(np)
      do ip=0,np
        x=xmin+dble(ip)*h
        do i=0,m

c.........Evaluacion de polinomios ortogonales
          if (i.eq.0) then
            b(0)=1.d+00
          else if (i.eq.1) then
            b(1)=b(0)
            b(0)=dl(i)*(x-dm(i))*b(1)
          else
            b(2)=b(1)
            b(1)=b(0)
            b(0)=dl(i)*(x-dm(i))*b(1)+ds(i)*b(2)
          endif

c.........Evaluacion del polinomio de minimos cuadrados
          j=i
          if (j.eq.0) then
            p(j)=       al(j)*b(0)
          else
            p(j)=p(j-1)+al(j)*b(0)
          endif

        enddo

        write(1,'(12e15.6)') x,(p(j),j=0,m)
      enddo

      close(unit=1)

      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

