      subroutine t2(embed_in,pair_in,den_in,rad,ndata)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Transformation from tabulated file format.
c     Only the differentiation of the function is added
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision (a-h,o-z)
      parameter (np=10000)
      dimension ax(np),bx(np),cx(np)
      dimension ay(np),by(np),cy(np)
      dimension da(np),db(np),dc(np)
      character trash
      character*80 embed_in,pair_in,den_in
      character*20 embed_out,pair_out,den_out
c-----------------------------------------------------------------------
      embed_out='out.embed' !output files
      pair_out='out.pair'
      den_out='out.den'
c-----------------------------------------------------------------------
! reads 1-st potential embed
      open(1,file=embed_in)
      do i=1,2
         read (1,*) trash
      end do
      do i=1,ndata
         read (1,*) ax(i),ay(i)
      end do
      close (1)
! reads 2-nd potential pair
      open(1,file=pair_in)
      do i=1,2
         read (1,*) trash
      end do
      do i=1,ndata
         read (1,*) bx(i),by(i)
      end do
      close (1)
! reads 3-rd potential rho
      open(1,file=den_in)
      do i=1,2
         read (1,*) trash
      end do
      do i=1,ndata
         read (1,*) cx(i),cy(i)
      end do
      close (1)
c ----------------------------------------------------------------------
! differentiation
      call derivacia2(ax,ay,ndata,da)
      call derivacia2(bx,by,ndata,db)
      call derivacia2(cx,cy,ndata,dc)
! output files
      open (2,file=embed_out)
      do i=1,ndata
         write (2,10) ax(i),ay(i),da(i)
      end do
      close (2)
      open (2,file=pair_out)
      do i=1,ndata
         write (2,10) bx(i),by(i),db(i)
      end do
      close (2)
      open (2,file=den_out)
      do i=1,ndata
         write (2,10) cx(i),cy(i),dc(i)
      end do
      close (2)
      rad=bx(ndata)
   10 format (F14.8,E24.9,E24.9)
      end
c
c ============================ subroutines =============================
c
      subroutine derivacia2 (x,y,n,dy)
      implicit double precision (a-h,o-z)
      dimension dy(n),x(n),y(n)
      do i=1,n
         dy(i)=0.0d0
      end do
      yp1=(y(2)-y(1))/(x(2)-x(1))
      ypn=(y(n)-y(n-1))/(x(n)-x(n-1))
      call spline2 (x,y,n,yp1,ypn,dy)
      end
c
c
c============================== downloaded subroutines =================
c
       SUBROUTINE spline2(x,y,n,yp1,ypn,y2)
! Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., yi = f(xi), with
! x1 < x2 < ... < xN, and given values yp1 and ypn for the first derivative of the interpolating
! function at points 1 and n, respectively, this routine returns an array y2(1:n) of
! length n which contains the second derivatives of the interpolating function at the tabulated
! points xi. If yp1 and/or ypn are equal to 1x10^30 or larger, the routine is signaled to set
! the corresponding boundary condition for a natural spline2, with zero second derivative on
! that boundary. Parameter: NMAX is the largest anticipated value of n.
!
! Please notice that the program spline2 is called only once to process an entire tabulated
! function in arrays xi and yi. Once this has been done, the values of the interpolated function
! for any value of x are obtained by calls (as many as desired) to a separate routine splint.
       INTEGER n,NMAX
       REAL*8 yp1,ypn,x(n),y(n),y2(n)
       PARAMETER (NMAX=500)
       INTEGER i,k
       REAL*8 p,qn,sig,un,u(n)
       if (yp1.gt..99e30) then
       y2(1)=0.
       u(1)=0.
       else
       y2(1)=-0.5
       u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
       endif
       do 11 i=2,n-1
       sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
       p=sig*y2(i-1)+2.
       y2(i)=(sig-1.)/p
       u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
   11 continue
       if (ypn.gt..99e30) then
       qn=0.
       un=0.
       else
       qn=0.5
       un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
       endif
       y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
       do 12 k=n-1,1,-1
       y2(k)=y2(k)*y2(k+1)+u(k)
   12 continue
       return
       END
c ======================================================================
c
c                     Last modified: 18.03.2013              
c                                 
c                 E. Vitkovska (atavaha@hotmail.com)
c                                
c=======================================================================
