module spline_rec
  implicit none
  
contains
  
  ! NAASA  1.1.007 DCOPY    FTN-A 05-02-78     THE UNIV OF MICH COMP CTR
      SUBROUTINE  DCOPY(N,DX,INCX,DY,INCY)
!
!     COPIES A VECTOR, X, TO A VECTOR, Y.
!     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
!     JACK DONGARRA, LINPACK, 6/17/77.
!
      real DX(1:N),DY(1:N)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
!
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20
!
!        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
!          NOT EQUAL TO 1
!
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
!
!        CODE FOR BOTH INCREMENTS EQUAL TO 1
!
!
!        CLEAN-UP LOOP
!
   20 M = MOD(N,7)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        DY(I) = DX(I)
        DY(I + 1) = DX(I + 1)
        DY(I + 2) = DX(I + 2)
        DY(I + 3) = DX(I + 3)
        DY(I + 4) = DX(I + 4)
        DY(I + 5) = DX(I + 5)
        DY(I + 6) = DX(I + 6)
   50 CONTINUE
      RETURN
   END 

   ! NAASA  2.1.038 DGTSL    FTN-A 05-02-78     THE UNIV OF MICH COMP CTR
      SUBROUTINE DGTSL(N,L,D,U,B,INFO)
      INTEGER N,INFO
      real L(1:N),D(1:N),U(1:N),B(1:N)
!
!     DGTSL GIVEN A GENERAL TRIDIAGONAL MATRIX AND A RIGHT HAND
!     SIDE WILL FIND THE SOLUTION.
!
!     ON ENTRY
!
!        N       INTEGER
!                IS THE ORDER OF THE TRIDIAGONAL MATRIX.
!
!        L       DOUBLE PRECISION(N)
!                IS THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX.
!                L(2) THROUGH L(N) SHOULD CONTAIN THE SUBDIAGONAL.
!                ON OUTPUT L IS DESTROYED.
!
!        D       DOUBLE PRECISION(N)
!                IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX.
!                ON OUTPUT D IS DESTROYED.
!
!        U       DOUBLE PRECISION(N)
!                IS THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX.
!                U(1) THROUGH U(N-1) SHOULD CONTAIN THE SUPERDIAGONAL.
!                ON OUTPUT U IS DESTROYED.
!
!        B       DOUBLE PRECISION(N)
!                IS THE RIGHT HAND SIDE VECTOR.
!
!     ON RETURN
!
!        B       IS THE SOLUTION VECTOR.
!
!        INFO    INTEGER
!                = 0 NORMAL VALUE.
!                = K IF THE K-TH ELEMENT OF THE DIAGONAL BECOMES
!                    EXACTLY ZERO.  THE SUBROUTINE RETURNS WHEN
!                    THIS IS DETECTED.
!
!     LINPACK. THIS VERSION DATED 07/14/77 .
!     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
!
!     NO EXTERNALS
!     FORTRAN DABS
!
!     INTERNAL VARIABLES
!
      INTEGER K,KB,KP1,NM1,NM2
      DOUBLE PRECISION T

         INFO = 0
!
!        CHECK FOR 1 X 1 CASE
!
         IF (N .NE. 1) GO TO 20
         IF (D(1) .EQ. 0.0D0) GO TO 10
            L(1) = D(1)
            B(1) = B(1)/D(1)
         GO TO 80
   10    CONTINUE
            INFO = 1
            L(1) = D(1)
!     ......EXIT
            GO TO 90
   20       CONTINUE
            L(1) = D(1)
            D(1) = U(1)
            U(1) = 0.0D0
            U(N) = 0.0D0
            NM1 = N - 1
!
            DO 50 K = 1, NM1
               KP1 = K + 1
!
!              FIND THE LARGEST OF THE TWO ROWS
!
               IF (ABS(L(KP1)) .LT. ABS(L(K))) GO TO 30
!
!                 INTERCHANGE ROW
!
                  T = L(KP1)
                  L(KP1) = L(K)
                  L(K) = T
                  T = D(KP1)
                  D(KP1) = D(K)
                  D(K) = T
                  T = U(KP1)
                  U(KP1) = U(K)
                  U(K) = T
                  T = B(KP1)
                  B(KP1) = B(K)
                  B(K) = T
   30          CONTINUE
!
!              ZERO ELEMENTS
!
               IF (L(K) .NE. 0.0D0) GO TO 40
                  INFO = K
!     ............EXIT
                  GO TO 90
   40          CONTINUE
               T = (-1.0D0/L(K))*L(KP1)
               L(KP1) = D(KP1) + T*D(K)
               D(KP1) = U(KP1) + T*U(K)
               U(KP1) = 0.0D0
               B(KP1) = B(KP1) + T*B(K)
   50       CONTINUE
!
!           BACK SOLVE
!
            NM2 = N - 2
            B(N) = B(N)/L(N)
            B(NM1) = (B(NM1) - D(NM1)*B(N))/L(NM1)
            IF (N .EQ. 2) GO TO 70
               DO 60 KB = 1, NM2
                  K = NM2 - KB + 1
                  B(K) = (B(K) - D(K)*B(K+1) - U(K)*B(K+2))/L(K)
   60          CONTINUE
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
!
      RETURN
   END 
   
!       SUBROUTINE GETDATA (X,Y,N,MAXN)

!       DOUBLE PRECISION X(MAXN), Y(MAXN)

!       CHARACTER*128 FNAME

!       PRINT *,'Name of input file:'
!       READ (*,'(A128)') FNAME
!       OPEN (UNIT=20, FILE=FNAME, STATUS='OLD')
!       I=0
! 10    I=I+1
!       READ (20,*,END=30) X(I),Y(I)
!       IF (I.LT.MAXN) GOTO 10
!       STOP 'I .LT. MAXN in GETDATA'

! 30    CLOSE (20)
!       N=I-1

!       RETURN
!       END
   
   !> Na vstupu predpoklada soubor bodu (X(I), Y(I)) (kazda dvojice na samostatnem
   !> radku).
   !> Body jsou zadavany postupne v tom poradi, v jakem se budou prokladat.
   !>

   subroutine spline_reconstruct(nn, x, y, TK, XD, YD, layer) 
     integer, intent(in) :: nn   ! number of nodes is nn
     real, dimension(0:nn+1), intent(inout) :: x, y  ! values to be interpolated
     real, dimension(1:nn+1), intent(inout) :: TK, XD, YD   ! reconstructed spline
     integer, intent(in) :: layer
     real, dimension(:), allocatable :: L, D, U, LL, DD, UU
     real :: TS, S,TT,XI,YI,F0,F1,F2,F3
     integer :: i,j,k, n, INFO, NI

     
     
     !DOUBLE PRECISION X(0:MAXN), Y(0:MAXN), XD(MAXN), YD(MAXN)
     ! DOUBLE PRECISION L(MAXN),D(MAXN),U(MAXN)
     ! DOUBLE PRECISION LL(MAXN),DD(MAXN),UU(MAXN)
     ! DOUBLE PRECISION TK(MAXN),TS
     ! DOUBLE PRECISION T,S,TT,XI,YI,F0,F1,F2,F3

     !CHARACTER*128 FNAME

      !PRINT *,'Name of input file:'
      !READ (*,'(A128)') FNAME
      !OPEN (UNIT=20, FILE=FNAME, STATUS='OLD')
      !I=0
!1     !I=I+1
!      READ (20,*,END=3) X(I),Y(I)
!      IF (I.LT.MAXN) GOTO 1
!      STOP 'I .LT. MAXN in GETDATA'
!
!3     CLOSE (20)

     !N=I-1
     ! if (x(n) .eq. x(1) .and. y(n) .eq. y(1)) n = n - 1
     ! if (n .le. 1) stop 'N .le. 1 in SPLINE'

      x(0) = x(nn)
      y(0) = y(nn)
      x(nn+1) = x(1)
      y(nn+1) = y(1)

      N = nn+1

      allocate(  L(1:n), D(1:n), U(1:n), LL(1:n), DD(1:n), UU(1:n) )

      
!  Vypocet tecnych vektoru spline krivky
      TS=0.D0
      do  I=1,N-1
         TK(I)=SQRT((X(I)-X(I+1))**2+(Y(I)-Y(I+1))**2)
         TS=TS+TK(I)
      end do

      DO  I=2,N-1
         D(I)=2.D0/TK(I-1)+2.D0/TK(I)
      enddo

      D(1)=2.D0
      D(N)=2.D0

      DO  I=2,N-1
       L(I)=1.D0/TK(I-1)
      enddo
      L(N)=1.D0

      DO I=2,N-1
         U(I)=1.D0/TK(I)
      enddo
      U(1)=1.D0

      do  I=2,N-1
         XD(I)=3.0/TK(I)**2*X(I+1) &
              -(3.0/TK(I)**2-3.0/TK(I-1)**2)*X(I) &
              -3.0/TK(I-1)**2*X(I-1)
         
         YD(I)=3.0/TK(I)**2*Y(I+1) &
              -(3.0/TK(I)**2-3.0/TK(I-1)**2)*Y(I) &
              -3.0/TK(I-1)**2*Y(I-1)
      end do

      XD(1)=3.0*(X(2)-X(1))/TK(1)
      XD(N)=3.0*(X(N)-X(N-1))/TK(N-1)
      YD(1)=3.0*(Y(2)-Y(1))/TK(1)
      YD(N)=3.0*(Y(N)-Y(N-1))/TK(N-1)

      CALL DCOPY (N,L,1,LL,1)
      CALL DCOPY (N,D,1,DD,1)
      CALL DCOPY (N,U,1,UU,1)
      CALL DGTSL (N,LL,DD,UU,XD,INFO)
      CALL DGTSL (N,L,D,U,YD,INFO)

!  Vypocet Fergusonovych kubik

      !PRINT *,'Number of inserted points (put n+1) :'
      !READ (*,*) NI
      !PRINT *,'Name of output file:'
      !READ (*,'(A128)') FNAME

      !OPEN (20,FILE='smaz_spline', STATUS='UNKNOWN')

      NI = 5
      !write(20+layer,100) x(1),y(1)
      
      do i=1,N-2
         DO  J=1,NI-1
            !  T=FLOAT(J-1)/FLOAT(NI-1)*TS
            s= 1.0 * j / ni * TK(i)
            F0=2.0*(S/TK(I))**3-3.0*(S/TK(I))**2+1
            F1=-2.0*(S/TK(I))**3+3.0*(S/TK(I))**2
            F2=S**3/TK(I)**2-2.0*S**2/TK(I)+S
            F3=S**3/TK(I)**2-S**2/TK(I)
            XI=F0*X(I)+F1*X(I+1)+F2*XD(I)+F3*XD(I+1)
            YI=F0*Y(I)+F1*Y(I+1)+F2*YD(I)+F3*YD(I+1)
            
            !WRITE (20+layer,100) XI,YI
         enddo
         !         if(i .lt. N-1) then
         !write(20+layer,100) x(i+1), y(i+1)
         !            print *,i,n,x(i+1),y(i+1)
         !         endif
      enddo
      !CLOSE (20)


      deallocate(L, D, U, LL, DD, UU) 
      
100   FORMAT (2E16.6)
    end subroutine spline_reconstruct




end module spline_rec
