text:
cc Copyright (C) 2017: Travis Askham, Leslie Greengard,
cc Zydrunas Gimbutas
cc email: askhamwhat@gmail.com
cc
cc This software is being released under a modified FreeBSD license
cc (see licenses folder in home directory).

c      This file contains the basic subroutines for
c      forming and evaluating multipole (partial wave) expansions
c      in two dimensions.
c
c-----------------------------------------------------------------------
c         **************  IMPORTANT NOTE **************
c
c      THE CONVENTION IN THIS SUBROUTINE LIBRARY IS THAT ALL COMPUTATIONS
C      USING EXPANSIONS ARE INCREMENTAL.
c      That is, the multipole and local expansions
c      on output are incremented from their input values, not overwritten.
c      The same is true for translation routines.
c
c-----------------------------------------------------------------------
c
c      MBH2DCONVTOMP_VEC: converts charge, dipole, quadrupole, octopole
c      			  sources to multipolar
c
c      MBH2DFORMMPMP_VEC:  creates multipole expansions (outgoing) due to
c                       a collection of multipolar sources.
c      MBH2DMPEVALP_VEC:  computes potentials from multipole expansions
c                       at a collection of targets
c      MBH2DMPEVALG_VEC:  computes potentials/gradients from multipole
c                          expansions at a collection of targets
c      MBH2DMPEVALH_VEC:  computes potentials/gradients/Hessians
c                         from multipole expansions at a collection of
c                         targets
c      MBH2DFORMTAMP_VEC:  creates local expansions due to
c                       a collection of multipolar sources
c      MBH2DTAEVALP_VEC:  computes potentials from local expansions
c                       at a collection of targets
c      MBH2DTAEVALG_VEC:  computes potentials/gradients from local expansions
c                       at a collection of targets
c      MBH2DTAEVALH_VEC:  computes potentials/gradients/Hessians
c                       from local expansions at a collection of targets
c
c      MBH2DMPMP_VEC:     Translates center of multipole expansions
c      MBH2DLOCLOC_VEC:   Translates center of local expansions
c      MBH2DMPLOC_VEC:    Converts multipole expansions to local expansions
c
c     MBH2DMPZERO_VEC:   utility to initialize multipole coefficients
c                           to zero.
C
c-----------------------------------------------------------------------

      subroutine mbh2dconvtomp_vec(nd,beta,ns,ifcharge,charge,
     1     ifdipole,dipstr,dipvec,ifquad,quadstr,quadvec,ifoct,octstr,
     2     octvec,nterms,mbhmpole,ympole)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Converts up to octopole strength sources to a short, equivalent
c     multipolar expansion
c
c     input:
c
c     nd - integer, number of expansions at each source location
c     beta - real *8, modified biharmonic parameter
c     ns - integer, number of sources
c     ifcharge - integer, flag, equals 1 if charges are present
c     ifdipole - integer, flag, equals 1 if dipoles are present
c     ifquad - integer, flag, equals 1 if quadrupoles are present
c     ifoct - integer, flag, equals 1 if octopoles are present
c     charge - real *8 (ns) charge strengths
c     nterms - integer, 0:nterms is the second dimension of
c                the output mbhmpole,ympole arrays
c                requirements:
c            for charge nterms >= 0
c            for dipole nterms >= 1
c            for quadrupoles nterms >= 2
c            for octopoles nterms >= 3
c
c     output:
c
c     NOTE: these arrays are *incremented* on output
c
c     mbhmpole - complex *16 (nd,0:nterms,ns), difference kernel part
c                 of equivalent multipolar expansions
c     ympole - complex *16 (nd,0:nterms,ns), Yukawa part of equivalent
c                 multipolar expansions
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global variables
      integer ns, nterms, ifcharge, ifdipole, ifquad, ifoct, nd
      real *8 beta, charge(nd,*), dipstr(nd,*)
      real *8 dipvec(nd,2,*)
      real *8 quadstr(nd,*), quadvec(nd,3,*), octstr(nd,*)
      real *8 octvec(nd,4,*)
      complex *16 mbhmpole(nd,0:nterms,*), ympole(nd,0:nterms,*)
c     local variables
      integer i, j
      complex *16 ima
      real *8 pi
c
      data ima/(0.0d0,1.0d0)/
c

      pi = 4.0d0*datan(1.0d0)

      do i = 1,ns

         if (ifcharge .eq. 1) then
            do j = 1,nd
               mbhmpole(j,0,i) = mbhmpole(j,0,i) +
     1              charge(j,i)/(2.0d0*pi*beta**2)
            enddo
         endif

         if (ifdipole .eq. 1) then
            do j = 1,nd
               mbhmpole(j,1,i) = mbhmpole(j,1,i) +
     1              dipstr(j,i)*(dipvec(j,1,i) +
     1              ima*dipvec(j,2,i))/(2.0d0*pi*beta)
            enddo
         endif

         if (ifquad .eq. 1) then
            do j = 1,nd
               mbhmpole(j,2,i) = mbhmpole(j,2,i) +
     1              quadstr(j,i)*(quadvec(j,1,i) +
     1              quadvec(j,2,i)*ima -
     2              quadvec(j,3,i))/(4.0d0*pi)
               ympole(j,0,i) = ympole(j,0,i) +
     1              quadstr(j,i)*(quadvec(j,1,i) +
     1              quadvec(j,3,i))/(4.0d0*pi)
            enddo
         endif

         if (ifoct .eq. 1) then
            do j = 1,nd
               mbhmpole(j,3,i) = mbhmpole(j,3,i) +
     1              beta*octstr(j,i)*(octvec(j,1,i) +
     1              ima*octvec(j,2,i) - octvec(j,3,i) -
     2              ima*octvec(j,4,i))/(8.0d0*pi)
               ympole(j,1,i) = ympole(j,1,i) +
     1              beta*octstr(j,i)*(3.0d0*octvec(j,1,i) +
     1              ima*octvec(j,2,i) + octvec(j,3,i) +
     2              3.0d0*ima*octvec(j,4,i))/(8.0d0*pi)
            enddo
         endif

      enddo

      return
      end

      subroutine mbh2dmpevalp_vec(nd,beta,rscale,center,
     1     mbhmpole,ympole,nterms,ztarg,ntarg,
     2     pot)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     beta               : Yukawa parameter
c     rscale             : the scaling factor of the expansions
c     center(2)          : expansion center
c     mbhmpole(0:nterms) : coefficients for difference-type functions
c     ympole(nterms)     : coefficients for modified Bessel functions
c     nterms             : number of terms in expansions
c     ztarg(2,ntarg)     : coordinates of targets
c     ntarg              : number of targets
c
c     OUTPUT:
c
c     pot(ntarg)      : value of potential at targets
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhmpole(nd,0:nterms), ympole(nd,0:nterms)
      real *8 rscale, center(2), ztarg(2,ntarg), beta
      real *8 pot(nd,ntarg)
      integer nterms, ntarg, nd
c     local
      real *8 zdiff(2), r, theta
      real *8 ders(0:1), pih
      real *8, allocatable :: diffs(:), kvec(:)
      complex *16 z, eye, ztemp1, ztemp2
      complex *16, allocatable :: mptemp1(:),mptemp2(:)
      real *8 dc, ds, dc2, ds2, dcx, dsx, dcy, dsy
      integer ifder, ifders, l, j, i
      data eye /(0.0d0,1.0d0)/

      allocate(mptemp1(0:nterms+6),mptemp2(0:nterms+6),
     1     diffs(0:nterms+6),kvec(0:nterms+6))

      do i = 1,ntarg
         zdiff(1) = ztarg(1,i)-center(1)
	 zdiff(2) = ztarg(2,i)-center(2)
      	 call h2cart2polar(zdiff,r,theta)

c     get values of difference functions
         ifders = 0
         call diffslogbk_fast(r,beta,rscale,diffs,ifders,ders,kvec,
     1        nterms+2)

         mptemp1(0)=diffs(0)
         mptemp2(0)=kvec(0)
         ztemp2=exp(eye*theta)
         ztemp1=ztemp2
	 do j=1,nterms+2
   	    mptemp1(j)=dcmplx(diffs(j)*dreal(ztemp1),
     1           diffs(j)*dimag(ztemp1))
	    mptemp2(j)=dcmplx(kvec(j)*dreal(ztemp1),
     1           kvec(j)*dimag(ztemp1))
            ztemp1 = ztemp1*ztemp2
	 enddo

c     evaluate
         do j = 1,nd
            pot(j,i) = pot(j,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpole(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympole(j,0))
         enddo

         do l = 1,nterms
            do j=1,nd
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpole(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpole(j,l))
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp2(l))*dreal(ympole(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympole(j,l))
            enddo
         enddo
      enddo

      return
      end

      subroutine mbh2dmpevalg_vec(nd,beta,rscale,center,
     1     mbhmpole,ympole,nterms,ztarg,ntarg,
     2     pot,grad)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     beta               : Yukawa parameter
c     rscale             : the scaling factor of the expansions
c     center(2)          : expansion center
c     mbhmpole(0:nterms) : coefficients for difference-type functions
c     ympole(nterms)     : coefficients for modified Bessel functions
c     nterms             : number of terms in expansions
c     ztarg(2,ntarg)     : coordinates of targets
c     ntarg              : number of targets
c
c     OUTPUT:
c
c     pot(ntarg)      : value of potential at targets
c     grad(2,ntarg)   : value of gradient at targets
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhmpole(nd,0:nterms), ympole(nd,0:nterms)
      real *8 rscale, center(2), ztarg(2,ntarg), beta
      real *8 pot(nd,ntarg)
      real *8 grad(nd,2,ntarg)
      integer nterms, ntarg, nd
c     local
      real *8 zdiff(2), r, theta
      real *8 ders(0:1), pih
      real *8, allocatable :: diffs(:), kvec(:)
      complex *16 z, eye, ztemp1, ztemp2
      complex *16, allocatable :: mptemp1(:),mptemp2(:)
      complex *16, allocatable :: ympolex(:,:), ympoley(:,:)
      complex *16, allocatable :: mbhmpolex(:,:), mbhmpoley(:,:)
      real *8 dc, ds, dc2, ds2, dcx, dsx, dcy, dsy
      integer ifder, ifders, l, j, i
      data eye /(0.0d0,1.0d0)/
c
      allocate(ympolex(nd,0:nterms+1),ympoley(nd,0:nterms+1),
     1     mbhmpolex(nd,0:nterms+1),mbhmpoley(nd,0:nterms+1))
      do i=0,nterms+1
         do j = 1,nd
            ympolex(j,i)=0
            ympoley(j,i)=0
            mbhmpolex(j,i)=0
            mbhmpoley(j,i)=0
         enddo
      enddo

      do j=1,nd
         ympolex(j,1) = -beta/rscale*ympole(j,0)
         ympoley(j,1) = -beta/rscale*ympole(j,0)*eye

         mbhmpolex(j,1) = -beta/rscale*mbhmpole(j,0)
         mbhmpoley(j,1) = -beta/rscale*mbhmpole(j,0)*eye
      enddo

      do i=1,nterms
         do j = 1,nd
            ympolex(j,i-1) = ympolex(j,i-1)
     1           -beta/2.0d0*ympole(j,i)*rscale
            ympolex(j,i+1) = ympolex(j,i+1)
     1           -beta/2.0d0*ympole(j,i)/rscale
            ympoley(j,i-1) = ympoley(j,i-1)
     1           +beta/2.0d0*ympole(j,i)*rscale*eye
            ympoley(j,i+1) = ympoley(j,i+1)
     1           -beta/2.0d0*ympole(j,i)/rscale*eye

            ympolex(j,i-1) = ympolex(j,i-1)
     1           -beta/2.0d0*mbhmpole(j,i)*rscale

            mbhmpolex(j,i+1) = mbhmpolex(j,i+1)
     1           -beta/2.0d0*mbhmpole(j,i)/rscale
            ympoley(j,i-1) = ympoley(j,i-1)
     1           +beta/2.0d0*mbhmpole(j,i)*rscale*eye
            mbhmpoley(j,i+1) = mbhmpoley(j,i+1)
     1           -beta/2.0d0*mbhmpole(j,i)/rscale*eye
         enddo
      enddo

      allocate(mptemp1(0:nterms+6),mptemp2(0:nterms+6),
     1     diffs(0:nterms+6),kvec(0:nterms+6))

      do i = 1,ntarg
         zdiff(1) = ztarg(1,i)-center(1)
	 zdiff(2) = ztarg(2,i)-center(2)
      	 call h2cart2polar(zdiff,r,theta)

c     get values of difference functions
         ifders = 0
         call diffslogbk_fast(r,beta,rscale,diffs,ifders,ders,kvec,
     1        nterms+2)

         mptemp1(0)=diffs(0)
         mptemp2(0)=kvec(0)
         ztemp2=exp(eye*theta)
         ztemp1=ztemp2
	 do j=1,nterms+2
   	    mptemp1(j)=dcmplx(diffs(j)*dreal(ztemp1),
     1           diffs(j)*dimag(ztemp1))
	    mptemp2(j)=dcmplx(kvec(j)*dreal(ztemp1),
     1           kvec(j)*dimag(ztemp1))
            ztemp1 = ztemp1*ztemp2
	 enddo

c     evaluate
         do j = 1,nd
            pot(j,i) = pot(j,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpole(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympole(j,0))
         enddo

         do l = 1,nterms
            do j=1,nd
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpole(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpole(j,l))
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp2(l))*dreal(ympole(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympole(j,l))
            enddo
         enddo
         do j = 1,nd
            grad(j,1,i) = grad(j,1,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpolex(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympolex(j,0))
            grad(j,2,i) = grad(j,2,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpoley(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympoley(j,0))
         enddo
         do l = 1,nterms+1
            do j = 1,nd
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpolex(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpolex(j,l))
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp2(l))*dreal(ympolex(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympolex(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpoley(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpoley(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp2(l))*dreal(ympoley(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympoley(j,l))
            enddo
         enddo
      enddo

      return
      end

      subroutine mbh2dmpevalh_vec(nd,beta,rscale,center,
     1     mbhmpole,ympole,nterms,ztarg,ntarg,
     2     pot,grad,hess)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     beta               : Yukawa parameter
c     rscale             : the scaling factor of the expansions
c     center(2)          : expansion center
c     mbhmpole(0:nterms) : coefficients for difference-type functions
c     ympole(nterms)     : coefficients for modified Bessel functions
c     nterms             : number of terms in expansions
c     ztarg(2,ntarg)     : coordinates of targets
c     ntarg              : number of targets
c
c     OUTPUT:
c
c     pot(ntarg)      : value of potential at targets
c     grad(2,ntarg)   : value of gradient at targets
c     hess(3,ntarg)   : value of Hessian at targets
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhmpole(nd,0:nterms), ympole(nd,0:nterms)
      real *8 rscale, center(2), ztarg(2,ntarg), beta
      real *8 pot(nd,ntarg)
      real *8 grad(nd,2,ntarg)
      real *8 hess(nd,3,ntarg)
      integer nterms, ntarg, nd
c     local
      real *8 zdiff(2), r, theta
      real *8 ders(0:1), pih
      real *8, allocatable :: diffs(:), kvec(:)
      complex *16 z, eye, ztemp1, ztemp2
      complex *16, allocatable :: mptemp1(:),mptemp2(:)
      complex *16, allocatable :: ympolex(:,:), ympoley(:,:)
      complex *16, allocatable :: mbhmpolex(:,:), mbhmpoley(:,:)
      complex *16, allocatable :: ympolexx(:,:), ympolexy(:,:),
     1     ympoleyy(:,:)
      complex *16, allocatable :: mbhmpolexx(:,:),mbhmpolexy(:,:),
     1     mbhmpoleyy(:,:)
      real *8 dc, ds, dc2, ds2, dcx, dsx, dcy, dsy
      integer ifder, ifders, l, j, i
      data eye /(0.0d0,1.0d0)/
c
      allocate(ympolex(nd,0:nterms+1),ympoley(nd,0:nterms+1),
     1     mbhmpolex(nd,0:nterms+1),mbhmpoley(nd,0:nterms+1))
      do i=0,nterms+1
         do j = 1,nd
            ympolex(j,i)=0
            ympoley(j,i)=0
            mbhmpolex(j,i)=0
            mbhmpoley(j,i)=0
         enddo
      enddo

      do j=1,nd
         ympolex(j,1) = -beta/rscale*ympole(j,0)
         ympoley(j,1) = -beta/rscale*ympole(j,0)*eye

         mbhmpolex(j,1) = -beta/rscale*mbhmpole(j,0)
         mbhmpoley(j,1) = -beta/rscale*mbhmpole(j,0)*eye
      enddo

      do i=1,nterms
         do j = 1,nd
            ympolex(j,i-1) = ympolex(j,i-1)
     1           -beta/2.0d0*ympole(j,i)*rscale
            ympolex(j,i+1) = ympolex(j,i+1)
     1           -beta/2.0d0*ympole(j,i)/rscale
            ympoley(j,i-1) = ympoley(j,i-1)
     1           +beta/2.0d0*ympole(j,i)*rscale*eye
            ympoley(j,i+1) = ympoley(j,i+1)
     1           -beta/2.0d0*ympole(j,i)/rscale*eye

            ympolex(j,i-1) = ympolex(j,i-1)
     1           -beta/2.0d0*mbhmpole(j,i)*rscale

            mbhmpolex(j,i+1) = mbhmpolex(j,i+1)
     1           -beta/2.0d0*mbhmpole(j,i)/rscale
            ympoley(j,i-1) = ympoley(j,i-1)
     1           +beta/2.0d0*mbhmpole(j,i)*rscale*eye
            mbhmpoley(j,i+1) = mbhmpoley(j,i+1)
     1           -beta/2.0d0*mbhmpole(j,i)/rscale*eye
         enddo
      enddo
c
      allocate(ympolexx(nd,0:nterms+2),ympolexy(nd,0:nterms+2),
     1     ympoleyy(nd,0:nterms+2),
     1     mbhmpolexx(nd,0:nterms+2),mbhmpolexy(nd,0:nterms+2),
     1     mbhmpoleyy(nd,0:nterms+2))

      do i=0,nterms+2
         do j=1,nd
            ympolexx(j,i)=0
            ympolexy(j,i)=0
            ympoleyy(j,i)=0
            mbhmpolexx(j,i)=0
            mbhmpolexy(j,i)=0
            mbhmpoleyy(j,i)=0
         enddo
      enddo

      do j = 1,nd
         ympolexx(j,1) = -beta/1.0d0/rscale*dreal(ympolex(j,0))
         ympolexy(j,1) = -beta/1.0d0/rscale*dreal(ympolex(j,0))*eye
         ympoleyy(j,1) = -beta/1.0d0/rscale*dreal(ympoley(j,0))*eye

         mbhmpolexx(j,1) = -beta/1.0d0/rscale*dreal(mbhmpolex(j,0))
         mbhmpolexy(j,1) =
     1        -beta/1.0d0/rscale*dreal(mbhmpolex(j,0))*eye
         mbhmpoleyy(j,1) =
     1        -beta/1.0d0/rscale*dreal(mbhmpoley(j,0))*eye
      enddo

      do i=1,nterms+1
         do j=1,nd
            ympolexx(j,i-1) = ympolexx(j,i-1)
     1  	 -beta/2.0d0*ympolex(j,i)*rscale
            ympolexx(j,i+1) = ympolexx(j,i+1)
     1           -beta/2.0d0*ympolex(j,i)/rscale
            ympolexy(j,i-1) = ympolexy(j,i-1)
     1           +beta/2.0d0*ympolex(j,i)*rscale*eye
            ympolexy(j,i+1) = ympolexy(j,i+1)
     1           -beta/2.0d0*ympolex(j,i)/rscale*eye
            ympoleyy(j,i-1) = ympoleyy(j,i-1)
     1           +beta/2.0d0*ympoley(j,i)*rscale*eye
            ympoleyy(j,i+1) = ympoleyy(j,i+1)
     1           -beta/2.0d0*ympoley(j,i)/rscale*eye

            ympolexx(j,i-1) = ympolexx(j,i-1)
     1           -beta/2.0d0*mbhmpolex(j,i)*rscale
            mbhmpolexx(j,i+1) = mbhmpolexx(j,i+1)
     1           -beta/2.0d0*mbhmpolex(j,i)/rscale
            ympolexy(j,i-1) = ympolexy(j,i-1)
     1           +beta/2.0d0*mbhmpolex(j,i)*rscale*eye
            mbhmpolexy(j,i+1) = mbhmpolexy(j,i+1)
     1           -beta/2.0d0*mbhmpolex(j,i)/rscale*eye
            ympoleyy(j,i-1) = ympoleyy(j,i-1)
     1           +beta/2.0d0*mbhmpoley(j,i)*rscale*eye
            mbhmpoleyy(j,i+1) = mbhmpoleyy(j,i+1)
     1           -beta/2.0d0*mbhmpoley(j,i)/rscale*eye
         enddo
      enddo

      allocate(mptemp1(0:nterms+6),mptemp2(0:nterms+6),
     1     diffs(0:nterms+6),kvec(0:nterms+6))

      do i = 1,ntarg
         zdiff(1) = ztarg(1,i)-center(1)
	 zdiff(2) = ztarg(2,i)-center(2)
      	 call h2cart2polar(zdiff,r,theta)

c     get values of difference functions
         ifders = 0
         call diffslogbk_fast(r,beta,rscale,diffs,ifders,ders,kvec,
     1        nterms+2)

         mptemp1(0)=diffs(0)
         mptemp2(0)=kvec(0)
         ztemp2=exp(eye*theta)
         ztemp1=ztemp2
	 do j=1,nterms+2
   	    mptemp1(j)=dcmplx(diffs(j)*dreal(ztemp1),
     1           diffs(j)*dimag(ztemp1))
	    mptemp2(j)=dcmplx(kvec(j)*dreal(ztemp1),
     1           kvec(j)*dimag(ztemp1))
            ztemp1 = ztemp1*ztemp2
	 enddo

c     evaluate
         do j = 1,nd
            pot(j,i) = pot(j,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpole(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympole(j,0))
         enddo

         do l = 1,nterms
            do j=1,nd
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpole(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpole(j,l))
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp2(l))*dreal(ympole(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympole(j,l))
            enddo
         enddo
         do j = 1,nd
            grad(j,1,i) = grad(j,1,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpolex(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympolex(j,0))
            grad(j,2,i) = grad(j,2,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpoley(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympoley(j,0))
         enddo
         do l = 1,nterms+1
            do j = 1,nd
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpolex(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpolex(j,l))
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp2(l))*dreal(ympolex(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympolex(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpoley(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpoley(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp2(l))*dreal(ympoley(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympoley(j,l))
            enddo
         enddo
         do j=1,nd
            hess(j,1,i) = hess(j,1,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpolexx(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympolexx(j,0))
            hess(j,2,i) = hess(j,2,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpolexy(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympolexy(j,0))
            hess(j,3,i) = hess(j,3,i) +
     1           dreal(mptemp1(0))*dreal(mbhmpoleyy(j,0)) +
     1           dreal(mptemp2(0))*dreal(ympoleyy(j,0))
         enddo
         do l = 1,nterms+1
            do j = 1,nd
               hess(j,1,i) = hess(j,1,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpolexx(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpolexx(j,l))
               hess(j,1,i) = hess(j,1,i) +
     1              dreal(mptemp2(l))*dreal(ympolexx(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympolexx(j,l))
               hess(j,2,i) = hess(j,2,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpolexy(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpolexy(j,l))
               hess(j,2,i) = hess(j,2,i) +
     1              dreal(mptemp2(l))*dreal(ympolexy(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympolexy(j,l))
               hess(j,3,i) = hess(j,3,i) +
     1              dreal(mptemp1(l))*dreal(mbhmpoleyy(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhmpoleyy(j,l))
               hess(j,3,i) = hess(j,3,i) +
     1              dreal(mptemp2(l))*dreal(ympoleyy(j,l)) +
     1              dimag(mptemp2(l))*dimag(ympoleyy(j,l))
            enddo
         enddo
      enddo

      return
      end

      subroutine mbh2dtaevalp_vec(nd,beta,rscale,center,mbhloc,lloc,
     1   nterms,ztarg,ntarg,
     2   pot)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     beta               : Yukawa parameter
c     rscale             : the scaling factor
c     center(2)          : expansion center
c     mbhloc(0:nterms)   : coefficients for difference-type functions
c     lloc(0:nterms)     : coefficients for Taylor series
c     nterms             : number of terms in expansions
c     ztarg(2,ntarg)     : coordinates of targets
c     ntarg              : number of targets
c
c     OUTPUT:
c
c     pot(ntarg)      : value of potential at targets
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhloc(nd,0:nterms), lloc(nd,0:nterms)
      real *8 rscale, center(2), ztarg(2,ntarg), beta
      real *8 pot(nd,ntarg), grad(nd,2,ntarg), hess(nd,3,ntarg)
      integer nterms, ntarg, nd
c     local
      real *8 zdiff(2), r, theta, ivec(0:200)
      real *8 pih, ders(0:1)
      real *8, allocatable :: diffs(:), pow(:), dpow(:)
      complex *16 z, eye, ztemp1, ztemp2
      complex *16, allocatable, dimension (:,:) :: mbhlocx,mbhlocy,
     1     mbhlocxx,mbhlocxy,mbhlocyy,llocx,llocy,llocxx,llocxy,llocyy
      complex *16, allocatable :: mptemp1(:), mptemp2(:)
      real *8 dtemp
      integer ifder, ifders, l, j, i
      data eye /(0.0d0,1.0d0)/

      allocate(pow(0:nterms+6),dpow(0:nterms+6))
      allocate(mptemp1(0:nterms+6),mptemp2(0:nterms+6))
      allocate(diffs(0:nterms+6))

c     evaluate

      do i = 1,ntarg

         zdiff(1) = ztarg(1,i)-center(1)
         zdiff(2) = ztarg(2,i)-center(2)
         call h2cart2polar(zdiff,r,theta)

c     get values of (beta*r)^k
         call mbh2d_rk(pow,dpow,r,beta,rscale,nterms+2)

c     get values of difference functions
         ifders = 0
         call diffszkik_fast(r,beta,rscale,diffs,ifders,ders,ivec,
     1        nterms+2)

         mptemp1(0)=diffs(0)
         mptemp2(0)=pow(0)
         ztemp2=exp(eye*theta)
         ztemp1=ztemp2
         do j=1,nterms+2
            mptemp1(j)=dcmplx(diffs(j)*dreal(ztemp1),
     1           diffs(j)*dimag(ztemp1))
            mptemp2(j)=dcmplx(pow(j)*dreal(ztemp1),pow(j)*dimag(ztemp1))
            ztemp1 = ztemp1*ztemp2
         enddo
         do j = 1,nd
            pot(j,i) = pot(j,i) + dreal(mptemp1(0))*dreal(mbhloc(j,0))
     1           + dreal(mptemp2(0))*dreal(lloc(j,0))
         enddo
         do l = 1,nterms
            do j = 1,nd
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp1(l))*dreal(mbhloc(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhloc(j,l))
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp2(l))*dreal(lloc(j,l)) +
     1              dimag(mptemp2(l))*dimag(lloc(j,l))
            enddo
         enddo
      enddo

      return
      end

      subroutine mbh2dtaevalg_vec(nd,beta,rscale,center,mbhloc,lloc,
     1   nterms,ztarg,ntarg,
     2   pot,grad)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     beta               : Yukawa parameter
c     rscale             : the scaling factor
c     center(2)          : expansion center
c     mbhloc(0:nterms)   : coefficients for difference-type functions
c     lloc(0:nterms)     : coefficients for Taylor series
c     nterms             : number of terms in expansions
c     ztarg(2,ntarg)     : coordinates of targets
c     ntarg              : number of targets
c
c     OUTPUT:
c
c     pot(ntarg)      : value of potential at targets
c     grad(2,ntarg)   : value of gradient at targets
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhloc(nd,0:nterms), lloc(nd,0:nterms)
      real *8 rscale, center(2), ztarg(2,ntarg), beta
      real *8 pot(nd,ntarg), grad(nd,2,ntarg), hess(nd,3,ntarg)
      integer nterms, ntarg, nd
c     local
      real *8 zdiff(2), r, theta, ivec(0:200)
      real *8 pih, ders(0:1)
      real *8, allocatable :: diffs(:), pow(:), dpow(:)
      complex *16 z, eye, ztemp1, ztemp2
      complex *16, allocatable, dimension (:,:) :: mbhlocx,mbhlocy,
     1     mbhlocxx,mbhlocxy,mbhlocyy,llocx,llocy,llocxx,llocxy,llocyy
      complex *16, allocatable :: mptemp1(:), mptemp2(:)
      real *8 dtemp
      integer ifder, ifders, l, j, i
      data eye /(0.0d0,1.0d0)/
      allocate(mbhlocx(nd,0:nterms+1),mbhlocy(nd,0:nterms+1),
     1     llocx(nd,0:nterms+1),llocy(nd,0:nterms+1))
      do i=0,nterms+1
         do j = 1,nd
            mbhlocx(j,i)=0
            mbhlocy(j,i)=0
            llocx(j,i)=0
            llocy(j,i)=0
         enddo
      enddo

      dtemp = 1.0d0/2.0d0

      do j = 1,nd
         mbhlocx(j,1) = dreal(mbhloc(j,0))*beta*rscale
         mbhlocy(j,1) = dreal(mbhloc(j,0))*beta*rscale*eye
         llocx(j,1) = dreal(mbhloc(j,0))*beta*dtemp*rscale
         llocy(j,1) = dreal(mbhloc(j,0))*beta*dtemp*rscale*eye
      enddo
      do i=1,nterms
         dtemp = dtemp/(2*(i+1))
         do j = 1,nd
            mbhlocx(j,i-1)=mbhlocx(j,i-1)+beta/2/rscale*mbhloc(j,i)
            mbhlocx(j,i+1)=mbhlocx(j,i+1)+beta/2*rscale*mbhloc(j,i)
            mbhlocy(j,i-1)=mbhlocy(j,i-1) -
     1           beta/2*(eye)/rscale*mbhloc(j,i)
            mbhlocy(j,i+1)=mbhlocy(j,i+1) +
     1           beta/2*(eye)*rscale*mbhloc(j,i)

            llocx(j,i-1)=llocx(j,i-1) + beta*i/rscale*lloc(j,i)
            llocy(j,i-1)=llocy(j,i-1) - beta*i*(eye)/rscale*lloc(j,i)

            llocx(j,i+1)=llocx(j,i+1) + beta/2*rscale*mbhloc(j,i)*dtemp
            llocy(j,i+1)=llocy(j,i+1) +
     1           beta/2*(eye)*rscale*mbhloc(j,i)*dtemp
         enddo
      enddo

      allocate(pow(0:nterms+6),dpow(0:nterms+6))
      allocate(mptemp1(0:nterms+6),mptemp2(0:nterms+6))
      allocate(diffs(0:nterms+6))

c     evaluate

      do i = 1,ntarg

         zdiff(1) = ztarg(1,i)-center(1)
         zdiff(2) = ztarg(2,i)-center(2)
         call h2cart2polar(zdiff,r,theta)

c     get values of (beta*r)^k
         call mbh2d_rk(pow,dpow,r,beta,rscale,nterms+2)

c     get values of difference functions
         ifders = 0
         call diffszkik_fast(r,beta,rscale,diffs,ifders,ders,ivec,
     1        nterms+2)

         mptemp1(0)=diffs(0)
         mptemp2(0)=pow(0)
         ztemp2=exp(eye*theta)
         ztemp1=ztemp2
         do j=1,nterms+2
            mptemp1(j)=dcmplx(diffs(j)*dreal(ztemp1),
     1           diffs(j)*dimag(ztemp1))
            mptemp2(j)=dcmplx(pow(j)*dreal(ztemp1),pow(j)*dimag(ztemp1))
            ztemp1 = ztemp1*ztemp2
         enddo
         do j = 1,nd
            pot(j,i) = pot(j,i) + dreal(mptemp1(0))*dreal(mbhloc(j,0))
     1           + dreal(mptemp2(0))*dreal(lloc(j,0))
         enddo
         do l = 1,nterms
            do j = 1,nd
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp1(l))*dreal(mbhloc(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhloc(j,l))
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp2(l))*dreal(lloc(j,l)) +
     1              dimag(mptemp2(l))*dimag(lloc(j,l))
            enddo
         enddo
         do j=1,nd
            grad(j,1,i) = grad(j,1,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocx(j,0)) +
     1           dreal(mptemp2(0))*dreal(llocx(j,0))
            grad(j,2,i) = grad(j,2,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocy(j,0)) +
     1           dreal(mptemp2(0))*dreal(llocy(j,0))
         enddo

         do l = 1,nterms+1
            do j = 1,nd
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocx(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocx(j,l))
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp2(l))*dreal(llocx(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocx(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocy(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocy(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp2(l))*dreal(llocy(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocy(j,l))
            enddo
         enddo
      enddo

      return
      end

      subroutine mbh2dtaevalh_vec(nd,beta,rscale,center,mbhloc,lloc,
     1   nterms,ztarg,ntarg,
     2   pot,grad,hess)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     beta               : Yukawa parameter
c     rscale             : the scaling factor
c     center(2)          : expansion center
c     mbhloc(0:nterms)   : coefficients for difference-type functions
c     lloc(0:nterms)     : coefficients for Taylor series
c     nterms             : number of terms in expansions
c     ztarg(2,ntarg)     : coordinates of targets
c     ntarg              : number of targets
c
c     OUTPUT:
c
c     pot(ntarg)      : value of potential at targets
c     grad(2,ntarg)   : value of gradient at targets
c     hess(3,ntarg)   : value of Hessian at targets
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhloc(nd,0:nterms), lloc(nd,0:nterms)
      real *8 rscale, center(2), ztarg(2,ntarg), beta
      real *8 pot(nd,ntarg), grad(nd,2,ntarg), hess(nd,3,ntarg)
      integer nterms, ntarg, nd
c     local
      real *8 zdiff(2), r, theta, ivec(0:200)
      real *8 pih, ders(0:1)
      real *8, allocatable :: diffs(:), pow(:), dpow(:)
      complex *16 z, eye, ztemp1, ztemp2
      complex *16, allocatable, dimension (:,:) :: mbhlocx,mbhlocy,
     1     mbhlocxx,mbhlocxy,mbhlocyy,llocx,llocy,llocxx,llocxy,llocyy
      complex *16, allocatable :: mptemp1(:), mptemp2(:)
      real *8 dtemp
      integer ifder, ifders, l, j, i
      data eye /(0.0d0,1.0d0)/
      allocate(mbhlocx(nd,0:nterms+1),mbhlocy(nd,0:nterms+1),
     1     llocx(nd,0:nterms+1),llocy(nd,0:nterms+1))
      do i=0,nterms+1
         do j = 1,nd
            mbhlocx(j,i)=0
            mbhlocy(j,i)=0
            llocx(j,i)=0
            llocy(j,i)=0
         enddo
      enddo

      dtemp = 1.0d0/2.0d0

      do j = 1,nd
         mbhlocx(j,1) = dreal(mbhloc(j,0))*beta*rscale
         mbhlocy(j,1) = dreal(mbhloc(j,0))*beta*rscale*eye
         llocx(j,1) = dreal(mbhloc(j,0))*beta*dtemp*rscale
         llocy(j,1) = dreal(mbhloc(j,0))*beta*dtemp*rscale*eye
      enddo
      do i=1,nterms
         dtemp = dtemp/(2*(i+1))
         do j = 1,nd
            mbhlocx(j,i-1)=mbhlocx(j,i-1)+beta/2/rscale*mbhloc(j,i)
            mbhlocx(j,i+1)=mbhlocx(j,i+1)+beta/2*rscale*mbhloc(j,i)
            mbhlocy(j,i-1)=mbhlocy(j,i-1) -
     1           beta/2*(eye)/rscale*mbhloc(j,i)
            mbhlocy(j,i+1)=mbhlocy(j,i+1) +
     1           beta/2*(eye)*rscale*mbhloc(j,i)

            llocx(j,i-1)=llocx(j,i-1) + beta*i/rscale*lloc(j,i)
            llocy(j,i-1)=llocy(j,i-1) - beta*i*(eye)/rscale*lloc(j,i)

            llocx(j,i+1)=llocx(j,i+1) + beta/2*rscale*mbhloc(j,i)*dtemp
            llocy(j,i+1)=llocy(j,i+1) +
     1           beta/2*(eye)*rscale*mbhloc(j,i)*dtemp
         enddo
      enddo
c
      allocate(mbhlocxx(nd,0:nterms+2),mbhlocxy(nd,0:nterms+2),
     1     mbhlocyy(nd,0:nterms+2),
     1     llocxx(nd,0:nterms+2),llocxy(nd,0:nterms+2),
     1     llocyy(nd,0:nterms+2))

      do i=0,nterms+2
         do j = 1,nd
            mbhlocxx(j,i)=0
            mbhlocxy(j,i)=0
            mbhlocyy(j,i)=0
            llocxx(j,i)=0
            llocxy(j,i)=0
            llocyy(j,i)=0
         enddo
      enddo

      dtemp = 1.0d0/2.0d0

      do j = 1,nd
         mbhlocxx(j,1) = dreal(mbhlocx(j,0))*beta*rscale
         mbhlocxy(j,1) = dreal(mbhlocx(j,0))*beta*rscale*eye
         mbhlocyy(j,1) = dreal(mbhlocy(j,0))*beta*rscale*eye
         llocxx(j,1) = dreal(mbhlocx(j,0))*beta*dtemp*rscale
         llocxy(j,1) = dreal(mbhlocx(j,0))*beta*dtemp*rscale*eye
         llocyy(j,1) = dreal(mbhlocy(j,0))*beta*dtemp*rscale*eye
      enddo

      do i=1,nterms+1
         dtemp = dtemp/(2*(i+1))
         do j = 1,nd
            mbhlocxx(j,i-1)=mbhlocxx(j,i-1)+beta/2/rscale*mbhlocx(j,i)
            mbhlocxx(j,i+1)=mbhlocxx(j,i+1)+beta/2*rscale*mbhlocx(j,i)
            mbhlocxy(j,i-1)=mbhlocxy(j,i-1)-
     1           beta/2*(eye)/rscale*mbhlocx(j,i)
            mbhlocxy(j,i+1)=mbhlocxy(j,i+1)+
     1           beta/2*(eye)*rscale*mbhlocx(j,i)
            mbhlocyy(j,i-1)=mbhlocyy(j,i-1)-
     1           beta/2*(eye)/rscale*mbhlocy(j,i)
            mbhlocyy(j,i+1)=mbhlocyy(j,i+1)+
     1           beta/2*(eye)*rscale*mbhlocy(j,i)

            llocxx(j,i-1)=llocxx(j,i-1)+beta*i/rscale*llocx(j,i)
            llocxy(j,i-1)=llocxy(j,i-1)-beta*i*(eye)/rscale*llocx(j,i)
            llocyy(j,i-1)=llocyy(j,i-1)-beta*i*(eye)/rscale*llocy(j,i)

            llocxx(j,i+1)=llocxx(j,i+1)+beta/2*rscale*mbhlocx(j,i)*dtemp
            llocxy(j,i+1)=llocxy(j,i+1)+
     1           beta/2*(eye)*rscale*mbhlocx(j,i)*dtemp
            llocyy(j,i+1)=llocyy(j,i+1)+
     1           beta/2*(eye)*rscale*mbhlocy(j,i)*dtemp
         enddo
      enddo

      allocate(pow(0:nterms+6),dpow(0:nterms+6))
      allocate(mptemp1(0:nterms+6),mptemp2(0:nterms+6))
      allocate(diffs(0:nterms+6))

c     evaluate

      do i = 1,ntarg

         zdiff(1) = ztarg(1,i)-center(1)
         zdiff(2) = ztarg(2,i)-center(2)
         call h2cart2polar(zdiff,r,theta)

c     get values of (beta*r)^k
         call mbh2d_rk(pow,dpow,r,beta,rscale,nterms+2)

c     get values of difference functions
         ifders = 0
         call diffszkik_fast(r,beta,rscale,diffs,ifders,ders,ivec,
     1        nterms+2)

         mptemp1(0)=diffs(0)
         mptemp2(0)=pow(0)
         ztemp2=exp(eye*theta)
         ztemp1=ztemp2
         do j=1,nterms+2
            mptemp1(j)=dcmplx(diffs(j)*dreal(ztemp1),
     1           diffs(j)*dimag(ztemp1))
            mptemp2(j)=dcmplx(pow(j)*dreal(ztemp1),pow(j)*dimag(ztemp1))
            ztemp1 = ztemp1*ztemp2
         enddo
         do j = 1,nd
            pot(j,i) = pot(j,i) + dreal(mptemp1(0))*dreal(mbhloc(j,0))
     1           + dreal(mptemp2(0))*dreal(lloc(j,0))
         enddo
         do l = 1,nterms
            do j = 1,nd
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp1(l))*dreal(mbhloc(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhloc(j,l))
               pot(j,i) = pot(j,i) +
     1              dreal(mptemp2(l))*dreal(lloc(j,l)) +
     1              dimag(mptemp2(l))*dimag(lloc(j,l))
            enddo
         enddo
         do j=1,nd
            grad(j,1,i) = grad(j,1,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocx(j,0)) +
     1           dreal(mptemp2(0))*dreal(llocx(j,0))
            grad(j,2,i) = grad(j,2,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocy(j,0)) +
     1           dreal(mptemp2(0))*dreal(llocy(j,0))
         enddo

         do l = 1,nterms+1
            do j = 1,nd
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocx(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocx(j,l))
               grad(j,1,i) = grad(j,1,i) +
     1              dreal(mptemp2(l))*dreal(llocx(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocx(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocy(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocy(j,l))
               grad(j,2,i) = grad(j,2,i) +
     1              dreal(mptemp2(l))*dreal(llocy(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocy(j,l))
            enddo
         enddo
         do j = 1,nd
            hess(j,1,i) = hess(j,1,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocxx(j,0))
     1           + dreal(mptemp2(0))*dreal(llocxx(j,0))
            hess(j,2,i) = hess(j,2,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocxy(j,0))
     1           + dreal(mptemp2(0))*dreal(llocxy(j,0))
            hess(j,3,i) = hess(j,3,i) +
     1           dreal(mptemp1(0))*dreal(mbhlocyy(j,0))
     1           + dreal(mptemp2(0))*dreal(llocyy(j,0))
         enddo
         do l = 1,nterms+1
            do j = 1,nd
               hess(j,1,i) = hess(j,1,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocxx(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocxx(j,l))
               hess(j,1,i) = hess(j,1,i) +
     1              dreal(mptemp2(l))*dreal(llocxx(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocxx(j,l))
               hess(j,2,i) = hess(j,2,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocxy(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocxy(j,l))
               hess(j,2,i) = hess(j,2,i) +
     1              dreal(mptemp2(l))*dreal(llocxy(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocxy(j,l))
               hess(j,3,i) = hess(j,3,i) +
     1              dreal(mptemp1(l))*dreal(mbhlocyy(j,l)) +
     1              dimag(mptemp1(l))*dimag(mbhlocyy(j,l))
               hess(j,3,i) = hess(j,3,i) +
     1              dreal(mptemp2(l))*dreal(llocyy(j,l)) +
     1              dimag(mptemp2(l))*dimag(llocyy(j,l))
            enddo
         enddo
      enddo

      return
      end

      subroutine mbh2dformmpmp_vec(nd,beta,rscale,source,ns,
     1     mbhmpolesrc,ympolesrc,ntermsrc,center,nterms,mbhmpole,ympole)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     form a multipole expansion based on multipolar sources
c
c     nd - integer, number of expansions per source
c     beta - real *8, modified biharmonic parameter
c     rscale - real *8, scaling factor for output expansions
c     source - real *8 (2,ns) source locations
c     ns - integer, number of sources
c     mbhmpolesrc - complex *16 (nd,0:ntermsrc,ns) unscaled multipolar
c                   charge of difference kernel type at each source
c     ympolesrc - complex *16 (nd,0:ntermsrc,ns) unscaled multipolar
c                   charge of Yukawa type at each source
c     center - real *8 (2), center about which to form the expansion
c     nterms - integer, number of terms in output expansion
c
c     output:
c
c     mbhmpole - complex *16 (nd,0:nterms), difference kernel type
c                expansion for source charges
c     ympole - complex *16 (nd,0:nterms), Yukawa type
c                expansion for source charges
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global variables
      integer ns, ntermsrc, nterms, nd
      real *8 beta, rscale, source(2,*)
      real *8 center(2)
      complex *16 mbhmpolesrc(nd,0:ntermsrc,ns),
     1     ympolesrc(nd,0:ntermsrc,ns)
      complex *16 mbhmpole(nd,0:nterms), ympole(nd,0:nterms)
c     local variables
      real *8 rscale1,rsc
      integer i,j,l
c
      rscale1=rscale

      do i = 1,ns

         rsc=1
         do j = 0,ntermsrc
            do l = 1,nd
               mbhmpolesrc(l,j,i)=mbhmpolesrc(l,j,i)/rsc
               ympolesrc(l,j,i)=ympolesrc(l,j,i)/rsc
            enddo
            rsc=rsc*rscale1
         enddo
         call mbh2dmpmp_vec(nd,beta,rscale1,source(1,i),
     1        mbhmpolesrc(1,0,i),ympolesrc(1,0,i),ntermsrc,
     2        rscale,center,mbhmpole,ympole,nterms)

         rsc=1
         do j = 0,ntermsrc
            do l = 1,nd
               mbhmpolesrc(l,j,i)=mbhmpolesrc(l,j,i)/rsc
               ympolesrc(l,j,i)=ympolesrc(l,j,i)/rsc
            enddo
            rsc=rsc/rscale1
         enddo

      enddo

      return
      end

      subroutine mbh2dformtamp_vec(nd,beta,rscale,source,ns,
     1     mbhmpolesrc,ympolesrc,ntermsrc,center,nterms,
     2     mbhloc,lloc)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     form a local expansion based on multipolar sources
c
c     nd - integer, number of expansions per source
c     beta - real *8, modified biharmonic parameter
c     rscale - real *8, scaling factor for output expansions
c     source - real *8 (2,ns) source locations
c     ns - integer, number of sources
c     mbhmpolesrc - complex *16 (nd,0:ntermsrc,ns) unscaled multipolar
c                   charge of difference kernel type at each source
c     ympolesrc - complex *16 (nd,0:ntermsrc,ns) unscaled multipolar
c                   charge of Yukawa type at each source
c     center - real *8 (2), center about which to form the expansion
c     nterms - integer, number of terms in output expansion
c
c     output:
c
c     mbhloc - complex *16 (nd,0:nterms), difference kernel type
c                expansion for source charges
c     lloc - complex *16 (nd,0:nterms), power series type
c                expansion for source charges
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global variables
      integer ns, ntermsrc, nterms, nd
      real *8 beta, rscale, source(2,*)
      real *8 center(2)
      complex *16 mbhmpolesrc(nd,0:ntermsrc,ns),
     1     ympolesrc(nd,0:ntermsrc,ns)
      complex *16 mbhloc(nd,0:nterms), lloc(nd,0:nterms)
c     local variables
      integer i, j, l
      real *8 rscale1,rsc

      rscale1=rscale

      do i = 1,ns

         rsc=1
         do j = 0,ntermsrc
            do l = 1,nd
               mbhmpolesrc(l,j,i)=mbhmpolesrc(l,j,i)/rsc
               ympolesrc(l,j,i)=ympolesrc(l,j,i)/rsc
            enddo
            rsc=rsc*rscale1
         enddo

         call mbh2dmploc_vec(nd,beta,rscale1,source(1,i),
     1        mbhmpolesrc(1,0,i),ympolesrc(1,0,i),ntermsrc,
     2        rscale,center,mbhloc,lloc,nterms)

         rsc=1
         do j = 0,ntermsrc
            do l = 1,nd
               mbhmpolesrc(l,j,i)=mbhmpolesrc(l,j,i)/rsc
               ympolesrc(l,j,i)=ympolesrc(l,j,i)/rsc
            enddo
            rsc=rsc/rscale1
         enddo

      enddo

      return
      end

      subroutine mbh2dmpmp_vec(nd,beta,rscale1,center1,
     1     mbhmpole1,ympole1,nterms1,rscale2,center2,mbhmpole2,ympole2,
     2     nterms2)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     translate multipole expansions
c
c     INPUT:
c
c     nd                : number of densities
c     rscale1           : scaling for original expansion
c     center1           : center for original expansion
c     mbhmpole1(nd,0:nterms1): original coeffs, difference-type functions
c     ympole1(nd,0:nterms1): original coeffs, modified Bessel functions
c     nterms1           : number of terms in original expansion
c     rscale2           : scaling for new expansion
c     center2           : center for new expansion
c     nterms2           : number of terms in original expansion
c     beta              : the modified biharmonic parameter
c
c     OUTPUT:
c
c     mbhmpole2(nd,0:nterms2): new coeffs, difference-type functions
c     ympole2(nd,0:nterms2): new coeffs, modified Bessel functions
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      real *8 rscale1, center1(2), beta
      real *8 rscale2, center2(2)
      complex *16 mbhmpole1(nd,0:nterms1), ympole1(nd,0:nterms1)
      complex *16 mbhmpole2(nd,0:nterms2), ympole2(nd,0:nterms2)
      integer nterms1, nterms2, nd
c     local
      integer i, j, ifder, nterms, m, l
      real *8 zdiff(2), r, theta, pi, ders(0:1)
      real *8, allocatable :: diffs(:), ival(:), dfac(:), dfac2(:)
      real *8, allocatable :: pow(:), dpow(:)
      real *8 twojm1
      real *8 rsj, rsi5, rsi7, rsi, fs2, rtemp
      complex *16, allocatable :: jtemp(:), difftemp(:), powtemp(:)
      complex *16 zk, z, ima, ztemp1, zmul, ztemp
      data ima /(0.0d0,1.0d0)/

      pi = 4.0d0*datan(1.0d0)
      zk = ima*beta

      nterms = nterms1+nterms2

      allocate(jtemp(0:nterms+6),difftemp(0:nterms+6),
     1     powtemp(0:nterms+6),
     2     diffs(0:nterms+6),ival(0:nterms+6),pow(0:nterms+6),
     3     dpow(0:nterms+6))

c     get local difference and BesselI vals

      zdiff(1)=center2(1)-center1(1)
      zdiff(2)=center2(2)-center1(2)
      call h2cart2polar(zdiff,r,theta)
      theta=theta-pi
      z=zk*r
      ifder=0

      call diffszkik_fast(r,beta,rscale1,diffs,ifder,ders,ival,
     1     nterms)

      rtemp = r
      call mbh2d_rksc(pow,dpow,rtemp,beta,rscale1,nterms)

c
      jtemp(0) = ival(0)
      difftemp(0) = diffs(0)
      powtemp(0) = pow(0)
      zmul=exp(ima*theta)
      ztemp1= zmul

      do j = 1,nterms
         jtemp( j) = ztemp1*ival(j)
         difftemp(j) = ztemp1*diffs(j)
         powtemp(j) = ztemp1*pow(j)
         ztemp1= ztemp1*zmul
      enddo

c     shift the expansion

      rsj = rscale1

      do l=1,nd
         ympole2(l,0) = ympole2(l,0) + dreal(ympole1(l,0))*jtemp(0)
         ympole2(l,0) = ympole2(l,0) + dreal(mbhmpole1(l,0))*difftemp(0)
         mbhmpole2(l,0) = mbhmpole2(l,0) + dreal(mbhmpole1(l,0))
      enddo

      do j = 1,nterms1
         do l = 1,nd
            ympole2(l,0) = ympole2(l,0)
     1           +(dreal(ympole1(l,j))*dreal(jtemp(j))
     1           +dimag(ympole1(l,j))*dimag(jtemp(j)))*rsj**2
            ympole2(l,0) = ympole2(l,0)
     1           +(dreal(mbhmpole1(l,j))*dreal(jtemp(j))
     1           +dimag(mbhmpole1(l,j))*dimag(jtemp(j)))*rsj**2
         enddo
         rsj=rsj*rscale1
      enddo

c
      rsi=rscale1

      rsi5=rscale1/rscale2
      do i = 1,nterms2
         do l = 1,nd
            ympole2(l,i) = ympole2(l,i)
     1           + 2.0d0*(dreal(ympole1(l,0))*dreal(jtemp(i))
     2           +ima*dreal(ympole1(l,0))*dimag(jtemp(i)))*rsi5
            ympole2(l,i) = ympole2(l,i)
     1           + 2.0d0*(dreal(mbhmpole1(l,0))*dreal(difftemp(i))
     2           +ima*dreal(mbhmpole1(l,0))*dimag(difftemp(i)))*rsi5
            mbhmpole2(l,i) = mbhmpole2(l,i)
     1           + 2.0d0*(dreal(mbhmpole1(l,0))*dreal(powtemp(i))
     2           +ima*dreal(mbhmpole1(l,0))*dimag(powtemp(i)))*rsi5
         enddo
         rsj=rscale1
         twojm1 = 1.0d0
         do j = 1,min(nterms1,i)
            do l = 1,nd
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(ympole1(l,j))*dreal(jtemp(i-j))
     1              -dimag(ympole1(l,j))*dimag(jtemp(i-j))
     2              +ima*(dreal(ympole1(l,j))*dimag(jtemp(i-j))
     3              +dimag(ympole1(l,j))*dreal(jtemp(i-j))))*rsi5
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(ympole1(l,j))*dreal(jtemp(i+j))
     1              +dimag(ympole1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(ympole1(l,j))*dimag(jtemp(i+j))
     3              -dimag(ympole1(l,j))*dreal(jtemp(i+j))))*rsj**2*rsi5
               mbhmpole2(l,i) = mbhmpole2(l,i)
     1              +(dreal(mbhmpole1(l,j))*dreal(powtemp(i-j))
     1              -dimag(mbhmpole1(l,j))*dimag(powtemp(i-j))
     2              +ima*(dreal(mbhmpole1(l,j))*dimag(powtemp(i-j))
     3              +dimag(mbhmpole1(l,j))*dreal(powtemp(i-j))))
     5              *rsi5
               ztemp = difftemp(i-j)
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(mbhmpole1(l,j))*dreal(ztemp)
     1              -dimag(mbhmpole1(l,j))*dimag(ztemp)
     2              +ima*(dreal(mbhmpole1(l,j))*dimag(ztemp)
     3              +dimag(mbhmpole1(l,j))*dreal(ztemp)))*rsi5
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(mbhmpole1(l,j))*dreal(jtemp(i+j))
     1              +dimag(mbhmpole1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(mbhmpole1(l,j))*dimag(jtemp(i+j))
     3              -dimag(mbhmpole1(l,j))*dreal(jtemp(i+j))))
     5              *rsj**2*rsi5
            enddo
            rsj=rsj*rscale1
            twojm1 = twojm1*2.0d0
         enddo
         rsj=rscale1**(i+1)
         fs2=rsi5*rscale1**2
         do j = i+1,nterms1
            do l=1,nd
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(mbhmpole1(l,j))*dreal(jtemp(j-i))
     1              +dimag(mbhmpole1(l,j))*dimag(jtemp(j-i))
     2              -ima*(dreal(mbhmpole1(l,j))*dimag(jtemp(j-i))
     3              -dimag(mbhmpole1(l,j))*dreal(jtemp(j-i))))*fs2
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(mbhmpole1(l,j))*dreal(jtemp(i+j))
     1              +dimag(mbhmpole1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(mbhmpole1(l,j))*dimag(jtemp(i+j))
     3              -dimag(mbhmpole1(l,j))*dreal(jtemp(i+j))))
     5              *rsj**2*rsi5
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(ympole1(l,j))*dreal(jtemp(j-i))
     1              +dimag(ympole1(l,j))*dimag(jtemp(j-i))
     2              -ima*(dreal(ympole1(l,j))*dimag(jtemp(j-i))
     3              -dimag(ympole1(l,j))*dreal(jtemp(j-i))))*fs2
               ympole2(l,i) = ympole2(l,i)
     1              +(dreal(ympole1(l,j))*dreal(jtemp(i+j))
     1              +dimag(ympole1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(ympole1(l,j))*dimag(jtemp(i+j))
     3              -dimag(ympole1(l,j))*dreal(jtemp(i+j))))
     5              *rsj**2*rsi5
            enddo
            rsj=rsj*rscale1
            fs2=fs2*rscale1**2
         enddo
         rsi=rsi*rscale1
         rsi5=rsi5*rscale1/rscale2
      enddo

      return
      end

      subroutine mbh2dmploc_vec(nd,beta,rscale1,center1,mbhmpole,ympole,
     1     nterms1,rscale2,center2,mbhloc,lloc,nterms2)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     rscale1           : scaling for original expansion
c     center1           : center for original expansion
c     mbhmpole(nd,0:nterms1): original coeffs, difference-type functions
c     ympole(nd,0:nterms1): original coeffs, modified Bessel functions
c     nterms1           : number of terms in original expansion
c     rscale2           : scaling for new expansion
c     center2           : center for new expansion
c     nterms2           : number of terms in original expansion
c     zcirc2(2,nterms2) : precomp'd evenly-spaced points on unit circle
c     beta              : the modified biharmonic parameter

c     work(*)           : work array. recommended length 16*nterms+500
c
c     OUTPUT:
c
c     mbhloc(nd,0:nterms2): new coeffs, difference-type functions
c     lloc(nd,0:nterms2): new coeffs, power series type
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      complex *16 mbhmpole(nd,0:nterms1), ympole(nd,0:nterms1)
      complex *16 mbhloc(nd,0:nterms2), lloc(nd,0:nterms2)
      real *8 rscale1, rscale2, center1(2), center2(2), beta
      integer nterms1, nterms2, nd
c     local
      real *8 zdiff(2), r, theta, pi, rsi, rsj, rsi5, rsi7, rsj2
      real *8 rsi52, rsi2, ders(0:1)
      real *8, allocatable :: diffs(:), kvec(:), pow(:), dpow(:),
     1     dfac2(:)
      complex *16, allocatable :: ktemp(:), difftemp(:), powtemp(:)
      complex *16 ztemp1, zmul, ima
      integer nterms, i, j, l, isign, ifders
      data ima /(0.0d0,1.0d0)/

      pi = 4.0d0*datan(1.0d0)

      nterms = nterms1+nterms2

      allocate(diffs(0:nterms+6),kvec(0:nterms+6),pow(0:nterms+6),
     1     dpow(0:nterms+6),dfac2(0:nterms+6),ktemp(0:nterms+6),
     2     difftemp(0:nterms+6),powtemp(0:nterms+6))

      dfac2(0) = 1.0d0
      do i = 1,nterms
         dfac2(i) = dfac2(i-1)/(2.0d0*i)
      enddo

      zdiff(1) = center2(1)-center1(1)
      zdiff(2) = center2(2)-center1(2)
      call h2cart2polar(zdiff,r,theta)

      theta = theta-pi

c     get values of (beta*r)^-k
      call mbh2d_rmk(pow,dpow,r,beta,rscale1,nterms)

c     get values of difference functions
      ifders = 0
      call diffslogbk_fast(r,beta,rscale1,diffs,ifders,ders,kvec,
     1     nterms)

c     form terms for shifting

      ktemp(0) = kvec(0)
      difftemp(0) = diffs(0)
      powtemp(0) = pow(0)
      zmul=-exp(ima*theta)
      ztemp1= zmul
      do j = 1,nterms
         ktemp( j) = ztemp1*kvec(j)
         difftemp(j) = ztemp1*diffs(j)
         powtemp(j) = ztemp1*pow(j)
         ztemp1= ztemp1*zmul
      enddo

c     shift

      do l = 1,nd
         mbhloc(l,0) = mbhloc(l,0) + mbhmpole(l,0)*ktemp(0)
         mbhloc(l,0) = mbhloc(l,0) + ympole(l,0)*ktemp(0)
         lloc(l,0) = lloc(l,0) + mbhmpole(l,0)*(ktemp(0)+pow(0))
         lloc(l,0) = lloc(l,0) + ympole(l,0)*ktemp(0)
      enddo
      do j = 1,nterms1
         do l = 1,nd
            mbhloc(l,0) = mbhloc(l,0) +
     1           (dreal(mbhmpole(l,j))*dreal(ktemp(j))
     1           +dimag(mbhmpole(l,j))*dimag(ktemp(j)))
            mbhloc(l,0) = mbhloc(l,0) +
     1           (dreal(ympole(l,j))*dreal(ktemp(j))
     1           +dimag(ympole(l,j))*dimag(ktemp(j)))
            lloc(l,0) = lloc(l,0) +
     1           (dreal(mbhmpole(l,j))*dreal(difftemp(j))
     1           +dimag(mbhmpole(l,j))*dimag(difftemp(j)))
            lloc(l,0) = lloc(l,0) +
     1           (dreal(ympole(l,j))*dreal(ktemp(j))
     1           +dimag(ympole(l,j))*dimag(ktemp(j)))
         enddo
      enddo
c
      rsi=rscale1
      rsi2=rscale1**2
      rsi5=rscale2/rscale1
      rsi52=rsi5*rscale2/rscale1
      isign = -1
      do i = 1,nterms2
         do l = 1,nd
            mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1           2.0d0*(dreal(mbhmpole(l,0))*dreal(ktemp(i)) +
     2           ima*dreal(mbhmpole(l,0))*dimag(ktemp(i)))
            mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1           2.0d0*(dreal(ympole(l,0))*dreal(ktemp(i))
     2           +ima*dreal(ympole(l,0))*dimag(ktemp(i)))
            lloc(l,i) = lloc(l,i) + isign*rsi5*
     1           2.0d0*(dreal(mbhmpole(l,0))*dreal(difftemp(i))
     2           + ima*dreal(mbhmpole(l,0))*dimag(difftemp(i)))*dfac2(i)
            lloc(l,i) = lloc(l,i) + isign*rsi5*
     1           2.0d0*(dreal(ympole(l,0))*dreal(ktemp(i))
     2           +ima*dreal(ympole(l,0))*dimag(ktemp(i)))*dfac2(i)
         enddo
         rsj=rscale1
         rsj2=rscale1**2
         do j = 1,min(nterms1,i)
            do l = 1,nd
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(ktemp(i-j))
     1              -dimag(mbhmpole(l,j))*dimag(ktemp(i-j))
     2              +ima*(dreal(mbhmpole(l,j))*dimag(ktemp(i-j))
     3              +dimag(mbhmpole(l,j))*dreal(ktemp(i-j))))*rsj2
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(ktemp(i+j))
     1              +dimag(mbhmpole(l,j))*dimag(ktemp(i+j))
     2              +ima*(dreal(mbhmpole(l,j))*dimag(ktemp(i+j))
     3              -dimag(mbhmpole(l,j))*dreal(ktemp(i+j))))
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(i-j))
     1              -dimag(ympole(l,j))*dimag(ktemp(i-j))
     2              +ima*(dreal(ympole(l,j))*dimag(ktemp(i-j))
     3              +dimag(ympole(l,j))*dreal(ktemp(i-j))))*rsj2
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(i+j))
     1              +dimag(ympole(l,j))*dimag(ktemp(i+j))
     2              +ima*(dreal(ympole(l,j))*dimag(ktemp(i+j))
     3              -dimag(ympole(l,j))*dreal(ktemp(i+j))))
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(ktemp(i-j))
     1              -dimag(mbhmpole(l,j))*dimag(ktemp(i-j))
     2              +ima*(dreal(mbhmpole(l,j))*dimag(ktemp(i-j))
     3              +dimag(mbhmpole(l,j))*dreal(ktemp(i-j))))*rsj2*
     5              dfac2(i)
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(difftemp(i+j))
     1              +dimag(mbhmpole(l,j))*dimag(difftemp(i+j))
     2              +ima*(dreal(mbhmpole(l,j))*dimag(difftemp(i+j))
     3              -dimag(mbhmpole(l,j))*dreal(difftemp(i+j))))*
     5              dfac2(i)
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(i-j))
     1              -dimag(ympole(l,j))*dimag(ktemp(i-j))
     2              +ima*(dreal(ympole(l,j))*dimag(ktemp(i-j))
     3              +dimag(ympole(l,j))*dreal(ktemp(i-j))))*rsj2*
     5              dfac2(i)
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(i+j))
     1              +dimag(ympole(l,j))*dimag(ktemp(i+j))
     2              +ima*(dreal(ympole(l,j))*dimag(ktemp(i+j))
     3              -dimag(ympole(l,j))*dreal(ktemp(i+j))))*dfac2(i)
            enddo
            rsj=rsj*rscale1
            rsj2=rsj2*rscale1**2
         enddo
         do j = i+1,nterms1
            do l = 1,nd
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(ktemp(j-i))
     1              +dimag(mbhmpole(l,j))*dimag(ktemp(j-i))
     2              -ima*(dreal(mbhmpole(l,j))*dimag(ktemp(j-i))
     3              -dimag(mbhmpole(l,j))*dreal(ktemp(j-i))))*rsi2
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(ktemp(i+j))
     1              +dimag(mbhmpole(l,j))*dimag(ktemp(i+j))
     2              +ima*(dreal(mbhmpole(l,j))*dimag(ktemp(i+j))
     3              -dimag(mbhmpole(l,j))*dreal(ktemp(i+j))))
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(j-i))
     1              +dimag(ympole(l,j))*dimag(ktemp(j-i))
     2              -ima*(dreal(ympole(l,j))*dimag(ktemp(j-i))
     3              -dimag(ympole(l,j))*dreal(ktemp(j-i))))*rsi2
               mbhloc(l,i) = mbhloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(i+j))
     1              +dimag(ympole(l,j))*dimag(ktemp(i+j))
     2              +ima*(dreal(ympole(l,j))*dimag(ktemp(i+j))
     3              -dimag(ympole(l,j))*dreal(ktemp(i+j))))
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(ktemp(j-i))
     1              +dimag(mbhmpole(l,j))*dimag(ktemp(j-i))
     2              -ima*(dreal(mbhmpole(l,j))*dimag(ktemp(j-i))
     3              -dimag(mbhmpole(l,j))*dreal(ktemp(j-i))))*rsi2*
     5              dfac2(i)
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(mbhmpole(l,j))*dreal(difftemp(i+j))
     1              +dimag(mbhmpole(l,j))*dimag(difftemp(i+j))
     2              +ima*(dreal(mbhmpole(l,j))*dimag(difftemp(i+j))
     3              -dimag(mbhmpole(l,j))*dreal(difftemp(i+j))))*
     5              dfac2(i)
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(j-i))
     1              +dimag(ympole(l,j))*dimag(ktemp(j-i))
     2              -ima*(dreal(ympole(l,j))*dimag(ktemp(j-i))
     3              -dimag(ympole(l,j))*dreal(ktemp(j-i))))*rsi2*
     5              dfac2(i)
               lloc(l,i) = lloc(l,i) + isign*rsi5*
     1              (dreal(ympole(l,j))*dreal(ktemp(i+j))
     1              +dimag(ympole(l,j))*dimag(ktemp(i+j))
     2              +ima*(dreal(ympole(l,j))*dimag(ktemp(i+j))
     3              -dimag(ympole(l,j))*dreal(ktemp(i+j))))*dfac2(i)
            enddo
         enddo
         isign = -isign
         rsi=rsi*rscale1
         rsi2=rsi2*rscale1**2
         rsi5=rsi5*rscale2/rscale1
      enddo

      return
      end

      subroutine mbh2dlocloc_vec(nd,beta,rscale1,center1,mbhloc1,lloc1,
     1     nterms1,rscale2,center2,mbhloc2,lloc2,nterms2,carray,ldc)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     rscale1           : scaling for original expansion
c     center1           : center for original expansion
c     mbhloc1(nd,0:nterms1): original coeffs, difference-type functions
c     lloc1(nd,0:nterms1): original coeffs, modified Bessel functions
c     nterms1           : number of terms in original expansion
c     rscale2           : scaling for new expansion
c     center2           : center for new expansion
c     nterms2           : number of terms in original expansion
c     beta              : the modified biharmonic parameter
c     carray(0:ldc,0:ldc) : precomputed array of binomial coefficients
c     ldc               : dimension of carray
c
c     OUTPUT:
c
c     mbhloc2(nd,0:nterms2): new coeffs, difference-type functions
c     lloc2(nd,0:nterms2): new coeffs, power series-like functions
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      real *8 rscale1, center1(2), beta, rscale2, center2(2)
      complex *16 mbhloc1(nd,0:nterms1), lloc1(nd,0:nterms1)
      complex *16 mbhloc2(nd,0:nterms2), lloc2(nd,0:nterms2)
      real *8 carray(0:ldc,0:ldc)
      integer ldc
      integer nterms1, nterms2, nd
c     local
      complex *16 z,ima, zmul,zinv,ztemp1,ztemp2
      complex *16, allocatable :: jtemp(:), difftemp(:), powtemp(:)
      real *8 zdiff(2), r, theta, pi, done
      real *8 rsj, rsi, rsi7, rsi5, fs2, ders(0:1)
      real *8, allocatable :: ival(:), diffs(:), pow(:), dpow(:),
     1     pow2(:), dfac2(:)
      integer i, j, m, l, nterms, ifder
c
      data ima/(0.0d0,1.0d0)/
c
      done=1
      pi=4*atan(done)
c
      nterms = nterms1+nterms2
c

      allocate(jtemp(0:nterms+6),difftemp(0:nterms+6),
     1     powtemp(0:nterms+6),ival(0:nterms+6),diffs(0:nterms+6),
     1     pow(0:nterms+6),dpow(0:nterms+6),pow2(0:nterms+6),
     1     dfac2(0:nterms+6))

      zdiff(1)=center2(1)-center1(1)
      zdiff(2)=center2(2)-center1(2)
      call h2cart2polar(zdiff,r,theta)
      theta=theta-pi

      ifder = 0
      call diffszkik_fast(r,beta,rscale1,diffs,ifder,ders,ival,
     1     nterms)

      call mbh2d_rk(pow,dpow,r,beta,rscale1,nterms)
      call mbh2d_rksc(pow2,dpow,r,beta,rscale1,nterms)

      dfac2(0) = 1.0d0
      do i = 1,nterms
         dfac2(i) = dfac2(i-1)/(2.0d0*i)
      enddo
c
      jtemp(0) = ival(0)
      powtemp(0) = pow(0)
      difftemp(0) = diffs(0)
      zmul=-exp(ima*theta)
      ztemp1= zmul
      do j = 1,nterms
         jtemp( j) = ztemp1*ival(j)
         difftemp(j) = ztemp1*diffs(j)
         powtemp(j) = ztemp1*pow(j)
         ztemp1= ztemp1*zmul
      enddo
c
      do l = 1,nd
         mbhloc2(l,0) = mbhloc2(l,0) + mbhloc1(l,0)*jtemp(0)
         lloc2(l,0) = lloc2(l,0) + mbhloc1(l,0)*difftemp(0)
         lloc2(l,0) = lloc2(l,0) + lloc1(l,0)*powtemp(0)
      enddo

      rsj=rscale1
      do j = 1,nterms1
         do l = 1,nd
            mbhloc2(l,0) = mbhloc2(l,0)+
     1           (dreal(mbhloc1(l,j))*dreal(jtemp(j))
     1           +dimag(mbhloc1(l,j))*dimag(jtemp(j)))
            lloc2(l,0) = lloc2(l,0)+
     1           (dreal(mbhloc1(l,j))*dreal(difftemp(j))
     1           +dimag(mbhloc1(l,j))*dimag(difftemp(j)))
            lloc2(l,0) = lloc2(l,0)+
     1           (dreal(lloc1(l,j))*dreal(powtemp(j))
     1           +dimag(lloc1(l,j))*dimag(powtemp(j)))
         enddo
      enddo
c
      rsi=rscale1
      rsi7=rscale2
      rsi5=rscale2/rscale1
      do i = 1,nterms2
         do l = 1,nd
            mbhloc2(l,i) = mbhloc2(l,i)
     1           + 2.0d0*(dreal(mbhloc1(l,0))*dreal(jtemp(i))
     2           +ima*dreal(mbhloc1(l,0))*dimag(jtemp(i)))*rsi7*rsi
            lloc2(l,i) = lloc2(l,i)
     1           + 2.0d0*(dreal(mbhloc1(l,0))*dreal(jtemp(i))
     2           +ima*dreal(mbhloc1(l,0))*dimag(jtemp(i)))*rsi7*rsi*
     1           dfac2(i)
         enddo
         fs2=rsi5*rscale1**2
         if( nterms1 .le. i-1 ) fs2=fs2*rscale1**(2*(i-1-nterms1))
         do j = min(nterms1,i-1),1,-1
            do l = 1,nd
               mbhloc2(l,i) = mbhloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(i-j))
     1              -dimag(mbhloc1(l,j))*dimag(jtemp(i-j))
     2              +ima*(dreal(mbhloc1(l,j))*dimag(jtemp(i-j))
     3              +dimag(mbhloc1(l,j))*dreal(jtemp(i-j))))*fs2
               mbhloc2(l,i) = mbhloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(i+j))
     1              +dimag(mbhloc1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(mbhloc1(l,j))*dimag(jtemp(i+j))
     3              -dimag(mbhloc1(l,j))*dreal(jtemp(i+j))))*rsi*rsi7
               lloc2(l,i) = lloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(i-j))
     1              -dimag(mbhloc1(l,j))*dimag(jtemp(i-j))
     2              +ima*(dreal(mbhloc1(l,j))*dimag(jtemp(i-j))
     3              +dimag(mbhloc1(l,j))*dreal(jtemp(i-j))))*fs2*
     5              dfac2(i)
               lloc2(l,i) = lloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(i+j))
     1              +dimag(mbhloc1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(mbhloc1(l,j))*dimag(jtemp(i+j))
     3              -dimag(mbhloc1(l,j))*dreal(jtemp(i+j))))
     4              *rsi*rsi7*dfac2(i)
            enddo
            fs2=fs2*rscale1**2
         enddo
         do j = i,nterms1
            do l = 1,nd
               mbhloc2(l,i) = mbhloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(j-i))
     1              +dimag(mbhloc1(l,j))*dimag(jtemp(j-i))
     2              -ima*(dreal(mbhloc1(l,j))*dimag(jtemp(j-i))
     3              -dimag(mbhloc1(l,j))*dreal(jtemp(j-i))))*rsi5
               lloc2(l,i) = lloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(difftemp(j-i))
     1              +dimag(mbhloc1(l,j))*dimag(difftemp(j-i))
     2              -ima*(dreal(mbhloc1(l,j))*dimag(difftemp(j-i))
     3              -dimag(mbhloc1(l,j))*dreal(difftemp(j-i))))*rsi5
     1              *dfac2(i)
               lloc2(l,i) = lloc2(l,i)+
     1              (dreal(lloc1(l,j))*dreal(powtemp(j-i))
     1              +dimag(lloc1(l,j))*dimag(powtemp(j-i))
     2              -ima*(dreal(lloc1(l,j))*dimag(powtemp(j-i))
     3              -dimag(lloc1(l,j))*dreal(powtemp(j-i))))*carray(j,i)
     1              *rsi5
               mbhloc2(l,i) = mbhloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(i+j))
     1              +dimag(mbhloc1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(mbhloc1(l,j))*dimag(jtemp(i+j))
     3              -dimag(mbhloc1(l,j))*dreal(jtemp(i+j))))*rsi*rsi7
               lloc2(l,i) = lloc2(l,i)+
     1              (dreal(mbhloc1(l,j))*dreal(jtemp(i+j))
     1              +dimag(mbhloc1(l,j))*dimag(jtemp(i+j))
     2              +ima*(dreal(mbhloc1(l,j))*dimag(jtemp(i+j))
     3              -dimag(mbhloc1(l,j))*dreal(jtemp(i+j))))
     4              *rsi*rsi7*dfac2(i)

            enddo
         enddo
         rsi=rsi*rscale1
         rsi7=rsi7*rscale2
         rsi5=rsi5*rscale2/rscale1
      enddo
      return
      end

      subroutine mbh2d_rk(pow,dpow,r,beta,rscale,nterms)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     r               : radius
c     rscale          : scaling factor
c     nterms          : number of terms in the expansion
c
c     OUTPUT:
c
c     pow(0:nterms)   : pow(i) = (r/rscale)^i
c     dpow(0:nterms)  : dpow(i) = i*(r/rscale)^(i-1)/rscale
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      real *8 pow(0:nterms), dpow(0:nterms), r, rscale, beta
      integer nterms
c     local
      real *8 dtemp1, dtemp2
      integer i

      dtemp2 = (r*beta)/rscale
      dtemp1 = 1.0d0

      pow(0) = 1.0d0
      dpow(0) = 0.0d0

      do i = 1,nterms
         dpow(i) = i*beta*dtemp1/rscale
         dtemp1 = dtemp1*dtemp2
         if (dtemp1 .gt. 1.0d250) dtemp1 = 0.0d0
         pow(i) = dtemp1
      enddo

      return
      end

      subroutine mbh2d_rksc(pow,dpow,r,beta,rscale,nterms)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     r               : radius
c     rscale          : scaling factor
c     nterms          : number of terms in the expansion
c
c     OUTPUT:
c
c     pow(0:nterms)   : pow(i) = (r*beta/(2*rscale))^i/(i!)
c     dpow(0:nterms)  : dpow(i) = i*(r/rscale)^(i-1)/rscale
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      real *8 pow(0:nterms), dpow(0:nterms), r, rscale, beta
      integer nterms
c     local
      real *8 dtemp1, dtemp2
      integer i

      dtemp2 = (r*beta)/(2.0d0*rscale)
      dtemp1 = 1.0d0

      pow(0) = 1.0d0
      dpow(0) = 0.0d0

      do i = 1,nterms
         dpow(i) = beta*dtemp1/(2.0d0*rscale)
         dtemp1 = dtemp1*dtemp2/i
         if (dtemp1 .gt. 1.0d250) dtemp1 = 0.0d0
         pow(i) = dtemp1
      enddo

      return
      end

      subroutine mbh2d_rmk(pow,dpow,r,beta,rscale,nterms)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     INPUT:
c
c     r               : radius
c     rscale          : scaling factor
c     nterms          : number of terms in the expansion
c
c     OUTPUT:
c
c     pow(0:nterms)   : pow(i) = (rscale/r)^(-i)
c     dpow(0:nterms)  : dpow(i) = -i*(rscale/r)^(-i-1)*rscale
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c     global
      real *8 pow(0:nterms), dpow(0:nterms), r, rscale, beta
      integer nterms
c     local
      real *8 dtemp1, dtemp2
      integer i

      dtemp2 = rscale/(r*beta)
      dtemp1 = dtemp2

      pow(0) = dlog(r)
      dpow(0) = 1.0d0/r

      do i = 1,nterms
         pow(i) = dtemp1
         if (dtemp1 .lt. 1.0d-200) dtemp1 = 0.0d0
         dtemp1 = dtemp1*dtemp2
         dpow(i) = -i*dtemp1*rscale/beta
      enddo

      return
      end

c
c

      subroutine mbh2d_init_carray(carray,ldc)
      implicit real *8 (a-h,o-z)
      real *8 carray(0:ldc,0:ldc)

      do l = 0,ldc
         carray(l,0) = 1.0d0
      enddo
      do m=1,ldc
         carray(m,m) = 1.0d0
         do l=m+1,ldc
            carray(l,m)=carray(l-1,m)+carray(l-1,m-1)
         enddo
      enddo
c
      return
      end


text:
C    Copyright (c) 2014, Sandia Corporation.
C    Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
C    the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C        * Redistributions of source code must retain the above copyright
C          notice, this list of conditions and the following disclaimer.
C
C        * Redistributions in binary form must reproduce the above
C          copyright notice, this list of conditions and the following
C          disclaimer in the documentation and/or other materials provided
C          with the distribution.
C
C        * Neither the name of Sandia Corporation nor the names of its
C          contributors may be used to endorse or promote products derived
C          from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C $Id: chkrgn.f,v 1.4 2004/01/21 05:18:39 gdsjaar Exp $
C $Log: chkrgn.f,v $
C Revision 1.4  2004/01/21 05:18:39  gdsjaar
C Initialized several variables identified by valgrind.
C
C Revision 1.3  2001/11/05 13:26:51  gdsjaar
C  Fixed array boundary problem in region check code.
C
C Revision 1.2  1990/11/30 11:25:08  gdsjaar
C Added initialization and MDSTAT calls
C
c Revision 1.1.1.1  90/11/30  11:04:38  gdsjaar
c FASTQ Version 2.0X
c
c Revision 1.1  90/11/30  11:04:37  gdsjaar
c Initial revision
c
C
CC* FILE: [.QMESH]CHKRGN
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/31/90
CC* MODIFICATION: ADDED ARGUMENTS TO CALL TO CHKHOL TO PASS MINIMUM
CC**              ELEMENT SIZE (SIZMIN) AND GETSIZ PARAMETERS OF
CC**              EMIN AND EMAX
C
      SUBROUTINE CHKRGN (IA, L, MP, ML, MS, MR, MSC, N24, IPOINT, COOR,
     &   IPBOUN, ILINE, LTYPE, NINT, FACTOR, LCON, ILBOUN, ISBOUN,
     &   ISIDE, NLPS, IFLINE, ILLIST, IREGN, NSPR, IFSIDE, ISLIST,
     &   NPPF, IFPB, LISTPB, NLPF, IFLB, LISTLB, NSPF, IFSB, LISTSB,
     &   IFHOLE, NHPR, IHLIST, LINKP, LINKL, LINKS, LINKR, LINKSC,
     &   LINKPB, LINKLB, LINKSB, RSIZE, SCHEME, DEFSCH, NPREGN, NPSBC,
     &   NPNODE, MAXNP, MAXNL, MAX3, X, Y, NID, LISTL, NNPS, ANGLE,
     &   MARKED, MXND, MXNPER, MXNL, MAXNBC, MAXSBC, AMESUR, XNOLD,
     &   YNOLD, NXKOLD, MMPOLD, LINKEG, LISTEG, BMESUR, MLINK, NPROLD,
     &   NPNOLD, NPEOLD, NNXK, REMESH, REXMIN, REXMAX, REYMIN, REYMAX,
     &   IDIVIS, SIZMIN, EMAX, EMIN, NOROOM, ERRCHK, ERR)
C***********************************************************************
C
C  SUBROUTINE CHKRGN - CHECK THAT A REGION MAY BE MESHED
C
C***********************************************************************
C
      DIMENSION IA(1)
      DIMENSION IPOINT(MP), COOR(2, MP), IPBOUN(MP)
      DIMENSION ILINE(ML), LTYPE(ML), NINT(ML), FACTOR(ML), LCON(3, ML)
      DIMENSION ILBOUN(ML), ISBOUN(ML), ISIDE(MS), NLPS(MS)
      DIMENSION IFLINE(MS), ILLIST(MS*3)
      DIMENSION IREGN(MR), NSPR(MR), IFSIDE(MR), ISLIST(MR*4)
      DIMENSION SCHEME(MSC), RSIZE (MR)
      DIMENSION NPPF(MP), IFPB(MP), LISTPB(2, MP)
      DIMENSION NLPF(ML), IFLB(ML), LISTLB(2, ML)
      DIMENSION NSPF(ML), IFSB(ML), LISTSB(2, ML)
      DIMENSION LINKP(2, MP), LINKL(2, ML), LINKS(2, MS)
      DIMENSION LINKR(2, MR), LINKSC(2, MR), LINKPB(2, MP)
      DIMENSION LINKLB(2, ML), LINKSB(2, ML)
      DIMENSION X(MAXNP), Y(MAXNP), NID(MAXNP)
      DIMENSION LISTL(MAXNL), MARKED(3, MAXNL)
      DIMENSION NNPS(MAX3), ANGLE(MAXNP)
      DIMENSION IFHOLE(MR), NHPR(MR), IHLIST(MR*2)
C
      DIMENSION IDUMMY(1)
C
      DIMENSION AMESUR(NPEOLD), XNOLD(NPNOLD), YNOLD(NPNOLD)
      DIMENSION NXKOLD(NNXK, NPEOLD), MMPOLD(3, NPROLD)
      DIMENSION LINKEG(2, MLINK), LISTEG(4 * NPEOLD), BMESUR(NPNOLD)
C
      CHARACTER*72 SCHEME, DEFSCH, SCHSTR
C
      LOGICAL NOROOM, EVEN, ERR, NORM, CCW, REAL, ADDLNK, REMESH
      LOGICAL PENTAG, TRIANG, TRNSIT, HALFC, COUNT, ERRCHK
C
      ipntr = 0
      addlnk = .false.
      COUNT = .TRUE.
      IF (REMESH) THEN
         EVEN = .TRUE.
      ELSE
         EVEN = .FALSE.
      ENDIF
      REAL = .FALSE.
C
C  CHECK TO MAKE SURE CONNECTING DATA FOR THE REGION EXISTS
C  AND FILL IN ANY BLANK INTERVALS ACCORDING TO THE GIVEN SIZE
C  FOR THE REGION AND THE LINE'S LENGTH
C
      CALL DATAOK (MP, ML, MS, MR, L, IREGN(L), COOR, ILINE, LTYPE,
     &   NINT, LCON, NLPS, IFLINE, ILLIST, NSPR, IFSIDE, ISLIST, LINKP,
     &   LINKL, LINKS, RSIZE(L), ERRCHK, ERR)
      IF (ERR) THEN
         WRITE (*, 10000) IREGN(L)
         ADDLNK = .TRUE.
         IMINUS = -L
         CALL LTSORT (MR, LINKR, IREGN(L), IMINUS, ADDLNK)
         ADDLNK = .FALSE.
C
C  CALCULATE THE PERIMETER OF THE REGION
C
      ELSE
         KNBC = 0
         KSBC = 0
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/31/90
CC* MODIFICATION: ADDED ARGUMENTS TO CALL TO PERIM TO PASS MINIMUM
CC**              ELEMENT SIZE (SIZMIN) AND GETSIZ PARAMETERS OF
CC**              EMIN AND EMAX
C
         CALL PERIM (MP, ML, MS, NSPR(L), MAXNL, MAXNP, 1, 1, KNBC,
     &      KSBC, IREGN(L), IPOINT, COOR, IPBOUN, ILINE, LTYPE, NINT,
     &      FACTOR, LCON, ILBOUN, ISBOUN, ISIDE, NLPS, IFLINE, ILLIST,
     &      ISLIST(IFSIDE(L)), NPPF, IFPB, LISTPB, NLPF, IFLB, LISTLB,
     &      NSPF, IFSB, LISTSB, LINKP, LINKL, LINKS, LINKPB, LINKLB,
     &      LINKSB, X, Y, NID, NPER, LISTL, NL, IDUMMY, MARKED, EVEN,
     &      REAL, ERR, CCW, COUNT, NOROOM, AMESUR, XNOLD, YNOLD, NXKOLD,
     &      MMPOLD, LINKEG, LISTEG, BMESUR, MLINK, NPROLD, NPNOLD,
     &      NPEOLD, NNXK, REMESH, REXMIN, REXMAX, REYMIN, REYMAX,
     &      IDIVIS, SIZMIN, EMAX, EMIN)
         IF ((NPER .LE. 0) .OR. (ERR)) THEN
            WRITE (*, 10010) IREGN(L)
            ADDLNK = .TRUE.
            IMINUS = -L
            CALL LTSORT (MR, LINKR, IREGN(L), IMINUS, ADDLNK)
            ADDLNK = .FALSE.
            GO TO 120
         END IF
C
C  WHEN CHECKING THE MAXIMUMS - ADD ENOUGH FOR ONE MORE INTERVAL
C  ON THE LINE AS THIS LINE MAY BE INCREMENTED BY ONE IF THE
C  PERIMETER IS ODD
C
         MAXNBC = MAX(MAXNBC, KNBC + 3)
         MAXSBC = MAX(MAXSBC, KSBC + 3)
         MXNL = MAX(MXNL, NL)
C
C  GET THE REGION SCHEME
C
         CALL LTSORT (MR, LINKSC, ABS(IREGN(L)), IPNTR, ADDLNK)
         IF ((IREGN(L) .LE. N24) .AND. (IPNTR .GT. 0)) THEN
            SCHSTR = SCHEME(IPNTR)
         ELSE
            SCHSTR = DEFSCH
         END IF
C
C  SEE IF A TRIANGULAR, PENTAGON, SEMICIRCLE, OR A TRANSITION
C  REGION HAS BEEN FLAGGED
C
         PENTAG = .FALSE.
         TRNSIT = .FALSE.
         TRIANG = .FALSE.
         CALL STRCUT (SCHSTR)
         CALL STRLNG (SCHSTR, LENSCH)
         DO 100 J = 1, LENSCH
            IF ((SCHSTR(J:J) .EQ. 'T') .OR.
     &         (SCHSTR(J:J) .EQ. 't')) THEN
               IF (NPER .GE. 6) THEN
                  TRIANG = .TRUE.
               ELSE
                  CALL MESAGE ('TRIANGULAR REGION MESH NOT')
                  CALL MESAGE ('POSSIBLE WITH PERIMETER < 6')
                  CALL MESAGE ('REGULAR PROCESSING ASSUMED')
               END IF
               GO TO 110
            ELSE IF ((SCHSTR(J:J) .EQ. 'U') .OR.
     &         (SCHSTR(J:J) .EQ. 'u')) THEN
               IF (NPER .GE. 10) THEN
                  PENTAG = .TRUE.
               ELSE
                  CALL MESAGE ('PENTAGON REGION MESH NOT')
                  CALL MESAGE ('POSSIBLE WITH PERIMETER < 10')
                  CALL MESAGE ('REGULAR PROCESSING ASSUMED')
               END IF
               GO TO 110
            ELSE IF ((SCHSTR(J:J) .EQ. 'B') .OR.
     &         (SCHSTR(J:J) .EQ. 'b')) THEN
               IF (NPER .GE. 8) THEN
                  TRNSIT = .TRUE.
                  HALFC = .FALSE.
               ELSE
                  CALL MESAGE ('TRANSITION REGION GENERATION NOT')
                  CALL MESAGE ('POSSIBLE WITH NO. IN PERIMETER < 8')
                  CALL MESAGE ('REGULAR PROCESSING WILL BE ATTEMPTED')
               END IF
               GO TO 110
            ELSE IF ((SCHSTR(J:J) .EQ. 'C') .OR.
     &         (SCHSTR(J:J) .EQ. 'c'))
     &         THEN
               IF (NPER .GE. 8) THEN
                  TRNSIT = .TRUE.
                  HALFC = .TRUE.
               ELSE
                  CALL MESAGE ('SEMICIRCLE REGION GENERATION NOT')
                  CALL MESAGE ('POSSIBLE WITH NO. IN PERIMETER < 8')
                  CALL MESAGE ('REGULAR PROCESSING WILL BE ATTEMPTED')
               END IF
               GO TO 110
            END IF
  100    CONTINUE
  110    CONTINUE
C
C  SET UP THE TRIANGLE DIVISIONS, AND FIND THE CENTER POINT
C
         IF (TRIANG) THEN
            CALL GETM3 (ML, MS, MAX3, NSPR(L), ISLIST(IFSIDE(L)), NINT,
     &         IFLINE, NLPS, ILLIST, LINKL, LINKS, X, Y, NID, NNPS,
     &         ANGLE, NPER, M1A, M1B, M2A, M2B, M3A, M3B, XCEN, YCEN,
     &         CCW, ERR)
C
C  CHECK FOR MAXIMUM DIMENSIONS NEEDED FOR EACH REGION
C  ASSUMING THAT 10 NECKLACES WILL BE ADEQUATE
C
            MXTEST = ((M1A + 1)*(M3B + 1)) + ((M1B + 1)*(M2A + 1))
     &         + ((M2B + 1)*(M3A + 1)) + (10*(NPER + 1)) + (NPER*2)
            MXND = MAX(MXTEST, MXND)
            MXNPER = MAX(MXNPER, (NPER + 2) * 2)
C
C  SET UP THE TRANSITION DIVISIONS, AND FIND THE CENTER POINT
C
         ELSE IF (TRNSIT) THEN
            CALL GETTRN (ML, MS, MAX3, NSPR(L), ISLIST(IFSIDE(L)), NINT,
     &         IFLINE, NLPS, ILLIST, LINKL, LINKS, X, Y, NID, NNPS,
     &         ANGLE, NPER, I1, I2, I3, I4, I5, I6, I7, I8, XCEN1,
     &         YCEN1, XCEN2, YCEN2, XMID1, YMID1, XMID2, YMID2, CCW,
     &         HALFC, ERR)
C
C  CHECK FOR MAXIMUM DIMENSIONS NEEDED FOR EACH REGION
C  ASSUMING THAT 10 NECKLACES WILL BE ADEQUATE
C
            MXTEST = ((I2 - I1)*(NPER - I8))
     &         + ((I3 - I2)*(NPER - I8)) + ((I3 - I2)*(I2 - I2))
     &         + ((I4 - I3)*(I7 - I6)) + ((I5 - I4)*(I6 - I5))
     &         + ((I5 - I4)*(I7 - I6)) + (10*(NPER + 1)) + (NPER*2)
            MXND = MAX(MXTEST, MXND)
            MXNPER = MAX(MXNPER, (NPER + 2) * 2)
C
C  SET UP THE PENTAGON DIVISIONS, AND FIND THE CENTER POINT
C
         ELSE IF (PENTAG) THEN
            CALL GETM5 (IA, ML, MS, MAX3, NSPR(L), ISLIST(IFSIDE(L)),
     &         NINT, IFLINE, NLPS, ILLIST, LINKL, LINKS, X, Y, NID,
     &         NNPS, ANGLE, NPER, M1A, M1B, M2, M3A, M3B, M4A, M4B,
     &         M5, MC, XCEN, YCEN, CCW, ERR)
C
C  CHECK FOR MAXIMUM DIMENSIONS NEEDED FOR THE REGION
C  ASSUMING THAT 10 NECKLACES WILL BE ADEQUATE
C
            MXTEST = (M1B*M2) + (M4A*M3B) + (M4B*M5) + (10*(NPER + 1))
            MXND = MAX(MXTEST, MXND)
            MXNPER = MAX(MXNPER, (NPER + 2) *2)
C
C  CALCULATE THE BASE OF THE RECTANGLE FOR THE REGION
C
         ELSE
            CALL GETM1 (ML, MS, MAX3, NSPR(L), ISLIST(IFSIDE(L)), NINT,
     &         IFLINE, NLPS, ILLIST, LINKL, LINKS, X, Y, NID, NNPS,
     &         ANGLE, NPER, SCHSTR, M1, CCW, NORM, REAL, ERR)
            IF (ERR) THEN
               WRITE (*, 10020) IREGN(L)
               ADDLNK = .TRUE.
               IMINUS = -L
               CALL LTSORT (MR, LINKR, IREGN(L), IMINUS, ADDLNK)
               ADDLNK = .FALSE.
               GO TO 120
            END IF
            M2 = NPER/2 - M1
C
C  CHECK FOR MAXIMUM DIMENSIONS NEEDED FOR EACH REGION
C  ASSUMING THAT 10 NECKLACES WILL BE ADEQUATE
C
            MXTEST = ((M1 + 1)*(M2 + 1)) + (10*(M1 + M2 + 2))
            MXND = MAX(MXTEST, MXND)
            MXNPER = MAX(MXNPER, NPER + 4)
         END IF
C
C  FLAG THE REGION AS BEING PROCESSABLE
C
         IREGN(L) = -IREGN(L)
C
C  MARK THE LINES AND POINTS IN THE REGION AS BEING USED
C
         CALL MKUSED (MAXNL, MP, ML, LISTL, IPOINT, NINT, LINKP, LINKL,
     &      LCON, NL)
C
C  CHECK ALL THE HOLES IN THE REGION
C
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/31/90
CC* MODIFICATION: ADDED ARGUMENTS TO CALL TO CHKHOL TO PASS MINIMUM
CC**              ELEMENT SIZE (SIZMIN) AND GETSIZ PARAMETERS OF
CC**              EMIN AND EMAX
C
         CALL CHKHOL (IA, L, MP, ML, MS, MR, MSC, IPOINT, COOR,
     &      IPBOUN, ILINE, LTYPE, NINT, FACTOR, LCON, ILBOUN, ISBOUN,
     &      ISIDE, NLPS, IFLINE, ILLIST, IREGN, NSPR, IFSIDE, ISLIST,
     &      NPPF, IFPB, LISTPB, NLPF, IFLB, LISTLB, NSPF, IFSB, LISTSB,
     &      IFHOLE, NHPR, IHLIST, LINKP, LINKL, LINKS, LINKR, LINKSC,
     &      LINKPB, LINKLB, LINKSB, RSIZE, NPREGN, NPSBC, NPNODE,
     &      MAXNP, MAXNL, MXNPER, KNBC, KSBC, X, Y, NID, LISTL, MARKED,
     &      MXNL, MAXNBC, MAXSBC, AMESUR, XNOLD, YNOLD, NXKOLD, MMPOLD,
     &      LINKEG, LISTEG, BMESUR, MLINK, NPROLD, NPNOLD, NPEOLD, NNXK,
     &      REMESH, REXMIN, REXMAX, REYMIN, REYMAX, IDIVIS, SIZMIN,
     &      EMAX, EMIN, NOROOM, ERRCHK, ERR)

      END IF
C
  120 CONTINUE
      RETURN
C
10000 FORMAT (' ** ERROR - DATA PROBLEMS FOR REGION:', I5, ' **')
10010 FORMAT (' ** ERROR - PERIMETER GENERATION ERRORS FOR REGION:'
     &   , I5, ' **')
10020 FORMAT (' ** ERROR - MAPPING BASE GENERATION ERRORS FOR REGION:'
     &   , I5, ' **')
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas