	PROGRAM polint
	USE types
	IMPLICIT NONE

        integer, parameter :: ngx=99  ! size of the grid
        integer, parameter :: n=2     ! order of interpolation
	
        REAL(SP), DIMENSION(ngx) :: xa,ya
	REAL(SP) x
	REAL(SP) :: y,dy
	INTEGER(I4B) :: m,n,ns
	REAL(SP), DIMENSION(n) :: c,d,den,ho

!       initialize coordinates and function values 
        forall(i=1:ngx) xa(i)=i-1
        forall(i=1:ngx:2) ya(i)=3*xa(i)+1 !initalize every second
!                                          the rest is to be inerpolated


!       for each point i=2:ngx-1:2 do the following

	c(1)=xa(i); c(2)=xa(i+2)                    !nodes
	d=c
        x=xa(i+1)                                   !middle point between nodes
	ho(1)=-1;ho(2)=1                            !distance to nodes
	ns=iminloc(abs(x-xa))                       !index of closest node
	y(i+1)=ya(i) or ya(i+2) dep. on ns          !first approx 

	ns=ns-1

	do m=1,n-1                                  !if n=2 no loop!
		den(1:n-m)=ho(1:n-m)-ho(1+m:n)      !den(1)=h(1)-h(2)
		if (any(den(1:n-m) == 0.0)) &
			call nrerror('polint: calculation failure')

		den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
                                                    !den(1)=(c(2)-d(1))/d(1)
		d(1:n-m)=ho(1+m:n)*den(1:n-m)       !d(1) = ho(2)*den(1) 
		c(1:n-m)=ho(1:n-m)*den(1:n-m)       !c(1) = ho(1)*den(1)
		if (2*ns < n-m) then                !if (0<1) then
			dy=c(ns+1)                    dy=c(1) 
		else                                !else 
			dy=d(ns)                      dy=d(0) 
			ns=ns-1                       ns=-1
		end if                              !endif
		y=y+dy                              !y=y+dy
	end do                                      !enddo 
	END program polint
