! Polinomul de interpolare Newton
! Calculeaza valorle in "nz" puncte "z"
!
! Fisier de intrare / iesire
! Date de intrare:
! 1) Titlu: <=80 caractere
! 2) kod_functie: #0: definita analitic cu rutina "fun.f90"; =0: prin valori citite.
! 3) n1: numar de noduri
! 4) x(1:n1): valorile "x"
! [NUMAI pt. kod_functie =0:
! 4-bis) f(1:n1): valorile "f(x)".] 
! 5) nz: numarul valorilor z ("x"-uri pt. valori interpolate)
! 6) z(1:nz): valori "z"
!



character(255) fname, dummy*80
real(4), allocatable :: x(:), f(:), z(:), c(:,:)

print '(a,\)', 'Fisier de date: '
read '(a)', fname

open(1, file =fname , status ='old')

read(1,'(a)') dummy

read(1,*) kod_functie	! 

read(1,*) n1	! nr. de noduri: n1 =n+1

n = n1 -1
Allocate (x(0:n), f(0:n), c(0:n, 0:n))
x =0; f =0; c =0.;

read(1,*) x

if(kod_functie ==0) then
read(1,*) f
endif

read(1,*) nz
allocate(z(nz))
z =0.

read(1,*) z

if(kod_functie /=0) then
do i =0, n
f(i) = fun(x(i))
enddo
endif

c =0.
c(:,0) =f

do k = 1, n		! coloana
do j = 0, n-k	! linia
c(j,k) = (c(j+1, k-1) - c(j, k-1))/(x(j+k) - x(j)) ! o.k.
enddo
enddo

! Rezultate
write(1, '(/, t4,a, t18,a)') 'x', 'f(x)'
do i=0,n
write(1,*) x(i), f(i)
enddo



write(1, '(/, a)') ' Coeficienti c_ij:'
do i =0, n
write(1,*) (c(i,j), j= 0, n-i)
enddo

! Calculeaza val. polinomului si eroraea:
if(kod_functie ==0) then 
write(1,'(/, t4,a, t18,a)') 'z', 'p(z)'
else
write(1,'(/, t4,a, t18,a, t39, a)') 'z', 'p(z)', 'Eroare relativa'
endif

do iz =1, nz
zz =z(iz)
q =1.; p =c(0,0)
	do i = 1, n
	q =q*(zz - x(i-1))
	p = p + c(0,i)*q
	enddo

IF(kod_functie /=0) THEN 
fz =fun(z(iz)); err =fz -p;

if(abs(fz) >0.2) then
write(1, '(t4,g13.6, t18,1pg14.7, t38,1pg10.3)') z(iz), p, err/fz
else
write(1, '(t4,g13.6, t18,1pg14.6, t38,1pg10.3, a)') z(iz), err, p, '  (* eroarea)'
endif

ELSE

write(1, '(t2,g13.6, t16,1pg14.7)') z(iz), p
ENDIF

enddo

deallocate(x, f, z, c)
END
