      program intpoldp
c INTerPOLating radial Density Profile data (in output log of ev8xi.f)
c (y/m/d) 
c 11/11/17 intpoldp.v1.f reads from log file. This new version can read
c          only data lines.
c 11/11/13 intpoldp.f was created by modifying intpolsq.for (N.Tajima)
c 92/08/14 intpolsq.for was created (N.Tajima)
      implicit real*8 (a-h,o-z)
c     parameter(ndstdi=5,ndstdo=6)
      parameter(maxi=4096)
      character buff*132,flin*60,flout*60
      dimension xa(maxi),ya(maxi,3)
      data ndi,ndo / 12, 13 /, idata0 /5/
c
      write(6,*) 'input data file name ? :'
      read(5,'(A)') flin
      open(ndi,file=flin,status='OLD',form='FORMATTED')
      i=idata0
      do
        read(ndi,'(A)',end=5,err=5) buff
          i=i+1
          if(i.gt.maxi) stop 1
          read(buff,*) xa(i),ya(i,1),ya(i,2),ya(i,3),npts
          if(xa(i).le.0.0d0 .or. xa(i).gt.1.0d+3 .or. npts.lt.0) then
             write(6,*) 'seems wrong input data file',i,xa(i),npts
             stop 2
          end if
      end do
    5 continue
      close(ndi)
      imax=i
      write(6,'(A,I5)') 'number of data lines = ',imax-idata0
      write(6,'(A,F9.5)') 'maximum radius = ',xa(imax)
      do i=1,idata0
        j=2*idata0-i+1 ! density is extended to negative radius
        xa(i)=-xa(j)   ! assuming it is an even function
        ya(i,1)=ya(j,1) ! for the sake of interpolation
        ya(i,2)=ya(j,2)
        ya(i,3)=ya(j,3)
      end do
c
      write(6,*) 'output file name ? :'
      read(5,'(A)') flout
      open(ndo,file=flout,status='UNKNOWN',form='FORMATTED')
c
      write(6,'(A,F9.5)') 'maximum output radius  (fm) ?   must be <='
     & ,xa(imax)
      read(5,*) xmax
      if(xmax .gt. xa(imax)) then
        write(6,'(F9.5,A,F9.5)') xmax,' is replaced with ',xa(imax)
        xmax=xa(imax)
      end if
      write(6,*) 'spacing of output radius grid (fm) ? '
      read(5,*) dx
c
      write(6,'(A,F9.5,A,F9.5)') '(maxinum r)=',xmax
     &                  ,'   (increment of r)=',dx
c
      kmax=nint(xmax/dx)
      do k=0,kmax
        x=k*dx
        call intpol(xa,ya(1,1),imax,x, y1)
        call intpol(xa,ya(1,2),imax,x, y2)
        call intpol(xa,ya(1,3),imax,x, y3)
c       err=y3-(y1+y2)
c       write(6,'(F9.5,1PE15.5)') x,err ! error is small
        write(ndo,'(F7.4,1P3E16.7)') x,y1,y2,y3
      end do
      close(ndo)
      end
c---------------------------------------------------------------------
      subroutine intpol(xall,yall,imax,x,y)
      implicit real*8 (a-h,o-z)
      parameter(nd=4)
      dimension xall(imax),yall(imax),xa(nd),ya(nd)
c
      i1=imax
      do i=1,imax
        if(xall(i).gt.x) then
          i1=i
          exit
        end if
      end do
c
      ndhf=nd/2
      ib=i1-ndhf
      ie=i1+ndhf-1
      if(ib.lt.1) then
        ib=1
        ie=ib+2*ndhf-1
      end if
      if(ie.gt.imax) then
        ie=imax
        ib=ie-2*ndhf+1
      end if
      ib=max(1   ,ib)
      ie=min(imax,ie)
      np=ie-ib+1
      do 20 i=ib,ie
      xa(i-ib+1)=xall(i)
      ya(i-ib+1)=yall(i)
   20 continue
      y=polylg(np,xa,ya,x)
      end
c--------------------------------------------------------------------c
C POLYLG : LEGENDRE POLYNOMIAL INTERPOLATION:90:8.3-90.8.3
      FUNCTION POLYLG(NA,XP,YP,X0)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(MXN=20,RMIN=1.0D-37)
      DIMENSION XP(*),YP(*)
      DIMENSION X(MXN),Y(MXN),C(MXN)
      SAVE
      DATA INICAL / 0 /
C
      N=NA
C
      IF(N.GT.0) THEN
        IF(N.GT.MXN) THEN
          WRITE(6,*) 'POLYLG:N=',N,'-->',MXN
          N=MXN
        END IF
        DO 10 I=1,N
          X(I)=XP(I)
          Y(I)=YP(I)
   10   CONTINUE
        DO 40 I=1,N
          S=1.0D0
          DO 20 J=1,I-1
            S=S*(X(I)-X(J))
   20     CONTINUE
          DO 30 J=I+1,N
            S=S*(X(I)-X(J))
   30     CONTINUE
          IF(S.EQ.(0.0D0)) THEN
            WRITE(6,*) 'POLYLG:SOME OF THE GIVEN X''S ARE THE SAME.'
            POLYLG=0.0D0
            INICAL=0
            RETURN
          END IF
          C(I)=Y(I)/S
   40   CONTINUE
        INICAL=1
      END IF
C
      IF(INICAL.EQ.0) THEN
        WRITE(6,*) 'POLYLG:N=0 WHEN POINTS HAVE NOT BEEN GIVEN:Y=0'
        POLYLG=0.0D0
        RETURN
      END IF
C
      S=0.0D0
      T=1.0D0
      DO 50 I=1,N
        DX=X0-X(I)
        IF((DX.GT.(-RMIN)).AND.(DX.LT.RMIN)) THEN
          POLYLG=Y(I)
          RETURN
        END IF
        S=S+C(I)/DX
        T=T*DX
   50 CONTINUE
      POLYLG=S*T
C
      RETURN
      END
