********************************************************************************
*                                                                              *
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### dointerp_mu ###	
c
c     Performs interpolations of functions in mu variable.
c
      subroutine dointerp_mu (nnit,nmumin_p,nmumax_p,nmumin,nmumax,
     &     fbefore,fafter) 

      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)	
      real*8 vmuq,coeffq(9),coeffq2(9,9),xmu,vpoly1q  
      include 'commons8.inc'
      include 'commons16.inc'
      dimension fbefore(nnit,*),fafter(nnit,*)

      rerror=2
      idebug1=0
      idebug2=0
      idebug3=0

c     interpolation for the first iord2 points in mu variable

      do imu=nmumin,nmumax
         if(vmu(imu).ge.vmu_p((nmumin_p-1)+iord2)) goto 115
      enddo
00115 continue
      nmu_first=imu
 
c     print *,"dointerp_mu: nmumin,nmu_first",nmumin,nmu_first
      
      do k=1,kend
         call lpcoeffq((nmumin_p-1)+iord2+1,k,coeffq)
         do i=1,kend
            coeffq2(i,k)=coeffq(i)
         enddo
      enddo
      do imu=nmumin,nmu_first
         xmu=vmu(imu)            
         do ini=1,nnit
            fafter(ini,imu)=0.d0
            do k=1,kend
               fafter(ini,imu)=fafter(ini,imu)+
     &              fbefore(ini,nmumin_p-1+k)
     &              *vpoly1q(xmu,coeffq2(1,k))
            enddo
            if (idebug1.eq.1) then 
               if (abs(fafter(ini,imu)-fbefore(ini,nmumin_p+2)).gt.
     &              abs(fbefore(ini,nmumin_p+2))*rerror) then
                  
                  write(*,'(2i5,d15.3,4x,3d15.3)') 
     &                 ini,imu,fafter(ini,imu),
     &                 (fbefore(ini,nmumin_p-1+k),k=1,kend)
               endif
            endif
            
         enddo
      enddo

c     Interpolation for the last iord2 points in mu variable.
c     Determine the location of the tail region in the new grid

      do imu=nmumin,nmumax
         if(vmu(imu).ge.vmu_p(nmumax_p)) goto 120
      enddo
00120 continue
      nmu_last=imu
      
      do k=1,kend
         call lpcoeffq(nmumax_p-iord2,k,coeffq)
         do i=1,kend
            coeffq2(i,k)=coeffq(i)
         enddo
      enddo
      do imu=nmu_last-iord2+1,nmumax
         xmu=vmu(imu)            
         do ini=1,nnit
            fafter(ini,imu)=0.d0
            do k=1,kend
               fafter(ini,imu)=fafter(ini,imu)+
     &              fbefore(ini,nmumax_p-iord+k)
     &              *vpoly1q(xmu,coeffq2(1,k))
            enddo
            if (idebug2.eq.1) then 
               if (abs(fafter(ini,imu)-fbefore(ini,nmumax_p-iord+2)).gt.
     &              abs(fbefore(ini,nmumax_p-iord+2))*rerror) then
                  
                  write(*,'(2i5,d15.3,4x,3d15.3)') 
     &                 ini,imu,fafter(ini,imu),
     &                 (fbefore(ini,nmumax_p-iord+k),k=1,kend)
               endif
            endif
         enddo
      enddo

c     interpolation for the inner points in this region
      
      do imu=nmu_first+1,nmu_last-iord2
         xmu=vmu(imu)            
         do i=1,nmumax_p
            if(vmu_p(i).ge.xmu) goto 130
         enddo
00130    continue
         imu_p=i
         do k=1,kend
            call lpcoeffq(imu_p,k,coeffq)
            do i=1,kend
               coeffq2(i,k)=coeffq(i)
            enddo
         enddo
         
         do ini=1,nnit
            fafter(ini,imu)=0.d0
            do k=1,kend
               fafter(ini,imu)=fafter(ini,imu)+
     &              fbefore(ini,imu_p-iord2-1+k)
     &              *vpoly1q(xmu,coeffq2(1,k))
            enddo
            if (idebug3.eq.1) then
               if (abs(fafter(ini,imu)-fbefore(ini,imu_p-iord2+1)).gt.
     &              abs(fbefore(ini,imu_p-iord2+1))*rerror) then
                  write(*,'(2i5,d15.3,4x,3d15.3,2i5)') 
     &                 ini,imu,fafter(ini,imu),
     &                 (fbefore(ini,imu_p-iord2-1+k),k=1,kend),
     &                 nmu_first,nmu_last
               endif
            endif
         enddo
      enddo
      return
      end