module mod

    implicit none
    private                  ! tutte le quantita' non esplicitamente public
                             ! sono private (visibili solo nel modulo)
    integer,parameter        :: k=selected_real_kind(14)
    public                   :: Poly, f, set_interval, k

    real(kind=k)                           :: a,b
    real(kind =k),dimension(:),allocatable :: nodo

  contains

! set_interval copia gli argomenti x e y nelle variabili (private) globali
! a e b del modulo: ogni sottoprogramma del modulo ha accesso ad a e b ma 
! non i programmi esterni al modulo
  subroutine set_interval(x,y)
    real(kind=k),intent(in) :: x,y
    a = x
    b = y
  end subroutine set_interval

! f(x) : funzione da interpolare 
  function f(x) result(y)
    real(kind =k),intent(in)  ::x
    real(kind =k)             ::y
    if ( x == 0) then 
         y = 0         ! prolungamento per continuita'  nell' origine
    else    
         y= exp (-1.0/x**2)
    end if
  end function f

  ! funzione  di Lagrange relativa al j-esimo nodo dell' interpolante
  ! di grado n
  function Lagrange(n,j,x)  result (u)
    real(kind =k),intent(in)               :: x  
    integer,intent(in)                     :: n,j
    real(kind =k)                          :: den,u,num
    integer                                :: i
    
    num=1.0
    den=1.0
    do i=0,n
    if (i == j) cycle                  !il prodotto deve escludere x==xi 
          num=num*(x-nodo(i))          !prodotto dei numeratori
          den=den*(nodo(j)-nodo(i))    !prodotto dei denominatori
    end do
    u=num/den
  end function Lagrange
  
! interpolante polinomiale di grado g nel punto x
! la function calcola i nodi
  function Poly(g,x) result(pol)
    integer, intent(in)          :: g      ! grado del polinomio (n.punti - 1)
    real(kind =k),intent(in)     :: x      ! 
    integer                      :: i,errall
    real(kind =k)                :: pol
  
    allocate(nodo(0:g),stat=errall)
    if(errall /= 0) stop 'errore alloc. nodo '
    do i=0,g
      nodo(i) = a + i*(b-a)/g
    end do
    pol=0.0
    do  i = 0, g
       pol = pol + f(nodo(i))*Lagrange(g,i,x) !formula di Lagrange
    end do
    deallocate(nodo)
  end function Poly
  
end module mod


program lagra
use mod
implicit none
integer        :: m, i, j, mi, ma, step
real(kind =k)  :: x, dx, a, b

a =-1.0
b = 1.0
call set_interval(a,b)  ! passa al modulo gli estremi dell' intervallo
do
   print*,"numero di punti (equispaziati) sui cui valutare funzione ed interpolanti per graficarle"
   print*,"(deve essere > 1)"
   read*,m
   if(m > 1)exit
end do
print*," grado min, max e step delle interpolanti"
print*," (verranno calcolate le int. di grado min,min+step, min+2*step..."
print*," per tutti i valori di n che garantiscono  min + n*step <= max )"
do
read*,mi,ma,step
if( mi>=1   .and. &
    ma>=mi  .and. &
    step>=1        )exit
print*," deve essere : min >=1, max>=min, step >=1 "
end do

dx = (b-a)/(m-1)
do i = 1,m
   x = a + (i-1)*dx
   write(unit = 1,fmt = *)x,f(x)
   write(unit = 2,fmt = *)x,(Poly(j,x),j=mi,ma,step)
end do
print*," il file fort.1 contiene i valori dell' ascissa e della funzione "
print*," il file fort.2 contiene i valori dell' ascissa e delle interpolanti su ",&
        2+(ma-mi)/step," colonne"
print*,"  "
print*," il file fort.3 contiene comandi per visualizzare le curve mediante gnuplot"
print*," uso:  gnuplot fort.3 "
print*," e poi invio per terminare "
write(unit=3,fmt=*)"set style data lin "
write(unit=3,fmt=*)"p 'fort.1' ,'fort.2' \"
i=3
do j=mi+1,ma,step
   write(unit=3,fmt=*)",'' u 1:",i,"\"
   i=i+1
end do
write(unit=3,fmt=*)" "
write(unit=3,fmt=*)"pause -1"
end program lagra