********************************************************************************
*                                                                              *
*  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 ### interpolq ### 
c
c     This routine works for ngrids<=3
c
      subroutine interpolq
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      _REAL_ vmuq,coeffq,xmu,vpoly1q
      include 'commons8.inc'
      include 'commons16.inc'
c
      dimension coeffq(9)

c  |--------------------------------------|
c  |            ifill=1                   |    ngrids=1 
c  |				          | 
c  |--------------------------------------|
c
c  |--------------|-----------------------|
c  |     ifill=2  | ifill=4               |    ngrids=2
c  |              |                       |
c  |--------------|-----------------------|
c
c  |--------------|-----------|-----------|
c  |	ifill=2	  | ifill=3   | ifill=4	  |    ngrids=3 
c  |              |           |           |
c  |--------------|-----------|-----------|

c     To relax a function defined on nni*nmu grid using the 8th-order central
c     difference formulae the function is immersed into (nni+8)*(nmu+8) grid
c     but the valued outside the boundary necessary to perform relaxation
c     within the boundary region have to be previded. The symmetry properties
c     are used to this end along ni=1, ni=nni and mu=1 lines. However the
c     boundaries between subgrids have to be treated differently.

c     At each subgrid boundary a set of the 8th-order Lagrange interpolation 
c     formulae are formed in order to be able to calculate missing outer
c     region data.

c     For each value of mup (p) a Lagrange polynomial of the 8th-order
c     employing 4 grid points to the left and right of imu (imu is included)
c     is build (x_{0} corresponds to the grid boundary mu value) 

c     f^{(p)}(x) = \sum_{k=p-4}^{p+4} f(x_{k})
c                  \prod_{i=p-4,i\ne k}^{p+4} {(x-x_{i}) \over (x_{k}-x_{i})}

c     Polynomials for each k have to be calculated and stored.

c     cint?(k,i) store the interpolation coefficients; i=1,...,4 corresponds
c     to appropriate (left or right) inter(sub)grid mu values in increasing 
c     order, k numbers the coefficients. 
c     The interpolation formula is constructed using 9 grid points somehow 
c     distributed around the point at which interpolated value is sought 
c     (see below).

c     Determine the grid points which will be used for the construction of the
c     polynomial

c !!!!!!!! iord is changed

      iord=9
      iord2=iord/2
      kbeg=1
      kend=9

      if (idbg(80).ne.0) then
         iord=3
         iord2=iord/2
         kbeg=4
         kend=6
      endif 

      if (ngrids.eq.2) then

c        ifill=2 
c        boundary between 1-2

         do i=1,4
              do k=1,9
                 cint2(k,i)=0.d0
                 cint4(k,i)=0.d0
              enddo
           enddo

           hmu1=hmu(1)
           hmu2=hmu(2)	

           if (hmu2.lt.hmu1) then	  
              write(iout6,*) '2nd subgrid should be sparser than the 1st'
              stop 'interpolq: ifill=2'
           endif	  	

c          for each i=1,...,4 determine the smallest imup for which 
c          |vmu(imup)-vmu(iemu(1))| > i*hmu(1)
c          and construct polynomial around imup+1-4 using 4 grid points to the	
c          left and right

           do i=1,4
              do imu=iemu(1)+1,iemu(2)
                 if ((vmu(imu)-vmu(iemu(1))).ge.dble(i)*hmu1) then	  
                    mup=imu+1-iord2
                    goto 90
                 endif				  	     
              enddo
00090         continue	
              iadint2(i)=mup-4
              do k=kbeg,kend
                 call lpcoeffq(mup,k,coeffq)
                 xmu=dble(vmu(iemu(1))+dble(i)*hmu1)
                 cint2(k,i)=vpoly1q(xmu,coeffq)
              enddo
           enddo				  

c          ifill=4 

c          to the left of the boundary (only 3 values are actually needed but
c          4 values would be generated)
c          for each i=1,...,4 determine the largest imup for which 
c          |vmu(imup)-vmu(iemu(1))| > i*hmu(2)
c          and construct polynomial around such imup using 4 grid points to the
c          left and right

c          boundary between 1-2	

           do i=1,4
              do imu=iemu(1)-1,1,-1
                 if ((vmu(iemu(1))-vmu(imu)).ge.dble(i)*hmu2) then	  
                    mup=imu
                    goto 100
                 endif				  		  
              enddo
00100         continue
              iadint4(5-i)=mup-4
              do k=kbeg,kend
                 call lpcoeffq(mup,k,coeffq)
                 xmu=dble(vmu(iemu(1))-dble(i)*hmu2)
c                reverse order of addressing columns, cf. fill4		
                 cint4(k,5-i)=vpoly1q(xmu,coeffq)
              enddo
           enddo				  

           if (idbg(81).ne.0) then
              write(iout6,*) 'iadint2',iadint2
              write(iout6,*) 'cint2'
              do i=1,4
                 write(iout6,*) i,(cint2(k,i),k=1,9)
              enddo

              write(iout6,*) 'iadint4',iadint4
              write(iout6,*) 'cint4'
              do i=1,4
                 write(iout6,*) i,(cint4(k,i),k=1,9)
              enddo
           endif
           return
	endif	  

c       33333333333333333333333333333333333333333333333
c       ifill=3 (ngrids=3)	

	if (ngrids.eq.3) then

c          boundary between 1-2

           hmu1=hmu(1)
           hmu2=hmu(2)	

           if (hmu2.lt.hmu1) then	  
              write(iout6,*) '2nd subgrid should be sparser than the 1st'
              stop 'interpolq: ifill=2'
           endif

           do i=1,4
                do k=1,9
                   cint2 (k,i)=0.d0
                   cint3l(k,i)=0.d0
                   cint3r(k,i)=0.d0
                   cint4 (k,i)=0.d0
                enddo
             enddo	  	

c            for each i=1,...,4 determine the smallest imup for which 
c            |vmu(imup)-vmu(iemu(1))| > i*hmu(1)
c            and construct polynomial around imup+1-4 using 4 grid points 
c            to the left and right

             do i=1,4
                do imu=iemu(1)+1,iemu(2)
                   if ((vmu(imu)-vmu(iemu(1))).ge.dble(i)*hmu1) then	  
                      mup=imu+1-iord2
                      goto 91
                   endif				  	     
                enddo
00091           continue	
                iadint2(i)=mup-4
                do k=kbeg,kend
                   call lpcoeffq(mup,k,coeffq)
                   xmu=dble(vmu(iemu(1))+dble(i)*hmu1)
                   cint2(k,i)=vpoly1q(xmu,coeffq)
                enddo
             enddo				  

c            ifill=4 

c            to the left of the boundary (only 3 values are actually needed but
c            4 values would be generated)
c            for each i=1,...,4 determine the largest imup for which 
c            |vmu(imup)-vmu(iemu(1))| > i*hmu(2)
c            and construct polynomial around such imup using 4 grid points to the
c            left and right

c            boundary between 2-3	

             hmu3=hmu(3)
             if (hmu3.lt.hmu2) then	  
                write(iout6,*) '3rd subgrid should be sparser than the 2nd'
                stop 'interpolq: ifill=3'
             endif	  	
             
             do i=1,4
                do imu=iemu(2)-1,1,-1
                   if ((vmu(iemu(2))-vmu(imu)).ge.dble(i)*hmu3) then	  
                      mup=imu
                      goto 101
                   endif				  		  
                enddo
00101           continue
                iadint4(5-i)=mup-4
                do k=kbeg,kend
                   call lpcoeffq(mup,k,coeffq)
                   xmu=dble(vmu(iemu(2))-dble(i)*hmu3)
c                  reverse order of addressing columns, cf. fill4		
                   cint4(k,5-i)=vpoly1q(xmu,coeffq)
                enddo
             enddo				  

c            to the left of the boundary (only 3 values are needed but 4 ...)

c            for each i=1,...,4 determine the largest imup for which 
c            |vmu(imup)-vmu(iemu(1))| > i*hmu(2)
c            and construct polynomial around such imup using 4 grid points to the
c            left and right

c            boundary between 1-2	

             do i=1,4
                do imu=iemu(1)-1,1,-1
                   if ((vmu(iemu(1))-vmu(imu)).ge.dble(i)*hmu2) then	  
                      mup=imu
                      goto 110
                   endif				  		  
                enddo
00110           continue
                iadint3l(5-i)=mup-4
                do k=kbeg,kend
                   call lpcoeffq(mup,k,coeffq)
                   xmu=dble(vmu(iemu(1))-dble(i)*hmu2)
c                  reverse order of addressing columns, cf. fill3
                   cint3l(k,5-i)=vpoly1q(xmu,coeffq)
                enddo
             enddo				  

c            boundary between 2-3

c            for each i=1,...,4 determine the smallest imup for which 
c            |vmu(imup)-vmu(iemu(2))| > i*hmu(2)
c            and construct polynomial around imup+1-4 using 4 grid points to the  	
c            left and right

             do i=1,4
                do imu=iemu(2)+1,iemu(3)
                   if ((vmu(imu)-vmu(iemu(2))).ge.dble(i)*hmu2) then	  
                      mup=imu+1-iord2
                      goto 120
                   endif
                enddo
00120           continue	
                iadint3r(i)=mup-4
                do k=kbeg,kend
                   call lpcoeffq(mup,k,coeffq)
                   xmu=dble(vmu(iemu(2))+dble(i)*hmu2)
	      cint3r(k,i)=vpoly1q(xmu,coeffq)
           enddo
        enddo	

        if (idbg(81).ne.0) then
           write(iout6,*) 'iadint2',iadint2
           write(iout6,*) 'cint2'
           do i=1,4
              write(iout6,*) i,(cint2(k,i),k=1,9)
           enddo

           write(iout6,*) 'iadint3l',iadint3l
           write(iout6,*) 'cint3l'
           do i=1,4
              write(iout6,*) i,(cint3l(k,i),k=1,9)
           enddo

           write(iout6,*) 'iadint3r',iadint3r
           write(iout6,*) 'cint3r'
           do i=1,4
              write(iout6,*) i,(cint3r(k,i),k=1,9)
           enddo

           write(iout6,*) 'iadint4',iadint4
           write(iout6,*) 'cint4'
           do i=1,4
              write(iout6,*) i,(cint4(k,i),k=1,9)
           enddo
	endif
        return
      endif	  

      if (ngrids.gt.3) then
         write(iout6,*) 'No more than 3 subgrids are allowd'
         stop 'interpolq'
      endif
      return
      end

c ### vpoly1q ###
c  
c     This function uses the Horner scheme to calculate value of the polynomial 
c     stored in array a at a particular point
c
      function vpoly1q (x,a)
      implicit _REAL_ (a-h,o-z)
      _REAL_ vmuq
      include 'commons16.inc'
      dimension a(9)

      vpoly1q=0.d0
      do i=iord,2,-1
         vpoly1q=(vpoly1q+a(i))*x
      enddo
      vpoly1q=vpoly1q+a(1)
      return
      end

c ### lagrpolq ###
c
      subroutine lagrpolq (dc1,dc2)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      _REAL_ vmuq,coeffq,coeff1,coeff2,vpolyq
      include 'commons8.inc'
      include 'commons16.inc'

      dimension dc1(9,7,9),dc2(9,7,9),coeffq(9),coeff1(9),coeff2(9)

c     derivative coeff. from 8th-order Lagrange interpolation formula
c     are stored in dc[1-2](ngbound,imu,k), where ngbound is the number of 
c     the grid boundaries (1 for 1-2, 2 for 2-3 etc), imu is one of the 7 
c     intergrid mu values and k points to a derivative coefficient.

c     the grid boundaries (1 for 1-2, 2 for 2-3 etc), imu is one of the 7 
c     intergrid mu values and k points to a derivative coefficient.

c     For each value of mup (p) a Lagrange polynomial of the 8th-order
c     employing 4 grid points to the left and right of imu (imu is included)
c     is build (x_{0} corresponds to the grid boundary mu value) 

c     f^{(p)}(x) = \sum_{k=p-4}^{p+4} f(x_{k})
c                  \prod_{i=p-4,i\ne k}^{p+4} {(x-x_{i}) \over (x_{k}-x_{i})}

c     Polynomials for each k have to be calculated and their first and second
c     derivatives obtained. The values of the derivative polynomials at grid 
c     point x_{p} (mup) are the differentiation coefficients sought.

c     in order to quarantee the real*8 in the Lagrange coefficients
c     the quadruple precision have to be used during their generation

      iord=9
      kbeg=1
      kend=9

      do i=1,mxnmu
         vmuq(i)=dble(vmu(i))
      enddo

c     loop over grid boundaries first (ngrids>1)

      do ig=1,ngrids-1
         ib=iemu(ig)-3	
         ie=iemu(ig)+3
c        loop over mup values (mup=ib...ie)
c        loop over k	     
         imup=0
         do mup=ib,ie
            imup=imup+1
	    do k=kbeg,kend
               call lpcoeffq(mup,k,coeffq)
               call lpderq(coeffq,coeff1,coeff2)
               dc1(ig,imup,k)=vpolyq(mup,coeff1)
               dc2(ig,imup,k)=vpolyq(mup,coeff2)
            enddo 
         enddo
      enddo
      return
      end

c ### lpcoeffq ###
c
c     This routine calculates coefficients of the (sub)Lagrange 
c     polynomial for a grid point k
c     \prod_{i=1,i\ne k}^{9} {(x-x_{i}) \over (x_{k}-x_{i})}
c
c     x_{k}= vmu(mup-5+k), k=1,...,9 
c
c     Last modification: 2.01.01
c     
      subroutine lpcoeffq (mup,k,coeffq)
      implicit _REAL_ (a-h,o-z)
      include 'commons16.inc'
      dimension a(12),b(12),coeffq(9)

c  calculate denominator product

      denom=1.d0
      ib=mup-(iord/2+1)
      do i=kbeg,kend
         if (i.ne.k) denom=denom*(vmuq(ib+k)-vmuq(ib+i))
      enddo

      do i=1,12
         a(i)=0.d0
         b(i)=0.d0	 
      enddo

c     calculate nominator product: 
c     a(1)+a(2)*x+a(3)*x^2+...+a(9)*x^8
c     storing coefficients in a and b

c     if (k.ne.kbeg) then
c     a(1)=-vmuq(ib+kbeg)
c     else
c     a(1)=-vmuq(ib+kbeg+1)	  
c     endif
c     a(2)=1.d0	

      a(1)=1.d0

c  multiply polynomial a by (x+c1), c1=-vmuq(ib+i)

      ic2=1	
      do i=kbeg,kend	  
         if (i.ne.k) then
  	    ic2=ic2+1
            c1=-vmuq(ib+i)	  
	    b(1)=a(1)*c1
	    do ic1=2,ic2
               b(ic1)=a(ic1)*c1+a(ic1-1)
	    enddo
            b(ic2+1)=a(ic2)	
            do j=1,ic2+1
               a(j)=b(j)
            enddo
         endif	    
      enddo   
c     
      do i=1,iord
         coeffq(i)=a(i)/denom
      enddo
      return
      end

c ### lpderq ###
c     
c     This routine calculates coefficients of the first and second
c     derivative of the polynomial stored in a
c
      subroutine lpderq (a,a1,a2)
      implicit _REAL_ (a-h,o-z)
      include 'commons16.inc'
      dimension a(9),a1(9),a2(9)

      do i=1,iord-1
         a1(i)=a(i+1)*dble(i)
      enddo
      a1(iord)=0.d0

      do i=1,iord-2
         a2(i)=a1(i+1)*dble(i)
      enddo
      a2(iord-1)=0.d0
      a2(iord)  =0.d0
      return
      end

c ### vpolyq ###
c
c     This function uses the Horner scheme to calculate value of the polynomial 
c     stored in array a at a particular point
c
      function vpolyq (mup,a)
      implicit _REAL_ (a-h,o-z)
      include 'commons16.inc'
      dimension a(9)

      x=vmuq(mup)  
      vpolyq=0.d0
      do i=iord,2,-1
         vpolyq=(vpolyq+a(i))*x
      enddo
      vpolyq=vpolyq+a(1)
      return
      end