cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  This is the 3D fast multipole method (FMM) code for the Yukawa 
c  potential (also called modified Helmholtz, screened Coulomb, 
c  linearized Poisson-Boltzman, etc.). It is based on the new version 
c  of FMM first introduced by Greengard and Rokhlin for the Laplace
c  equation in 1997 (see Ref. 2 below). It computes the screened 
c  Coulombic interaction between n particles.
c
c  Copyright (C) Jingfang Huang 
c
c  This program is free software; you can redistribute it and/or
c  modify it under the terms of the GNU General Public License
c  as published by the Free Software Foundation; either version 2
c  of the License, or (at your option) any later version.
c
c  This program is distributed in the hope that it will be useful,
c  but WITHOUT ANY WARRANTY; without even the implied warranty of
c  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c  GNU General Public License for more details.
c
c  License details are available in license.txt, you may also write to 
c
c    The Free Software Foundation, Inc., 
c    51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
c
c  Reference:
c    1. Greengard, Leslie F.; Huang, Jingfang: A new version of the 
c       fast multipole method for screened Coulomb interactions in 
c       three dimensions. J. Comput. Phys. 180 (2002), no.2, 642--658. 
c 
c    2. Greengard, Leslie; Rokhlin, Vladimir: A new version of the 
c       fast multipole method for the Laplace equation in three dimensions.  
c       Acta numerica, 1997,  229--269, Acta Numer., 6, Cambridge Univ. 
c       Press, Cambridge, 1997. 
c
c  For suggestions, comments, and bug reports, please contact 
c
c    Jingfang Huang
c    CB# 3250, Phillips Hall
c    Department of Mathematics, UNC
c    Chapel Hill, NC 27599-3250, USA.
c    Email: huang@amath.unc.edu
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  Subroutines for FMMYUK.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YFORMMP(BETA,X0Y0Z0,ZPARTS,CHARGE,NPARTS,
     1                  MPOLE,NTERMS,SCALE,P,C)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine forms the multipole expansion cause by the nparts
c      particles in the box.
c
c  on input:
c
c    beta : the frequency.
c    x0y0z0: center of the expansion
c    nparts: number of sources
c    zparts(3,nparts): array of coordinates of sources
c    charge(nparts): array of strengths of sources
c    nparts: the total number of particles.
c    nterms: order of desired expansion
c    scale: the scaling factor.
c
c  on output:
c
c    mpole: coefficients of multipole expansion
c
c  working space :
c
c    p: used for storing the associate legendre polynomials.
c
c  subroutine called : dsqrt(), in(), lgndr()
c  called from : brfrc()
c
c  note 1: this subroutine needs the precomputed variables c(,)
c          derived from entry frmini()
c
c       2: the multipole expansion is scaled to avoid over- and
c          under-flow.
c
c       3: only the n=0, ... ,nterms, m=0, ..., nterms coefficients
c          are calculated.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS,NPARTS
      REAL *8 SCALE
      REAL *8 BETA,X0Y0Z0(3),ZPARTS(3,NPARTS),CHARGE(NPARTS)
      REAL *8 P(0:NTERMS,0:NTERMS),C(0:NTERMS,0:NTERMS)
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 I,L,M,NCALC
      REAL *8 RK,BI(0:60)
      REAL *8 PRECIS,D,CP,PROJ
      REAL *8 RX,RY,RZ,RR,CTHETA
      COMPLEX *16 IMAG,EPHI(1:60),CPZ
      DATA IMAG/(0.0D0,1.0D0)/
      DATA PRECIS/1.0D-14/
c
      DO 4900 I = 1, NPARTS
        RX = ZPARTS(1,I) - X0Y0Z0(1)
        RY = ZPARTS(2,I) - X0Y0Z0(2)
        RZ = ZPARTS(3,I) - X0Y0Z0(3)
c
c-------compute  distance rr
c         ctheta = cos(theta)
c         ephi(1)=e^(i*phi)
c
        PROJ = RX*RX+RY*RY
        RR = PROJ+RZ*RZ
        PROJ = DSQRT(PROJ)
        D = DSQRT(RR)
c
c-------note: here is a hack. when computing cos(theta) as
c             rz/d, we have to be careful about the possibility
c             of d  being 0 (in which case ctheta is not well
c             defined - we arbitrarily set it to 1.)
c
        IF ( D .LE. PRECIS ) THEN
          CTHETA = 1.0D0
        ELSE
          CTHETA = RZ/D
        ENDIF
c
        IF ( PROJ .LE. PRECIS*D ) THEN
          EPHI(1) = 1.0D0
        ELSE
          EPHI(1) = RX/PROJ - IMAG*RY/PROJ
        ENDIF
c
c-------create array of powers of e^(-i*phi).
c
        DO 4100 L = 2,NTERMS
          EPHI(L) = EPHI(L-1)*EPHI(1)
4100    CONTINUE
c
c-------compute the sperical modified bessel function. note that bi is scaled.
c
        RK=D*BETA
        CALL IN(SCALE,RK,NTERMS,BI, NCALC)
        IF (NCALC .NE. NTERMS+1) THEN
          CALL PRIN2(' bi is, *', BI, NTERMS+1)
          CALL PRIN2(' rk is *', RK, 1)
          CALL PRIN2('scale is *', SCALE, 1)
          CALL PRINF(' the nterms+1 is *', NTERMS+1, 1)
          CALL PRINF(' the ncalc is *', NCALC, 1)
          PRINT *, 'error in forming multipole expansion.'
          STOP
        ENDIF
c
        MPOLE(0,0) = MPOLE(0,0) + CHARGE(I)*BI(0)
c
c-------compute legendre polynomials of argument cos(theta) = ctheta
c         and add contributions from legendre polynomials
c
        CALL LGNDR(NTERMS,CTHETA,P)
c
        DO 4300 L = 1,NTERMS
          CP = CHARGE(I)*P(L,0)*BI(L)*C(L,0)
          MPOLE(L,0) = MPOLE(L,0) + CP
4300    CONTINUE
c
c-------add contributions from associated legendre functions.
c
        DO 4500 L = 1,NTERMS
          DO 4400 M=1,L
            CP = CHARGE(I)*BI(L)*C(L,M)*P(L,M)
            MPOLE(L,M) = MPOLE(L,M) + CP*EPHI(M)
4400      CONTINUE
4500    CONTINUE
4900  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YBRTAEV(BETA,LOCAL,X0Y0Z0,POINT,NTERMS,RPOT,FIELD,
     1  SCALE,P)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    evaluates local expansion at arbitrary point.
c
c  on input:
c
c    beta : the frequency of the equation.
c    local: coefficients of local expansion(scaled)
c    x0y0z0: the center of the expansion
c    point: point of evaluation
c    nterms: order of expansion
c    p: work arrays to hold legendre polynomials
c       and associated legendre functions, respectively.
c
c  on output:
c    rpot: computed potential
c
c  note: the current version will compute the potential only,
c    and the computation of the field will be added later.
c
c   --------------------------------------------------
c  subroutine called :
c
c  called from : brfrc()
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 BETA,SCALE
      REAL *8 X0Y0Z0(3),P(0:NTERMS+1,0:NTERMS+1),POINT(3),RPOT
      REAL *8 FIELD(3)
      COMPLEX *16 LOCAL(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 I,L,M,NTP1,NCALC
      REAL *8 RPOTZ,RX,RY,RZ,PROJ,RR,CTHETA
      REAL *8 PRECIS,D,CP,BI(0:60)
      COMPLEX *16 FIELDTEMP(3)
      COMPLEX *16 PTEMP(0:60,0:60)
      COMPLEX *16 IMAG,EPHI(60),CPZ
c
      DATA IMAG/(0.0D0,1.0D0)/
      DATA PRECIS/1.0D-14/
c
c-----functions called.
c
      REAL *8 DREAL,DSQRT,DBLE
c
      NTP1 = NTERMS+1
c
c-----compute relevant functions of spherical coordinates
c       d = distance*beta, ctheta = cos(theta), ephi = exp(i*phi)
c
      RX = POINT(1) - X0Y0Z0(1)
      RY = POINT(2) - X0Y0Z0(2)
      RZ = POINT(3) - X0Y0Z0(3)
      PROJ = RX*RX+RY*RY
      RR = PROJ+RZ*RZ
      PROJ = DSQRT(PROJ)
      D = DSQRT(RR)
c
      IF (D .LE. PRECIS) THEN
        CTHETA = 0.0D0
      ELSE
        CTHETA = RZ/D
      ENDIF
      IF ( PROJ .LE. PRECIS*D ) THEN
        EPHI(1) = 1.0D0
      ELSE
        EPHI(1) = RX/PROJ + IMAG*RY/PROJ
      ENDIF
      D=D*BETA
c
c-----create array of powers of e^(i*m*phi).
c
      DO 5600 L = 1,NTERMS+2
        EPHI(L+1) = EPHI(L)*EPHI(1)
5600  CONTINUE
c
c-----compute the scaled modified bessel function.
c
      CALL IN(SCALE, D, NTERMS+1, BI, NCALC)
c
      IF (NCALC.NE.NTERMS+2) THEN
        PRINT *,'wrong output from besseli, stop',NTERMS+1,NCALC
        STOP
      ENDIF
c
c-----compute values of legendre functions
c
      CALL LGNDR(NTERMS+1,CTHETA,P)
c
c-----add contribution from 0 mode
c
      RPOT = RPOT + DREAL(LOCAL(0,0))*BI(0)
      PTEMP(0,0)=BI(0)
c
c-----add contributions from legendre polynomials
c
      DO 5700 L = 1,NTERMS
        CP=BI(L)*P(L,0)
        PTEMP(L,0)=CP
        RPOT = RPOT + DREAL(LOCAL(L,0))*CP
5700  CONTINUE
c
c-----add contributions from associated legendre functions.
c
      DO 5900 L = 1,NTERMS
c
c-------compute potential
c
c
c
c ------- print potential to binary file
c
c
        DO 5800 M=1,L
          PTEMP(L,M)=(BI(L)*P(L,M))*EPHI(M)
          RPOT = RPOT + 2.0D0*DREAL(LOCAL(L,M)*PTEMP(L,M))
5800    CONTINUE
5900  CONTINUE

 




c
c-----now calculate the force. see ben et al jcp paper for
c       the formulas.
c
c-----construct a few more terms of ptemp for force field.
c
      DO M=0,NTERMS+1
        PTEMP(NTERMS+1,M)=BI(NTERMS+1)*P(NTERMS+1,M)*EPHI(M)
      ENDDO
c
c-----contribution from (n=0,m=0) term
c
      RPOTZ=LOCAL(0,0)*SCALE
      CPZ=RPOTZ*PTEMP(1,1)
      FIELDTEMP(1)= DCONJG(CPZ)
      FIELDTEMP(2)= RPOTZ*PTEMP(1,0)
      FIELDTEMP(3)=-CPZ
c
c-----contribution from (n=1,m=0)
c
      RPOTZ=LOCAL(1,0)/3.0D0
      CPZ=(RPOTZ*SCALE)*PTEMP(2,1)
      FIELDTEMP(1)=FIELDTEMP(1)+DCONJG(CPZ)
      FIELDTEMP(2)=FIELDTEMP(2)+RPOTZ*(PTEMP(0,0)/SCALE+
     1  2.0D0*PTEMP(2,0)*SCALE)
      FIELDTEMP(3)=FIELDTEMP(3)-CPZ
c
c-----contribution from (n>1,m=0)
c
      DO L=2,NTERMS
        RPOTZ=LOCAL(L,0)/DBLE(2*L+1)
        CPZ=RPOTZ*(PTEMP(L-1,1)/SCALE-PTEMP(L+1,1)*SCALE)
        FIELDTEMP(1)=FIELDTEMP(1)-DCONJG(CPZ)
        FIELDTEMP(2)=FIELDTEMP(2)+RPOTZ*(DREAL(PTEMP(L-1,0))*DBLE(L)/
     1    SCALE+DBLE(L+1)*DREAL(PTEMP(L+1,0))*SCALE)
        FIELDTEMP(3)=FIELDTEMP(3)+CPZ
      ENDDO
c
c-----now update the force field.
c
      FIELD(1)=FIELD(1)-BETA*DREAL(FIELDTEMP(1)-FIELDTEMP(3))/2.0D0
      FIELD(2)=FIELD(2)-BETA*DREAL((FIELDTEMP(1)+FIELDTEMP(3))*IMAG)
     1  /2.0D0
      FIELD(3)=FIELD(3)+BETA*DREAL(FIELDTEMP(2))
c
c-----contribution from (n=1,m=1)
c
      CPZ=LOCAL(1,1)/3.0D0
      FIELDTEMP(1)=CPZ*
     1  ((DREAL(PTEMP(0,0))/SCALE-DREAL(PTEMP(2,0))*SCALE)*2.0D0)
      CPZ=CPZ*SCALE
      FIELDTEMP(2)=CPZ*PTEMP(2,1)
      FIELDTEMP(3)=-CPZ*PTEMP(2,2)
c
c-----contributions for n=2,nterms; m=1,n-2.
c
      DO L=2,NTERMS
        DO M=1,L-2
          CPZ=LOCAL(L,M)/DBLE(2*L+1)
          FIELDTEMP(1)=FIELDTEMP(1)+CPZ*(PTEMP(L-1,M-1)*
     1      (DBLE((L+M-1)*(L+M))/SCALE)-PTEMP(L+1,M-1)*
     2      (DBLE((L-M+1)*(L-M+2))*SCALE))
          FIELDTEMP(2)=FIELDTEMP(2)+CPZ*(PTEMP(L-1,M)*
     1      (DBLE(L+M)/SCALE)+PTEMP(L+1,M)*(DBLE(L-M+1)*SCALE))
          FIELDTEMP(3)=FIELDTEMP(3)+CPZ*(PTEMP(L-1,M+1)
     1      /SCALE-PTEMP(L+1,M+1)*SCALE)
        ENDDO
c
c-------m=n-1.
c
        CPZ=LOCAL(L,L-1)/DBLE(2*L+1)
        FIELDTEMP(1)=FIELDTEMP(1)+CPZ*(PTEMP(L-1,L-2)*
     1    DBLE((L+L-2)*(L+L-1))/SCALE-PTEMP(L+1,L-2)*
     2    6.0D0*SCALE)
        FIELDTEMP(2)=FIELDTEMP(2)+CPZ*(PTEMP(L-1,L-1)*
     1    DBLE(L+L-1)/SCALE+PTEMP(L+1,L-1)*2.0D0*SCALE)
        FIELDTEMP(3)=FIELDTEMP(3)-CPZ*PTEMP(L+1,L)*SCALE
c
c-------m=n.
c
        CPZ=LOCAL(L,L)/DBLE(2*L+1)
        FIELDTEMP(1)=FIELDTEMP(1)+CPZ*(PTEMP(L-1,L-1)*
     1    (DBLE((L+L-1)*(L+L))/SCALE)-PTEMP(L+1,L-1)*
     2    (2.0D0*SCALE))
        FIELDTEMP(2)=FIELDTEMP(2)+CPZ*PTEMP(L+1,L)*SCALE
        FIELDTEMP(3)=FIELDTEMP(3)-CPZ*PTEMP(L+1,L+1)*SCALE
c
      ENDDO
c
c-----now calculate the force.
c
      FIELD(1)=FIELD(1)-BETA*DREAL(FIELDTEMP(1)-FIELDTEMP(3))
      FIELD(2)=FIELD(2)-BETA*DREAL((FIELDTEMP(1)+FIELDTEMP(3))*IMAG)
      FIELD(3)=FIELD(3)+2.0D0*BETA*DREAL(FIELDTEMP(2))
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMPOLETOEXP(MPOLE,NTERMS,NLAMBS,NUMTETS,
     1                      NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c     this subroutine converts a multipole expansion mpole into the
c     corresponding exponential moment function mexp for the
c     both the +z direction and the -z direction.
c
c     u(x,y,z) = sum_{n=0}^{nterms} sum_{m=-n,n}
c                mpole(n,m) y_n^m(cos theta) e^{i m phi}/r^{n+1}
c
c              = (1/2pi) int_0^infty e^{-lambda z}
c                int_0^{2pi} e^{ilambda(xcos(alpha)+ysin(alpha))}
c                mexpup(lambda,alpha) dalpha dlambda
c
c     for +z direction and
c
c              = (1/2pi) int_0^infty e^{lambda z}
c                int_0^{2pi} e^{-ilambda(xcos(alpha)+ysin(alpha))}
c                mexpdown(lambda,alpha) dalpha dlambda
c
c     for -z direction.
c
c     note: the expression for the -z direction corresponds to the
c     mapping (x,y,z) -> (-x,-y,-z), i.e. reflection through the origin.
c     one could also use rotation about the y axis, for which
c     (x,y,z) -> (-x,y,-z) but we stick to the reflected convention.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c     note: the multipole expansion is assumed to have been rescaled
c           so that the box containing sources has unit dimension.
c
c     note: we only store mpole(n,m) for n,m >= 0, since mpole(n,-m)=
c           dconjg(mpole(n,m)). since we store the exponential
c           moment function in the fourier domain (w.r.t. the alpha
c           variable), we compute
c
c       m_lambda(m) = (i)**m sum_{n=m}^n c(n,m) mpole(n,m) lambda^n
c
c           for m >= 0 only, where c(n,m) = 1/sqrt((n+m)!(n-m)!).
c
c       for possible future reference, it should be noted that
c       it is not true that m_lamb(-m) = dconjg(m_lamb(m)).
c       inspection of the integral formula for y_n^{-m} shows that
c       m_lamb(-m) = dconjg(m_lamb(m)) * (-1)**m.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c     on input:
c
c     mpole(0:nterms,0:nterms): the multipole expansion
c
c     rlams(nlambs):  discretization points in lambda integral
c
c     nlambs:         number of discretization pts.
c
c     numtets(nlambs): number of fourier modes needed in expansion
c                    of alpha variable for each lambda value.
c                    note : the numtets is given by numthehalf().
c
c     nexptot =      sum_j numtets(j)
c
c     rlsc() : p_n^m for different lambda_k
c
c     on output:
c
c     mexpf(nexptot): fourier coefficients of the function
c                     mexp(lambda,alpha) for successive discrete
c                     lambda values. they are ordered as follows:
c
c                 mexpf(1,...,numtets(1)) = fourier modes
c                             for lambda_1
c                 mexpf(numtets(1)+1,...,numtets(2)) = fourier modes
c                             for lambda_2
c                 etc.
c     note by huangjf : in return, we will output
c       in mexpup, sum_{n=m}^{nterms} m_n^m*p_n^m. (all are scaled)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c     note by jingfang :
c     1.this subroutine will compute the inner sum.
c       instead of compute all the nterms modes, only
c       the necessary modes are calculated. the number of mode
c       needed are provided by the subroutine numthehalf().
c       note that the number is always less than nterms needed,
c       and the worst case is nterms+1?.
c
c     2.
c       subroutine called :
c       called from : mkudexp(), mknsexp(), mkewexp()
c
c     3. the down-list will have the same fourier modes as the
c        up-list if we only change the sign of the z. so we don't need to
c        compute them separately.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS,NLAMBS,NUMTETS(NLAMBS),NEXPTOT
      REAL *8 RLSC(0:NTERMS,0:NTERMS,NLAMBS)
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MEXPUP(NEXPTOT), MEXPDOWN(NEXPTOT)
c
c-----local variables.
c
      INTEGER *4 NTOT,NL,MTH,NCURRENT,NM
      REAL *8 SGN
      COMPLEX *16 ZTMP1,ZTMP2
c
c-----loop over multipole order to generate mexpup and mexpdown values.
c
      NTOT = 0
      DO NL = 1,NLAMBS
        SGN = -1.0D0
        DO MTH = 0,NUMTETS(NL)-1
          NCURRENT = NTOT+MTH+1
          ZTMP1 = 0.0D0
          ZTMP2 = 0.0D0
          SGN = -SGN
          DO NM = MTH,NTERMS,2
           ZTMP1 = ZTMP1 +
     1       RLSC(NM,MTH,NL)*MPOLE(NM,MTH)
          ENDDO
          DO NM = MTH+1,NTERMS,2
            ZTMP2 = ZTMP2 +
     1        RLSC(NM,MTH,NL)*MPOLE(NM,MTH)
          ENDDO
          MEXPUP(NCURRENT) = (ZTMP1 + ZTMP2)
          MEXPDOWN(NCURRENT) = SGN*(ZTMP1 - ZTMP2)
        ENDDO
        NTOT = NTOT+NUMTETS(NL)
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YRLSCINI(SCAL,BETA,RLSC,NLAMBS,RLAMS,NTERMS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    precomputes the coefficients for mpoletoexp for different levels.
c    note: for different level, this subroutine should be
c          called because of the change of the scaled frequency.
c    in return, the rlsc will contain the associated legendre polynomial
c    p_n^m( (beta+lamda_i)/beta )*scal^n
c
c  on input :
c    scal : the scaling factor for the p_n^m. otherwise p_n^m will
c           overflow.
c    beta : the scaled frequency/beta for the current level.
c    nlambs : the total number of lamdas for the first integral.
c    rlams : the nodes/lamdas.
c    nterms : the total number of terms in the multipole expansion.
c
c  on output :
c    rlsc(n,m,nlamda) : the required information for the whole level.
c
c  subroutine called :
c
c  called from : brfrc()
c
c  note : the rlsc() will be used in the subroutine mpoletoexp().
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      IMPLICIT NONE
c
      INTEGER *4 NLAMBS, NTERMS
      REAL *8 RLSC(0:NTERMS,0:NTERMS,NLAMBS)
      REAL *8 RLAMS(NLAMBS)
      REAL *8 SCAL, BETA
c
c-----local variables.
c
      INTEGER *4 NL
      REAL *8 U1
c
      DO NL = 1,NLAMBS
        U1=RLAMS(NL)/BETA+1.0D0
        CALL LGNDRGT1(SCAL, NTERMS, U1, RLSC(0,0,NL))
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YFTOPHYS(MEXPF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1                      NTHMAX,MEXPPHYS,FEXPE,FEXPO)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine evaluates the fourier expansion of the
c    exponential moment function m(\lambda,\alpha) at equispaced
c    nodes.
c
c  on input:
c
c    mexpf(*):     fourier coefficients of the function
c                  mexp(lambda,alpha) for discrete lambda values.
c                  they are ordered as follows:
c
c               mexpf(1,...,numfour(1)) = fourier modes for lambda_1
c               mexpf(numfour(1)+1,...,numfour(2)) = fourier modes
c                                              for lambda_2
c               etc.
c
c    nlambs:        number of discretization pts. in lambda integral
c    rlams(nlambs): discretization points in lambda integral.
c    numfour(j):   number of fourier modes in the expansion
c                      of the function m(\lambda_j,\alpha)
c    nthmax =      max_j numfour(j)
c    numphys : number of fourier modes in the plane wave expansion.
c    fexpe =      precomputed array of exponentials needed for
c                 fourier series evaluation. even terms.
c    fexpo =      precomputed array of exponentials needed for
c                 fourier series evaluation. odd terms.
c  note : we will keep these two terms because in
c         the helmholtz equation, we will need these.
c         however, in yukawa, it is not necessary to have
c         them separated.--huangjf
c
c  on output:
c    mexpphys(*):  discrete values of the moment function
c                  m(\lambda,\alpha), ordered as follows.
c
c        mexpphys(1),...,mexpphys(numphys(1)) = m(\lambda_1,0),...,
c             m(\lambda_1, 2*pi*(numphys(1)-1)/numphys(1)).
c        mexpphys(numphys(1)+1),...,mexpphys(numphys(2)) =
c             m(\lambda_2,0),...,
c                 m(\lambda_2, 2*pi*(numphys(2)-1)/numphys(2)).
c        etc.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c      this subroutine computes the outer sum, it is possible
c      to do this using fft. but the current version will not do that.
c
c  subroutine called :
c
c  called from : mkudexp(), mknsexp(), mkewexp()
c
c  note :
c    the current subroutine computes sum_{m=-numfour, numfour}
c      e^{im*alpha} * i^|m| *inner(m)
c
c    the constant will left to the pw_local.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLAMBS,NUMFOUR(NLAMBS),NUMPHYS(NLAMBS),NTHMAX
      REAL *8 RLAMS(NLAMBS)
c
      COMPLEX *16 MEXPF(1)
      COMPLEX *16 MEXPPHYS(1)
      COMPLEX *16 FEXPE(1)
      COMPLEX *16 FEXPO(1)
c
c-----local variables.
c
      INTEGER *4 I,IVAL,MM
      INTEGER *4 NFTOT,NPTOT,NEXTE,NEXTO
      REAL *8 SGN,PI,RTMP
      COMPLEX *16 IMA
c
c-----functions
c
      REAL *8 DATAN
c
      DATA IMA/(0.0D0,1.0D0)/
c
      PI=DATAN(1.0D0)*4.0D0
c
      NFTOT = 0
      NPTOT  = 0
      NEXTE = 1
      NEXTO = 1
      DO 2000 I=1,NLAMBS
        DO 1200 IVAL=1,NUMPHYS(I)/2
          MEXPPHYS(NPTOT+IVAL) = MEXPF(NFTOT+1)
          SGN=-2.0D0
          DO MM = 2,NUMFOUR(I),2
            SGN=-SGN
            RTMP = SGN*DBLE( FEXPE(NEXTE)*MEXPF(NFTOT+MM) )
            NEXTE = NEXTE + 1
            MEXPPHYS(NPTOT+IVAL) = MEXPPHYS(NPTOT+IVAL) +
     1        DCMPLX(0.0D0,RTMP)
          ENDDO
c
          SGN=2.0D0
          DO MM = 3,NUMFOUR(I),2
            SGN=-SGN
            RTMP = SGN*DBLE( FEXPO(NEXTO)*MEXPF(NFTOT+MM) )
            NEXTO = NEXTO + 1
            MEXPPHYS(NPTOT+IVAL) = MEXPPHYS(NPTOT+IVAL) +
     1        RTMP
          ENDDO
1200    CONTINUE
c
        NFTOT = NFTOT+NUMFOUR(I)
        NPTOT = NPTOT+NUMPHYS(I)/2
2000  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPHYSTOF(MEXPF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1                      NTHMAX,MEXPPHYS,FEXPBACK)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine converts the discretized exponential moment function
c    into its fourier expansion.
c    it calculates the inner sum of the exp->local expansion.
c      (/sum_{j=1}^{m(k)} w(k,j)*e^{-im*alpha_j})/m(k) for k=1, nlambs,
c      and m=0, numfour.
c      numfour is the total number of the fourier modes.
c      or in other words, those l_n^m <>0.
c      the summation is over the numphys.
c
c  on input:
c
c    mexpphys(*):  discrete values of the moment function
c                  m(\lambda,\alpha), ordered as follows.
c
c        mexpphys(1),...,mexpphys(numphys(1)) = m(\lambda_1,0),...,
c             m(\lambda_1, 2*pi*(numphys(1)-1)/numphys(1)).
c        mexpphys(numphys(1)+1),...,mexpphys(numphys(2)) =
c             m(\lambda_2,0),...,
c                 m(\lambda_2, 2*pi*(numphys(2)-1)/numphys(2)).
c        etc.
c
c    nlambs:        number of discretization pts. in lambda integral
c    rlams(nlambs): discretization points in lambda integral.
c    numfour(j):   number of fourier modes in the expansion
c                      of the function m(\lambda_j,\alpha)
c    nthmax =      max_j numfour(j)
c    fexpback : contains the precomputed e^{-im *alpha_j}
c
c  on output:
c
c    mexpf(*):     fourier coefficients of the function
c                  mexp(lambda,m) for discrete lambda values.
c                  they are ordered as follows:
c
c               mexpf(1,...,numfour(1)) = fourier modes for lambda_1
c               mexpf(numfour(1)+1,...,numfour(2)) = fourier modes
c                                              for lambda_2
c               etc.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLAMBS,NUMFOUR(NLAMBS),NUMPHYS(NLAMBS),NTHMAX
      REAL *8 RLAMS(NLAMBS)
c
      COMPLEX *16 MEXPF(1)
      COMPLEX *16 MEXPPHYS(1)
      COMPLEX *16 FEXPBACK(1)
c
c-----local variables.
c
      INTEGER *4 I,IVAL,MM
      INTEGER *4 NFTOT,NPTOT,NEXT,NALPHA,NALPHA2
      REAL *8 RTMP,PI,HALPHA
      COMPLEX *16 IMA,ZTMP
c
c-----functions
c
      REAL *8 DATAN
c
      DATA IMA/(0.0D0,1.0D0)/
c
      PI=DATAN(1.0D0)*4.0D0
      NFTOT = 0
      NPTOT = 0
      NEXT = 1
c
      DO 2000 I=1,NLAMBS
        NALPHA = NUMPHYS(I)
        NALPHA2 = NALPHA/2
        HALPHA=2.0D0*PI/DBLE(NALPHA)
c
c-------first mm=0 case.
c
        MEXPF(NFTOT+1) = 0.0D0
        DO IVAL=1,NALPHA2
          MEXPF(NFTOT+1) = MEXPF(NFTOT+1) +
     1      2.0D0*DBLE(MEXPPHYS(NPTOT+IVAL))
        ENDDO
        MEXPF(NFTOT+1) = MEXPF(NFTOT+1)/DBLE(NALPHA)
c
c-------even mm.  (w(k,j)+conj(w(k,j))*e^{-im alpha_j}
c
        DO MM = 3,NUMFOUR(I),2
          MEXPF(NFTOT+MM) = 0.0D0
          DO IVAL=1,NALPHA2
            RTMP = 2.0D0*DREAL(MEXPPHYS(NPTOT+IVAL))
            MEXPF(NFTOT+MM) = MEXPF(NFTOT+MM) +
     1        FEXPBACK(NEXT)*RTMP
            NEXT = NEXT+1
          ENDDO
          MEXPF(NFTOT+MM) = MEXPF(NFTOT+MM)/DBLE(NALPHA)
        ENDDO
c
c-------odd mm. (w(k,j)-conj(w(k,j))*e^{-im alpha_j}
c
        DO MM = 2,NUMFOUR(I),2
          MEXPF(NFTOT+MM) = 0.0D0
          DO IVAL=1,NALPHA2
            ZTMP = 2.0D0*DCMPLX(0.0D0,DIMAG(MEXPPHYS(NPTOT+IVAL)))
            MEXPF(NFTOT+MM) = MEXPF(NFTOT+MM) +
     1        FEXPBACK(NEXT)*ZTMP
            NEXT = NEXT+1
          ENDDO
          MEXPF(NFTOT+MM) = MEXPF(NFTOT+MM)/DBLE(NALPHA)
        ENDDO
c
        NFTOT = NFTOT+NUMFOUR(I)
        NPTOT = NPTOT+NUMPHYS(I)/2
2000  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMPSHIFT(IFL,MPOLE,MPOLEN,MARRAY,NTERMS,DC,RD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine shifts the center of a child box multipole
c      expansion to the parent center, via the rotation scheme.
c      we rotate the coordinate system, shift along the z-axis,
c      and then rotate back.
c      there are eight possible child locations, defined in the
c      calling sequence of this routine by the parameters ifl
c      and rd (see below).
c
c  on input:
c
c     integer *4  ifl    = flag which describes the quadrant in which
c                          the child box lies (1,2,3 or 4).
c     complex *16 mpole  = coefficients of original multipole exp.
c     integer *4  nterms = integer indicates the terms retained in the
c                          expansion.
c     real *8     dc     = precomputed array containing
c                          the shifting coefficients
c                          along the z-axis. this is precomputed by
c                          the subroutine mpshftcoef() at the beginning
c                          of different levels.
c     real *8     rd     = precomputed array containing rotation matrix
c                          about y-axis.
c                 there are two possible y rotations, depending on
c                 whether the child box lies in +z half space or the
c                 -z half space. they are referred to in the calling
c                 program as rdp and rdm, respectively.
c                 this is precomputed in the subroutine rotgen().
c
c     complex *16 marray = work array
c
c  on output:
c
c     complex *16 mpolen = coefficients of shifted multipole exp.
c
c     note 1 : the rotation part is the same as the old subroutine
c              of the laplace equation. the shifting along the z-axis
c              is changed to the new version. the rotation matrix
c              is precomputed at the very beginning and the shifting
c              matrix is computed at different levels.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 IFL,NTERMS
c
      REAL *8 DC(0:NTERMS,0:NTERMS,0:NTERMS)
      REAL *8 RD(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
c
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MPOLEN(0:NTERMS,0:NTERMS)
      COMPLEX *16 MARRAY(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M,JNEW,KNEW,MP
      REAL *8 ARG
      COMPLEX *16 EPHI(0:60),IMAG
      DATA IMAG/(0.0D0,1.0D0)/
c
c-----functions called.
c
      COMPLEX *16 DCONJG
      COMPLEX *16 DCMPLX
c
      EPHI(0)=1.0D0
      ARG = DSQRT(2.0D0)/2.0D0
      IF (IFL.EQ.1) THEN
        EPHI(1) = DCMPLX(-ARG,ARG)
      ELSE IF (IFL.EQ.2) THEN
        EPHI(1) = DCMPLX(ARG,ARG)
      ELSE IF (IFL.EQ.3) THEN
        EPHI(1) = DCMPLX(ARG,-ARG)
      ELSE IF (IFL.EQ.4) THEN
        EPHI(1) = DCMPLX(-ARG,-ARG)
      ELSE
        CALL PRINF('error in lolo with ifl = *',IFL,1)
      ENDIF
c
c-----create array of powers of e^(i*m*phi).
c
      DO L = 1,NTERMS+1
        EPHI(L+1) = EPHI(L)*EPHI(1)
      ENDDO
c
c-----a rotation of phi radians about the z-axis in the
c       original coordinate system.
c
      DO L=0,NTERMS
        DO M=0,L
          MPOLEN(L,M)=DCONJG(EPHI(M))*MPOLE(L,M)
        ENDDO
      ENDDO
c
c-----a rotation about the y'-axis  in the rotated system.
c
      DO L=0,NTERMS
        DO M=0,L
          MARRAY(L,M)=MPOLEN(L,0)*RD(L,0,M)
          DO MP=1,L
            MARRAY(L,M)=MARRAY(L,M)+MPOLEN(L,MP)*RD(L,MP,M)+
     1        DCONJG(MPOLEN(L,MP))*RD(L,MP,-M)
          ENDDO
        ENDDO
      ENDDO
c
c-----shift along z-axis.
c       note that everything is scaled.
c
      DO JNEW=0,NTERMS
        DO KNEW=0,JNEW
          MPOLEN(JNEW,KNEW)=0.0D0
          DO L=KNEW, NTERMS
            MPOLEN(JNEW,KNEW)=MPOLEN(JNEW,KNEW)+MARRAY(L,KNEW)*
     1        DC(KNEW,JNEW,L)
          ENDDO
        ENDDO
      ENDDO
c
c-----reverse rotation about the y'-axis.
c
      DO L=0,NTERMS
        DO M=0,L,2
          MARRAY(L,M)=MPOLEN(L,0)*RD(L,0,M)
          DO MP=1,L,2
            MARRAY(L,M)=MARRAY(L,M)-(MPOLEN(L,MP)*RD(L,MP,M)+
     1        DCONJG(MPOLEN(L,MP))*RD(L,MP,-M))
          ENDDO
          DO MP=2,L,2
            MARRAY(L,M)=MARRAY(L,M)+(MPOLEN(L,MP)*RD(L,MP,M)+
     1        DCONJG(MPOLEN(L,MP))*RD(L,MP,-M))
          ENDDO
        ENDDO
        DO M=1,L,2
          MARRAY(L,M)=-MPOLEN(L,0)*RD(L,0,M)
          DO MP=1,L,2
            MARRAY(L,M)=MARRAY(L,M)+(MPOLEN(L,MP)*RD(L,MP,M)+
     1        DCONJG(MPOLEN(L,MP))*RD(L,MP,-M))
          ENDDO
          DO MP=2,L,2
            MARRAY(L,M)=MARRAY(L,M)-(MPOLEN(L,MP)*RD(L,MP,M)+
     1        DCONJG(MPOLEN(L,MP))*RD(L,MP,-M))
          ENDDO
        ENDDO
      ENDDO
c
c-----rotate back phi radians about the z-axis in the above system.
c
      DO L=0,NTERMS
        DO M=0,L
          MPOLEN(L,M)=EPHI(M)*MARRAY(L,M)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YLCSHIFT(IFL,LOCAL,LOCALN,MARRAY,NTERMS,DC,RD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine shifts the local expansion of a parent cell
c     to the center of one of its children, via the rotation scheme.
c     that is, we rotate the coordinate system, shift along the z-axis,
c     and then rotate back.
c     there are eight possible child locations, defined in the
c     calling sequence of this routine by the parameters ifl
c     and rd (see below).
c
c on input:
c
c     integer *4  ifl    = flag which describes the quadrant in which
c                          the child box lies (1,2,3 or 4).
c     complex *16 local  = coefficients of original multipole exp.
c     integer *4  nterms = integer indicates the terms retained in the
c                          expansion.
c     real *8     dc     = precomputed array containing coefficients
c                          for the local translation along the z axis.
c                          this is precomputed by the subroutine
c                          lcshftcoef()
c     real *8     rd     = precomputed array containing rotation matrix
c                          about y-axis.
c                 there are two possible y rotations, depending on
c                 whether the child box lies in +z half space or the
c                 -z half space. they are referred to in the calling
c                 program as rdp and rdm, respectively.
c                          this is precomputed by the subroutine
c                          rotgen().
c     complex *16 marray = work array
c
c on output:
c
c     complex *16 localn = coefficients of shifted multipole exp.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 IFL,NTERMS
      REAL *8 DC(0:NTERMS,0:NTERMS,0:NTERMS)
      REAL *8 RD(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
c
      COMPLEX *16 LOCAL(0:NTERMS,0:NTERMS)
      COMPLEX *16 LOCALN(0:NTERMS,0:NTERMS)
      COMPLEX *16 MARRAY(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M,MP,JNEW,KNEW
      REAL *8 ARG
      COMPLEX *16 EPHI(0:60),IMAG
c
c-----functions.
c
      REAL *8 DSQRT
      COMPLEX *16 DCOMPLX,DCONJG
c
      DATA IMAG/(0.0D0,1.0D0)/
c
      EPHI(0)=1.0D0
      ARG = DSQRT(2.0D0)/2.0D0
      IF (IFL.EQ.1) THEN
        EPHI(1) = DCMPLX(ARG,-ARG)
      ELSE IF (IFL.EQ.2) THEN
        EPHI(1) = DCMPLX(-ARG,-ARG)
      ELSE IF (IFL.EQ.3) THEN
        EPHI(1) = DCMPLX(-ARG,ARG)
      ELSE IF (IFL.EQ.4) THEN
        EPHI(1) = DCMPLX(ARG,ARG)
      ELSE
        CALL PRINF('error in lolo with ifl = *',IFL,1)
      ENDIF
c
c----- create array of powers of r and e^(i*m*phi).
c
      DO L = 1,NTERMS+1
        EPHI(L+1) = EPHI(L)*EPHI(1)
      ENDDO
c
c----- a rotation of phi radians about the z-axis in the
c      original coordinate system.
c
      DO L=0,NTERMS
        DO M=0,L
          LOCALN(L,M)=DCONJG(EPHI(M))*LOCAL(L,M)
        ENDDO
      ENDDO
c
c-----a rotation about the y'-axis to align z' axis.
c
      DO L=0,NTERMS
        DO M=0,L
          MARRAY(L,M)=LOCALN(L,0)*RD(L,0,M)
          DO MP=1,L
            MARRAY(L,M)=MARRAY(L,M)+LOCALN(L,MP)*RD(L,MP,M)+
     1        DCONJG(LOCALN(L,MP))*RD(L,MP,-M)
          ENDDO
        ENDDO
      ENDDO
c
c-----shift along z'-axis.
c
      DO JNEW= 0,NTERMS
        DO KNEW=0,JNEW
          LOCALN(JNEW,KNEW)=0.0D0
          DO L=KNEW,NTERMS
            LOCALN(JNEW,KNEW)=LOCALN(JNEW,KNEW)+MARRAY(L,KNEW)*
     1        DC(KNEW,JNEW,L)
          ENDDO
        ENDDO
      ENDDO
c
c-----rotate back about the y'-axis.
c
      DO L=0,NTERMS
        DO M=0,L,2
          MARRAY(L,M)=LOCALN(L,0)*RD(L,0,M)
          DO MP=1,L,2
            MARRAY(L,M)=MARRAY(L,M)-(LOCALN(L,MP)*RD(L,MP,M)+
     1        DCONJG(LOCALN(L,MP))*RD(L,MP,-M))
          ENDDO
          DO MP=2,L,2
            MARRAY(L,M)=MARRAY(L,M)+(LOCALN(L,MP)*RD(L,MP,M)+
     1        DCONJG(LOCALN(L,MP))*RD(L,MP,-M))
          ENDDO
        ENDDO
        DO M=1,L,2
          MARRAY(L,M)=-LOCALN(L,0)*RD(L,0,M)
          DO MP=1,L,2
            MARRAY(L,M)=MARRAY(L,M)+(LOCALN(L,MP)*RD(L,MP,M)+
     1        DCONJG(LOCALN(L,MP))*RD(L,MP,-M))
          ENDDO
          DO MP=2,L,2
            MARRAY(L,M)=MARRAY(L,M)-(LOCALN(L,MP)*RD(L,MP,M)+
     1        DCONJG(LOCALN(L,MP))*RD(L,MP,-M))
          ENDDO
        ENDDO
      ENDDO
c
c-----rotate back about the z-axis.
c
      DO L=0,NTERMS
        DO M=0,L
          LOCALN(L,M)=EPHI(M)*MARRAY(L,M)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMKEXPS(BETA,RLAMS,NLAMBS,NUMPHYS,NEXPTOTP,XS,YS,ZS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine computes the tables of exponentials needed
c     for translating exponential representations of harmonic
c     functions, discretized via norman's quadratures.
c
c     u   = \int_0^\infty e^{-(lambda+beta) z}
c     \int_0^{2\pi} e^{i \dsqrt(lambda*lambda+2*beta*lambda)
c                         (x cos(u)+y sin(u))}
c           mexpphys(lambda,u) du dlambda
c
c     mexpphys(*):  discrete values of the moment function
c                   m(\lambda,u), ordered as follows.
c
c         mexpphys(1),...,mexpphys(numphys(1)) = m(\lambda_1,0),...,
c              m(\lambda_1, 2*pi*(numphys(1)-1)/numphys(1)).
c         mexpphys(numphys(1)+1),...,mexpphys(numphys(2)) =
c              m(\lambda_2,0),...,
c                  m(\lambda_2, 2*pi*(numphys(2)-1)/numphys(2)).
c         etc.
c         note : in the current version, only half of the modes are
c                stored because of the symmetry of u and u+pi.
c
c  on input:
c
c     beta : the scaled frequency.
c     rlams(nlambs)   discretization points in lambda integral
c     nlambs          number of discret. pts. in lambda integral
c     numphys(j)     number of nodes in u integral needed
c                    for corresponding lambda =  lambda_j.
c     nexptotp        sum_j numphys(j)
c
c  on output:
c
c        define w1=\lambda_j+beta, and w2= sqrt(lambda_j**2+2*beta*lambda_j)
c     xs(1,nexptotp)   e^{i w2 (cos(u_k)}  in above ordering
c     xs(2,nexptotp)   e^{i w2 (2 cos(u_k)}  in above ordering.
c     xs(3,nexptotp)   e^{i w2 (3 cos(u_k)}  in above ordering.
c     ys(1,nexptotp)   e^{i w2 (sin(u_k)}  in above ordering.
c     ys(2,nexptotp)   e^{i w2 (2 sin(u_k)}  in above ordering.
c     ys(3,nexptotp)   e^{i w2 (3 sin(u_k)}  in above ordering.
c     zs(1,nexptotp)   e^{-w1}     in above ordering.
c     zs(2,nexptotp)    e^{-2 w1}   in above ordering.
c     zs(3,nexptotp)    e^{-3 w1}   in above ordering.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLAMBS,NUMPHYS(NLAMBS),NEXPTOTP
c
      REAL *8 BETA
      REAL *8 RLAMS(NLAMBS)
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 NTOT,NL,MTH,NCURRENT
      REAL *8 W1,W2,PI,HU,U
      COMPLEX *16 IMA
c
c-----functions.
c
      REAL *8 DATAN,DSQRT,DEXP,DCOS,DSIN
      COMPLEX *16 CDEXP
      DATA IMA/(0.0D0,1.0D0)/
c
c-----loop over each lambda value
c
      PI = 4D0*DATAN(1.0D0)
      NTOT = 0
      DO NL = 1,NLAMBS
         W1= RLAMS(NL)+BETA
         W2=DSQRT( RLAMS(NL)*(RLAMS(NL)+BETA*2.0D0) )
         HU=2.0D0*PI/DBLE(NUMPHYS(NL))
         DO MTH = 1,NUMPHYS(NL)/2
            U = DBLE(MTH-1)*HU
            NCURRENT = NTOT+MTH
            ZS(1,NCURRENT) = DEXP( -W1 )
            ZS(2,NCURRENT) = ZS(1,NCURRENT)*ZS(1,NCURRENT)
            ZS(3,NCURRENT) = ZS(2,NCURRENT)*ZS(1,NCURRENT)
            XS(1,NCURRENT) = CDEXP(IMA*W2*DCOS(U))
            XS(2,NCURRENT) = XS(1,NCURRENT)*XS(1,NCURRENT)
            XS(3,NCURRENT) = XS(2,NCURRENT)*XS(1,NCURRENT)
            YS(1,NCURRENT) = CDEXP(IMA*W2*DSIN(U))
            YS(2,NCURRENT) = YS(1,NCURRENT)*YS(1,NCURRENT)
            YS(3,NCURRENT) = YS(2,NCURRENT)*YS(1,NCURRENT)
         ENDDO
         NTOT = NTOT+NUMPHYS(NL)/2
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMKFEXP(NLAMBS,NUMFOUR,NUMPHYS,FEXPE,FEXPO,FEXPBACK)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     precomputes the e^(im*alpha) needed in mp->pw->local.
c
c     this subroutine computes the tables of exponentials needed
c     for mapping from fourier to physical domain.
c     in order to minimize storage, they are organized in a
c     one-dimenional array corresponding to the order in which they
c     are accessed by subroutine ftophys.
c
c     size of fexpe, fexpo =          40000   for nlambs = 39
c     size of fexpe, fexpo =          15000   for nlambs = 30
c     size of fexpe, fexpo =           4000   for nlambs = 20
c     size of fexpe, fexpo =            400   for nlambs = 10
c
c
c  on input :
c
c       nlambs : the total number of nodes for the outer integral.
c       numfour : contains the number of nodes for the fourier
c                 representation.
c       numphys : contains the number of nodes for the inner integral
c                 for the plane wave expansion.
c
c  on output :
c
c       fexpe : the exponentials for the fourier modes. e^(im*alpha)
c               where m is all the fourier modes and alpha comes
c               from the physical modes. odd terms.
c       fexpo : the even terms.
c       fexpback : the exponentials used for the translation
c                  from plane wave to local.
c
c     functions called :
c
c     called from :
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4  NLAMBS,NUMPHYS(NLAMBS),NUMFOUR(NLAMBS)
c
      COMPLEX *16 FEXPE(1)
      COMPLEX *16 FEXPO(1)
      COMPLEX *16 FEXPBACK(1)
c
c-----local variables.
c
      INTEGER *4 I,J,MM,NEXTE,NEXTO,NEXT,NALPHA,NALPHA2
      REAL *8 HALPHA,ALPHA,PI
      COMPLEX *16 IMA
      DATA IMA/(0.0D0,1.0D0)/
c
      PI = 4.0D0*DATAN(1.0D0)
      NEXTE = 1
      NEXTO = 1
c
      DO I=1,NLAMBS
        NALPHA = NUMPHYS(I)
        NALPHA2 = NALPHA/2
        HALPHA=2.0D0*PI/DBLE(NALPHA)
c
        DO J=1,NALPHA2
          ALPHA=DBLE(J-1)*HALPHA
          DO MM = 2,NUMFOUR(I),2
            FEXPE(NEXTE) = CDEXP(IMA*DBLE(MM-1)*ALPHA)
            NEXTE = NEXTE + 1
          ENDDO
c
          DO MM = 3,NUMFOUR(I),2
            FEXPO(NEXTO) = CDEXP(IMA*DBLE(MM-1)*ALPHA)
            NEXTO = NEXTO + 1
          ENDDO
        ENDDO
      ENDDO
c
      NEXT = 1
      DO I=1,NLAMBS
        NALPHA = NUMPHYS(I)
        NALPHA2 = NALPHA/2
        HALPHA=2*PI/NALPHA
        DO MM = 3,NUMFOUR(I),2
          DO J=1,NALPHA2
            ALPHA=(J-1)*HALPHA
            FEXPBACK(NEXT)  = CDEXP(-IMA*(MM-1)*ALPHA)
            NEXT = NEXT + 1
          ENDDO
        ENDDO
c
        DO MM = 2,NUMFOUR(I),2
          DO J=1,NALPHA2
            ALPHA=(J-1)*HALPHA
            FEXPBACK(NEXT)  = CDEXP(-IMA*(MM-1)*ALPHA)
            NEXT = NEXT + 1
          ENDDO
        ENDDO
      ENDDO
c
c      print *, '# of terms of fexpe, fexpo, fexpback', nexte,nexto,next
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMPSHFTCOEF(SCAL,BETA,R0,NTERMS,C0,INFO)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine precomputes the coefficients for the
c      shifting of the multipole expansion.
c
c  on input :
c
c    scal : the scale factor of the child level. the parent
c           has the scaling factor of 2*scal
c    beta : the coefficient for the helmholtz equation.
c    r0 : the real shifting distance.
c    nterms : number of terms in the multipole expansion.
c
c  on output :
c
c    c0() : the coefficients of the shifting of the multipole expansion.
c    info : error messages.
c
c  workspace :
c    fact(0:200): assigned locally.
c    bj(0:200): assigned locally.
c
c  function called : in().
c
c  note : the new scaled multipole expansion coefficients
c         m_l^m= sum_{n=m}^{nterms} m_n^m *c(m,l,n)
c         where m_n^m is also the scaled multipole expansion
c         of the child (by scal^n).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS,INFO
      REAL *8 BETA
      REAL *8 C0(0:NTERMS, 0:NTERMS, 0:NTERMS)
      REAL *8 R0,SCAL
c
c-----local variables.
c
      INTEGER *4 I,MNEW,LNEW,NN,NP,NCALC
      REAL *8 FACT(0:200),R0K
      REAL *8 BJ(0:200)
c
      FACT(0)=1.0D0
      DO 1 I=1, 2*NTERMS+1
        FACT(I)=FACT(I-1)*DBLE(I)
1     CONTINUE
c
      R0K=R0*BETA
      CALL IN(R0K, R0K, 2*NTERMS+2, BJ, NCALC)
c
      DO 2 MNEW=0,NTERMS
        DO 3 LNEW=MNEW, NTERMS
          DO 10 NN=MNEW, NTERMS
            C0(MNEW, LNEW, NN)=0.0D0
            DO 20 NP=MNEW, MIN(NN,LNEW)
              C0(MNEW, LNEW, NN) =C0(MNEW, LNEW, NN)
     1   +SCAL**(NN-LNEW)*(2.0D0)**(-LNEW-NP)*(-1.0D0)**(LNEW+NN)*
     2        DBLE(2*LNEW+1)*FACT(LNEW-MNEW)/FACT(NP+MNEW)*
     3        FACT(NN+MNEW)*FACT(2*NP)/FACT(NP)/
     4        FACT(NP-MNEW)/FACT(LNEW-NP)
     5        /FACT(NN-NP)*BJ(LNEW+NN-NP)*(R0K)**(LNEW+NN-2*NP)
20          CONTINUE
10        CONTINUE
3       CONTINUE
2     CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YLCSHFTCOEF(SCAL,BETA,R0,NTERMS,C0,INFO)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine precomputes the coefficients for the
c      shifting of the local expansion.
c
c  on input:
c    scal: the scale factor.
c    beta: the coefficient for the helmholtz equation.
c    nterms: number of terms in the local expansion.
c
c  on output :
c    c0(): the coefficients of the shifting of the local expansion.
c    info: error messages.
c
c  workspace:
c    none
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS,INFO
      REAL *8 SCAL,BETA,R0
      REAL *8 C0(0:NTERMS, 0:NTERMS, 0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 I,MNEW,LNEW,NN,NP,NCALC
      REAL *8 R0K
      REAL *8 FACT(0:200)
      REAL *8 BJ(0:200)
c
      FACT(0)=1.0D0
      DO 1 I=1, 2*NTERMS+1
        FACT(I)=FACT(I-1)*DBLE(I)
1     CONTINUE
c
      R0K=R0*BETA
      CALL IN(R0K, R0K, 2*NTERMS+2, BJ, NCALC)
c
      DO 2 MNEW=0,NTERMS
        DO 3 LNEW=MNEW, NTERMS
          DO 10 NN=MNEW, NTERMS
            C0(MNEW, LNEW, NN)=0.0D0
            DO 20 NP=MNEW, MIN(NN,LNEW)
             C0(MNEW, LNEW, NN) =C0(MNEW, LNEW, NN)
     1         +SCAL**(LNEW-NN)*(2.0D0)**(-LNEW-NP)*
     2         DBLE(2*LNEW+1)*FACT(LNEW-MNEW)/FACT(NP+MNEW)*
     3         FACT(NN+MNEW)*FACT(2*NP)/FACT(NP)/
     4         FACT(NP-MNEW)/FACT(LNEW-NP)
     5         /FACT(NN-NP)*BJ(LNEW+NN-NP)*(R0K)**(LNEW+NN-2*NP)
20          CONTINUE
10        CONTINUE
3       CONTINUE
2     CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YFORMLC(BETA,X0Y0Z0,ZPARTS,CHARGE,NPARTS,
     1                  LOCAL,NTERMS,SCALE,P,C)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    this subroutine forms the local expansion caused by the nparts
c      particles in the box.
c
c  on input:
c    beta: the frequency.
c    x0y0z0: center of the expansion
c    nparts: number of sources
c    zparts(3,nparts): array of coordinates of sources
c    charge(nparts): array of strengths of sources
c    nparts: the total number of particles.
c    nterms: order of desired expansion
c    scale: the scaling factor.
c
c  on output:
c    local: coefficients of local expansion
c
c  working space :
c    p : used for storing the associate legendre polynomials.
c
c  subroutine called : dsqrt(), in(), lgndr()
c
c  called from : brfrc()
c
c  note 1: this subroutine needs the precomputed variables c(,)
c          derived from entry frmini()
c
c       2: the multipole expansion is scaled to avoid over- and
c          under-flow.
c
c       3: only the n=0, ... ,nterms, m=0, ..., nterms coefficients
c          are calculated.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS,NPARTS
c
      REAL *8 X0Y0Z0(3),ZPARTS(3,NPARTS),CHARGE(NPARTS)
      REAL *8 P(0:NTERMS,0:NTERMS)
      REAL *8 C(0:NTERMS,0:NTERMS),SCALE,BETA
c
      COMPLEX *16 LOCAL(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 I,L,M,NCALC
      REAL *8 PRECIS,D,CP,PROJ
      REAL *8 RX,RY,RZ,RR,CTHETA
      REAL *8 RK,BK(0:60)
      COMPLEX *16 IMAG,EPHI(1:60)
c
      DATA IMAG/(0.0D0,1.0D0)/
      DATA PRECIS/1.0D-14/
c
      DO 4900 I = 1, NPARTS
        RX = ZPARTS(1,I) - X0Y0Z0(1)
        RY = ZPARTS(2,I) - X0Y0Z0(2)
        RZ = ZPARTS(3,I) - X0Y0Z0(3)
c
c-------compute  distance rr
c         ctheta = cos(theta)
c         ephi(1)=e^(i*phi)
c
        PROJ = RX*RX+RY*RY
        RR = PROJ+RZ*RZ
        PROJ = DSQRT(PROJ)
        D = DSQRT(RR)
c
c-------note: here is a hack. when computing cos(theta) as
c         rz/d, we have to be careful about the possibility
c         of d  being 0 (in which case ctheta is not well
c         defined - we arbitrarily set it to 1.)
c
        IF ( D .LE. PRECIS ) THEN
          CTHETA = 1.0D0
        ELSE
          CTHETA = RZ/D
        ENDIF
c
        IF ( PROJ .LE. PRECIS*D ) THEN
          EPHI(1) = 1.0D0
        ELSE
          EPHI(1) = RX/PROJ - IMAG*RY/PROJ
        ENDIF
c
c-------create array of powers of e^(-i*phi).
c
        DO 4100 L = 2,NTERMS
          EPHI(L) = EPHI(L-1)*EPHI(1)
4100    CONTINUE
c
c-------compute the sperical modified bessel function. note that bk is scaled.
c
        RK=D*BETA
        CALL KN(SCALE,RK,NTERMS+1,BK, NCALC)
        IF (NCALC .NE. NTERMS+1) THEN
          CALL PRIN2(' bk is, *', BK, NTERMS+1)
          CALL PRIN2(' rk is *', RK, 1)
          CALL PRIN2('scale is *', SCALE, 1)
          CALL PRINF(' the nterms+1 is *', NTERMS+1, 1)
          CALL PRINF(' the ncalc is *', NCALC, 1)
          PRINT *, 'error in forming local expansions.'
          STOP
        ENDIF
c
        LOCAL(0,0) = LOCAL(0,0) + CHARGE(I)*BK(0)
c
c-------compute legendre polynomials of argument cos(theta) = ctheta
c         and add contributions from legendre polynomials
c
        CALL LGNDR(NTERMS,CTHETA,P)
c
        DO 4300 L = 1,NTERMS
          CP = CHARGE(I)*P(L,0)*BK(L)*C(L,0)
          LOCAL(L,0) = LOCAL(L,0) + CP
4300    CONTINUE
c
c-------add contributions from associated legendre functions.
c
        DO 4500 L = 1,NTERMS
          DO 4400 M=1,L
            CP = CHARGE(I)*BK(L)*C(L,M)*P(L,M)
            LOCAL(L,M) = LOCAL(L,M) + CP*EPHI(M)
4400      CONTINUE
4500    CONTINUE
4900  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YBRMPEV(BETA,MPOLE,X0Y0Z0,POINT,NTERMS,RPOT,FIELD,
     1  SCALE,P)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    evaluates multipole expansion at arbitrary point.
c
c  on input:
c    beta : the frequency of the equation.
c    mpole: coefficients of multipole expansion(scaled)
c    x0y0z0: the center of the expansion
c    point: point of evaluation
c    nterms: order of expansion
c    p: work arrays to hold legendre polynomials
c       and associated legendre functions, respectively.
c
c  on output:
c    rpot: computed potential
c    field: computed field
c
c  note: the current version will compute the potential only,
c            and the computation of the field will be added later.
c
c************************************************************************
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 BETA,X0Y0Z0(3),POINT(3),RPOT
      REAL *8 P(0:NTERMS+1,0:NTERMS+1)
      REAL *8 FIELD(3)
      REAL *8 SCALE
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 NTP1,I,L,M,NCALC
      REAL *8 BK(0:60)
      REAL *8 PRECIS,D,CP,PROJ
      REAL *8 RX,RY,RZ,RR,CTHETA,RPOTZ
      COMPLEX *16 FIELDTEMP(3)
      COMPLEX *16 PTEMP(0:60,0:60)
      COMPLEX *16 IMAG,EPHI(1:60),CPZ
c
      DATA IMAG/(0.0D0,1.0D0)/
      DATA PRECIS/1.0D-14/
c
c-----functions called.
c
      REAL *8 DREAL
c
      NTP1 = NTERMS+1
c
c-----compute relevant functions of spherical coordinates
c       d = distance*beta, ctheta = cos(theta), ephi = exp(i*phi)
c
      RX = POINT(1) - X0Y0Z0(1)
      RY = POINT(2) - X0Y0Z0(2)
      RZ = POINT(3) - X0Y0Z0(3)
      PROJ = RX*RX+RY*RY
      RR = PROJ+RZ*RZ
      PROJ = DSQRT(PROJ)
      D = DSQRT(RR)
c
      IF (D .LE. PRECIS) THEN
        CTHETA = 0.0D0
      ELSE
        CTHETA = RZ/D
      ENDIF
      IF ( PROJ .LE. PRECIS*D ) THEN
        EPHI(1) = 1.0D0
      ELSE
        EPHI(1) = RX/PROJ + IMAG*RY/PROJ
      ENDIF
      D=D*BETA
c
c-----create array of powers of e^(i*m*phi).
c
      DO 5600 L = 1,NTERMS+2
        EPHI(L+1) = EPHI(L)*EPHI(1)
5600  CONTINUE
c
c-----compute the scaled modified bessel function.
c
      CALL KN(SCALE, D, NTERMS+1, BK, NCALC)
c
      IF (NCALC.NE.NTERMS+1) THEN
        PRINT *,'wrong output from besselk, stop',NTERMS+1,NCALC
        STOP
      ENDIF
c
c-----compute values of legendre functions
c
      CALL LGNDR(NTERMS+1,CTHETA,P)
c
c-----add contribution from 0 mode
c
      RPOT = RPOT + MPOLE(0,0)*BK(0)
      PTEMP(0,0)=BK(0)
c
c-----add contributions from legendre polynomials
c
      DO 5700 L = 1,NTERMS
        PTEMP(L,0)= BK(L)*P(L,0)
        RPOT = RPOT + MPOLE(L,0)*PTEMP(L,0)
5700  CONTINUE
c
c-----add contributions from associated legendre functions.
c
      DO 5900 L = 1,NTERMS
c
c-------compute potential
c
        DO 5800 M=1,L
          PTEMP(L,M)=(BK(L)*P(L,M))*EPHI(M)
          RPOT = RPOT + 2.0D0*DREAL(MPOLE(L,M)*PTEMP(L,M))
5800    CONTINUE
5900  CONTINUE
c
c-----now calculate the force. see ben et al jcp paper for
c       the formulas.
c
c-----construct a few more terms of ptemp for force field.
c
      DO M=0,NTERMS+1
        PTEMP(NTERMS+1,M)=(BK(NTERMS+1)*P(NTERMS+1,M))*EPHI(M)
      ENDDO
c
c-----contribution from (n=0,m=0) term
c
      RPOTZ=MPOLE(0,0)/SCALE
      CPZ=RPOTZ*PTEMP(1,1)
      FIELDTEMP(1)=-DCONJG(CPZ)
      FIELDTEMP(2)=-RPOTZ*PTEMP(1,0)
      FIELDTEMP(3)= CPZ
c
c-----contribution from (n=1,m=0)
c
      RPOTZ=MPOLE(1,0)/3.0D0
      CPZ=(RPOTZ/SCALE)*PTEMP(2,1)
      FIELDTEMP(1)=FIELDTEMP(1)-DCONJG(CPZ)
      FIELDTEMP(2)=FIELDTEMP(2)-RPOTZ*(PTEMP(0,0)*SCALE+
     1  2.0D0*PTEMP(2,0)/SCALE)
      FIELDTEMP(3)=FIELDTEMP(3)+CPZ
c
c-----contribution from (n>1,m=0)
c
      DO L=2,NTERMS
        RPOTZ=MPOLE(L,0)/DBLE(2*L+1)
        CPZ=RPOTZ*(PTEMP(L-1,1)*SCALE-PTEMP(L+1,1)/SCALE)
        FIELDTEMP(1)=FIELDTEMP(1)+DCONJG(CPZ)
        FIELDTEMP(2)=FIELDTEMP(2)-RPOTZ*(DREAL(PTEMP(L-1,0))*DBLE(L)*
     1    SCALE+DBLE(L+1)*DREAL(PTEMP(L+1,0))/SCALE)
        FIELDTEMP(3)=FIELDTEMP(3)-CPZ
      ENDDO
c
c-----now update the force field.
c
      FIELD(1)=FIELD(1)-BETA*DREAL(FIELDTEMP(1)-FIELDTEMP(3))/2.0D0
      FIELD(2)=FIELD(2)-BETA*DREAL((FIELDTEMP(1)+FIELDTEMP(3))*IMAG)
     1  /2.0D0
      FIELD(3)=FIELD(3)+BETA*DREAL(FIELDTEMP(2))
c
c-----contribution from (n=1,m=1)
c
      CPZ=MPOLE(1,1)/3.0D0
      FIELDTEMP(1)=-CPZ*(PTEMP(0,0)*SCALE-PTEMP(2,0)/SCALE)*2.0D0
      FIELDTEMP(2)=-CPZ*PTEMP(2,1)/SCALE
      FIELDTEMP(3)= CPZ*PTEMP(2,2)/SCALE
c
c-----contributions for n=2,nterms; m=1,n-2.
c
      DO L=2,NTERMS
        DO M=1,L-2
          CPZ=MPOLE(L,M)/DBLE(2*L+1)
          FIELDTEMP(1)=FIELDTEMP(1)-CPZ*(PTEMP(L-1,M-1)*
     1      (DBLE((L+M-1)*(L+M))*SCALE)-PTEMP(L+1,M-1)*
     2      (DBLE((L-M+1)*(L-M+2))/SCALE))
          FIELDTEMP(2)=FIELDTEMP(2)-CPZ*(PTEMP(L-1,M)*
     1      (DBLE(L+M)*SCALE)+PTEMP(L+1,M)*(DBLE(L-M+1)/SCALE))
          FIELDTEMP(3)=FIELDTEMP(3)-CPZ*(PTEMP(L-1,M+1)
     1      *SCALE-PTEMP(L+1,M+1)/SCALE)
        ENDDO
c
c-------m=n-1.
c
        CPZ=MPOLE(L,L-1)/DBLE(2*L+1)
        FIELDTEMP(1)=FIELDTEMP(1)-CPZ*(PTEMP(L-1,L-2)*
     1    (DBLE((L+L-2)*(L+L-1))*SCALE)-PTEMP(L+1,L-2)*
     2    (6.0D0/SCALE))
        FIELDTEMP(2)=FIELDTEMP(2)-CPZ*(PTEMP(L-1,L-1)*
     1    (DBLE(L+L-1)*SCALE)+PTEMP(L+1,L-1)*(2.0D0/SCALE))
        FIELDTEMP(3)=FIELDTEMP(3)+CPZ*PTEMP(L+1,L)/SCALE
c
c-------m=n.
c
        CPZ=MPOLE(L,L)/DBLE(2*L+1)
        FIELDTEMP(1)=FIELDTEMP(1)-CPZ*(PTEMP(L-1,L-1)*
     1    (DBLE((L+L-1)*(L+L))*SCALE)-PTEMP(L+1,L-1)*
     2    (2.0D0/SCALE))
        FIELDTEMP(2)=FIELDTEMP(2)-CPZ*PTEMP(L+1,L)/SCALE
        FIELDTEMP(3)=FIELDTEMP(3)+CPZ*PTEMP(L+1,L+1)/SCALE
c
      ENDDO
c
c-----now calculate the force.
c
      FIELD(1)=FIELD(1)-BETA*DREAL(FIELDTEMP(1)-FIELDTEMP(3))
      FIELD(2)=FIELD(2)-BETA*DREAL((FIELDTEMP(1)+FIELDTEMP(3))*IMAG)
      FIELD(3)=FIELD(3)+2.0D0*BETA*DREAL(FIELDTEMP(2))
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YEXPTOLOCAL(BETA,RLSC,LOCAL,NTERMS,RLAMS,WHTS,NLAMBS,
     1  YTOP,NUMTETS,NTHMAX,NEXPTOT,MEXPUP,MEXPDOWN,SCALE)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    this subroutine converts the fourier representation of two
c    exponential moment functions into a local multipole expansion
c    (with respect to the same box center).
c      l_n^m= (see reference). and is scaled.
c
c    u(x,y,z) = \int_0^\infty e^{-\lambda z}
c                \int_0^{2\pi} e^{i\lambda(xcos(alpha)+ysin(alpha))}
c                mexpup(lambda,alpha) dalpha dlambda
c            +
c                \int_0^\infty e^{\lambda z}
c                \int_0^{2\pi} e^{i\lambda(xcos(alpha)+ysin(alpha))}
c                mexpdown(lambda,alpha) dalpha dlambda
c
c             = \sum_{n=0}^{nterms} \sum_{m=-n,n}
c                local(n,m) y_n^m(cos theta) e^{i m \phi} r^{n}
c
c  on input:
c    beta : the scaled beta used in the calculation of l_n^m.
c    rlsc : the precomputed and scaled p_n^m *sc^n
c    nterms : the total number of expansions.
c    mexpup(nexptot): fourier coefficients of the function
c                    mexpup for discrete lambda
c                    values. they are ordered as follows:
c
c                 mexpup(1,...,numtets(1)) = fourier modes
c                             for lambda_1
c                 mexpup(numtets(1)+1,...,numtets(2)) = fourier modes
c                             for lambda_2
c                 etc.
c    mexpdown(nexptot): as above for down expansion
c    rlams(nlambs): discretization points in lambda integral
c    whts(nlambs): quadrature weights in lambda integral
c    nlambs:      number of discretization pts. in lambda integral
c    numtets(j): number of fourier modes in expansion of alpha
c                variable for lambda_j.
c    nthmax:     max_j numtets(j)
c    nexptot:    sum_j numtets(j)
c
c  on output:
c    local(0:nterms,0:nterms): output multipole expansion of order
c                              nterms.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS,NLAMBS,NUMTETS(NLAMBS),NEXPTOT,NTHMAX
      INTEGER *4 IEXPU,IEXPD
c
      REAL *8 BETA,RLSC(0:NTERMS,0:NTERMS,NLAMBS)
      REAL *8 RLAMS(NLAMBS),WHTS(NLAMBS),SCALE
      REAL *8 YTOP(0:NTERMS,0:NTERMS)
c
      COMPLEX *16 LOCAL(0:NTERMS,0:NTERMS)
      COMPLEX *16 MEXPUP(NEXPTOT)
      COMPLEX *16 MEXPDOWN(NEXPTOT)
c
c-----local variables.
c
      INTEGER *4 I,NM,NTOT,NL,MMAX,MTH,NCURRENT
      REAL *8 PI,RSCALE
      COMPLEX *16 ZEYE(0:200)
      COMPLEX *16 MEXPPLUS(2000)
      COMPLEX *16 MEXPMINUS(2000)
c
c-----local functions.
c
      REAL *8 DATAN
      COMPLEX *16 DCMPLX
c
      PI=DATAN(1.0D0)*4.0D0
c
c-----compute necessary powers of -i
c
      ZEYE(0) = 1.0D0
      DO I = 1,NTHMAX
        ZEYE(I) = ZEYE(I-1)*DCMPLX(0.0D0,1.0D0)
      ENDDO
c
c-----initialize local expansion
c
      DO NM = 0,NTERMS
        DO MTH = 0,NTERMS
          LOCAL(NM,MTH) = 0.0D0
        ENDDO
      ENDDO
c
c-----compute sum and difference of mexpup and mexpdown
c
      DO NM = 1,NEXPTOT
        MEXPPLUS(NM) = MEXPDOWN(NM) + MEXPUP(NM)
        MEXPMINUS(NM) = MEXPDOWN(NM) - MEXPUP(NM)
      ENDDO
c
c-----loop over multipole order to generate mexp values.
c
      NTOT = 1
      DO NL = 1,NLAMBS
c
c-------add contributions to local expansion. first compute
c       p_n^m*w_k*mexplus/minus.
c
        DO NM = 0,NTERMS,2
          MMAX = NUMTETS(NL)-1
          IF (MMAX.GT.NM) MMAX = NM
          DO MTH = 0,MMAX
            NCURRENT = NTOT+MTH
            LOCAL(NM,MTH) = LOCAL(NM,MTH) + RLSC(NM,MTH,NL)*
     1        WHTS(NL)*MEXPPLUS(NCURRENT)
          ENDDO
        ENDDO
c
        DO NM = 1,NTERMS,2
          MMAX = NUMTETS(NL)-1
          IF (MMAX.GT.NM) MMAX = NM
          DO MTH = 0,MMAX
            NCURRENT = NTOT+MTH
            LOCAL(NM,MTH) = LOCAL(NM,MTH) + RLSC(NM,MTH,NL)*
     1        WHTS(NL)*MEXPMINUS(NCURRENT)
          ENDDO
        ENDDO
        NTOT = NTOT+NUMTETS(NL)
      ENDDO
c
c-----scale the expansions according to formula
c
      RSCALE=PI/BETA/2.0D0
c
      DO NM = 0,NTERMS
        DO MTH = 0,NM
          LOCAL(NM,MTH)=LOCAL(NM,MTH)*ZEYE(MTH)*RSCALE*YTOP(NM,MTH)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPROCESSUP(SCALE,LEXP1,IUALL,NUALL,IXUALL,
     1  IYUALL,IU1234,NU1234,IX1234,IY1234,MEXUALL,
     2  MEXU1234,XS,YS,ZS,NEXPTOTP,MNEXPTOTP)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine processes the up interaction lists.
c
c  on input:
c
c     iuall(nuall), ixuall(nuall), iyuall(nuall) are the boxes
c          receiving all child box data and the x and y offsets from
c          child 1, respectively.
c     iu1234(nu1234), ix1234(nu1234), iy1234(nu1234) are the boxes
c          receiving data from child boxes 1,2,3,4 and the x and y
c          offsets from child 1, respectively.
c
c     mexuall is the exponential expansion for all eight children.
c     mexu1234(nexptotp) is the exponential expansion for
c          children 1,2,3,4.
c     xs,ys,zs are the shift coefficients computed by subroutine
c          mkexps.
c
c  on output:
c
c     lexp1, which contains the local up expansion information for
c            all boxes, is incremented for each box in the
c            interaction lists.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NEXPTOTP,MNEXPTOTP
      INTEGER *4 NUALL, NU1234
      INTEGER *4 IUALL(NUALL),IXUALL(NUALL),IYUALL(NUALL)
      INTEGER *4 IU1234(NU1234),IX1234(NU1234),IY1234(NU1234)
c
      REAL *8 SCALE
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 LEXP1(MNEXPTOTP,1)
      COMPLEX *16 MEXUALL(NEXPTOTP)
      COMPLEX *16 MEXU1234(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 I,JJ
      COMPLEX *16 ZMUL
c
      DO I = 1,NUALL
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(3,JJ)
          IF (IXUALL(I).GT.0)
     1      ZMUL = ZMUL*XS(IXUALL(I),JJ)
          IF (IXUALL(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IXUALL(I),JJ))
          IF (IYUALL(I).GT.0)
     1      ZMUL = ZMUL*YS(IYUALL(I),JJ)
          IF (IYUALL(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IYUALL(I),JJ))
          LEXP1(JJ,IUALL(I)) = LEXP1(JJ,IUALL(I)) +
     1      MEXUALL(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NU1234
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX1234(I).GT.0)
     1      ZMUL = ZMUL*XS(IX1234(I),JJ)
          IF (IX1234(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX1234(I),JJ))
          IF (IY1234(I).GT.0)
     1      ZMUL = ZMUL*YS(IY1234(I),JJ)
          IF (IY1234(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY1234(I),JJ))
          LEXP1(JJ,IU1234(I)) = LEXP1(JJ,IU1234(I)) +
     1      MEXU1234(JJ)*ZMUL
        ENDDO
      ENDDO
c
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPROCESSDN(SCALE,LEXP2,IDALL,NDALL,IXDALL,
     1  IYDALL,ID5678,ND5678,IX5678,IY5678,MEXDALL,
     2  MEXD5678,XS,YS,ZS,NEXPTOTP,MNEXPTOTP)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine processes the down interaction lists.
c
c  on input:
c
c     idall(ndall), ixdall(ndall), iydall(ndall) are the boxes
c          receiving all child box data and the x and y offsets from
c          child 1, respectively.
c     id5678(nd5678), ix5678(nd5678), iy5678(nd5678) are the boxes
c          receiving data from child boxes 5,6,7,8 and the x and y
c          offsets from child 1, respectively.
c
c     mexdall is the exponential expansion for all eight children.
c     mexd5678(nexptotp) is the exponential expansion for
c          children 5,6,7,8.
c     xs,ys,zs are the shift coefficients computed by subroutine
c          mkexps.
c
c  on output:
c
c     lexp2: which contains the local down expansion information for
c            all boxes, is incremented for each box in the
c            interaction lists.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NEXPTOTP,MNEXPTOTP
      INTEGER *4 NDALL,ND5678
      INTEGER *4 IDALL(NDALL),IXDALL(NDALL),IYDALL(NDALL)
      INTEGER *4 ID5678(ND5678),IX5678(ND5678),IY5678(ND5678)
c
      REAL *8 SCALE
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 LEXP2(MNEXPTOTP,1)
      COMPLEX *16 MEXDALL(NEXPTOTP)
      COMPLEX *16 MEXD5678(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 I, JJ
      COMPLEX *16 ZMUL
c
      DO I = 1,NDALL
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IXDALL(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IXDALL(I),JJ))
          IF (IXDALL(I).LT.0)
     1      ZMUL = ZMUL*XS(-IXDALL(I),JJ)
          IF (IYDALL(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IYDALL(I),JJ))
          IF (IYDALL(I).LT.0)
     1      ZMUL = ZMUL*YS(-IYDALL(I),JJ)
          LEXP2(JJ,IDALL(I)) = LEXP2(JJ,IDALL(I)) +
     1      MEXDALL(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,ND5678
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX5678(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX5678(I),JJ))
          IF (IX5678(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX5678(I),JJ)
          IF (IY5678(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY5678(I),JJ))
          IF (IY5678(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY5678(I),JJ)
          LEXP2(JJ,ID5678(I)) = LEXP2(JJ,ID5678(I)) +
     1      MEXD5678(JJ)*ZMUL
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPROCESSNO(SCALE,LEXP1,INALL,NNALL,IXNALL,IYNALL,
     1  IN1256,NN1256,IX1256,IY1256,
     2  IN12,NN12,IX12,IY12,IN56,NN56,IX56,IY56,
     3  MEXNALL,MEXN1256,MEXN12,MEXN56,XS,YS,ZS,NEXPTOTP,
     4  MNEXPTOTP)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine processes the north interaction lists.
c
c  on input:
c
c     inall(nnall), ixnall(nnall), iynall(nnall) are the boxes
c          receiving all child box data and the x and y offsets from
c          child 1, respectively.
c     the other north lists are similarly defined (see ymknolist).
c
c     mexnall is the exponential expansion for all eight children, etc.
c     xs,ys,zs are the shift coefficients computed by subroutine
c          mkexps.
c
c  on output:
c
c     lexp1, which contains the local north expansion information for
c            all boxes, is incremented for each box in the
c            interaction lists.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NEXPTOTP,MNEXPTOTP
      INTEGER *4 NNALL,NN1256,NN12,NN56
      INTEGER *4 INALL(NNALL),IXNALL(NNALL),IYNALL(NNALL)
      INTEGER *4 IN1256(NN1256),IX1256(NN1256),IY1256(NN1256)
      INTEGER *4 IN12(NN12),IX12(NN12),IY12(NN12)
      INTEGER *4 IN56(NN56),IX56(NN56),IY56(NN56)
c
      REAL *8 SCALE
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 LEXP1(MNEXPTOTP,*)
      COMPLEX *16 MEXNALL(NEXPTOTP)
      COMPLEX *16 MEXN1256(NEXPTOTP)
      COMPLEX *16 MEXN12(NEXPTOTP)
      COMPLEX *16 MEXN56(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 I, JJ
      COMPLEX *16 ZMUL
c
      DO I = 1,NNALL
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(3,JJ)
          IF (IXNALL(I).GT.0)
     1      ZMUL = ZMUL*XS(IXNALL(I),JJ)
          IF (IXNALL(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IXNALL(I),JJ))
          IF (IYNALL(I).GT.0)
     1      ZMUL = ZMUL*YS(IYNALL(I),JJ)
          IF (IYNALL(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IYNALL(I),JJ))
          LEXP1(JJ,INALL(I)) = LEXP1(JJ,INALL(I)) +
     1      MEXNALL(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NN1256
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX1256(I).GT.0)
     1      ZMUL = ZMUL*XS(IX1256(I),JJ)
          IF (IX1256(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX1256(I),JJ))
          IF (IY1256(I).GT.0)
     1      ZMUL = ZMUL*YS(IY1256(I),JJ)
          IF (IY1256(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY1256(I),JJ))
          LEXP1(JJ,IN1256(I)) = LEXP1(JJ,IN1256(I)) +
     1      MEXN1256(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NN12
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX12(I).GT.0)
     1      ZMUL = ZMUL*XS(IX12(I),JJ)
          IF (IX12(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX12(I),JJ))
          IF (IY12(I).GT.0)
     1      ZMUL = ZMUL*YS(IY12(I),JJ)
          IF (IY12(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY12(I),JJ))
          LEXP1(JJ,IN12(I)) = LEXP1(JJ,IN12(I)) +
     1      MEXN12(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NN56
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX56(I).GT.0)
     1      ZMUL = ZMUL*XS(IX56(I),JJ)
          IF (IX56(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX56(I),JJ))
          IF (IY56(I).GT.0)
     1      ZMUL = ZMUL*YS(IY56(I),JJ)
          IF (IY56(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY56(I),JJ))
          LEXP1(JJ,IN56(I)) = LEXP1(JJ,IN56(I)) +
     1      MEXN56(JJ)*ZMUL
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPROCESSSO(SCALE,LEXP2,ISALL,NSALL,IXSALL,IYSALL,
     1  IS3478,NS3478,IX3478,IY3478,
     1  IS34,NS34,IX34,IY34,IS78,NS78,IX78,IY78,
     2  MEXSALL,MEXS3478,MEXS34,MEXS78,XS,YS,ZS,NEXPTOTP,
     3  MNEXPTOTP)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine processes the south interaction lists.
c
c  on input:
c
c     isall(nsall), ixsall(nsall), iysall(nsall) are the boxes
c          receiving all child box data and the x and y offsets from
c          child 1, respectively.
c     the other south lists are similarly defined (see mksolist).
c
c     mexsall is the exponential expansion for all eight children, etc.
c     xs,ys,zs are the shift coefficients computed by subroutine
c          mkexps.
c
c  on output:
c
c     lexp2, which contains the local north expansion information for
c            all boxes, is incremented for each box in the
c            interaction lists.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NEXPTOTP,MNEXPTOTP
      INTEGER *4 NSALL,NS3478,NS34,NS78
      INTEGER *4 ISALL(NSALL),IXSALL(NSALL),IYSALL(NSALL)
      INTEGER *4 IS3478(NS3478),IX3478(NS3478),IY3478(NS3478)
      INTEGER *4 IS34(NS34),IX34(NS34),IY34(NS34)
      INTEGER *4 IS78(NS78),IX78(NS78),IY78(NS78)
      INTEGER *4 IEXP(8),IEXP1(1)
c
      REAL *8 SCALE
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 LEXP2(MNEXPTOTP,*)
      COMPLEX *16 MEXSALL(NEXPTOTP)
      COMPLEX *16 MEXS3478(NEXPTOTP)
      COMPLEX *16 MEXS34(NEXPTOTP)
      COMPLEX *16 MEXS78(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 I,JJ
      COMPLEX *16 ZMUL
c
      DO I = 1,NSALL
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IXSALL(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IXSALL(I),JJ))
          IF (IXSALL(I).LT.0)
     1      ZMUL = ZMUL*XS(-IXSALL(I),JJ)
          IF (IYSALL(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IYSALL(I),JJ))
          IF (IYSALL(I).LT.0)
     1      ZMUL = ZMUL*YS(-IYSALL(I),JJ)
          LEXP2(JJ,ISALL(I)) = LEXP2(JJ,ISALL(I)) +
     1      MEXSALL(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NS3478
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX3478(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX3478(I),JJ))
          IF (IX3478(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX3478(I),JJ)
          IF (IY3478(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY3478(I),JJ))
          IF (IY3478(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY3478(I),JJ)
          LEXP2(JJ,IS3478(I)) = LEXP2(JJ,IS3478(I)) +
     1      MEXS3478(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NS34
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX34(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX34(I),JJ))
          IF (IX34(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX34(I),JJ)
          IF (IY34(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY34(I),JJ))
          IF (IY34(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY34(I),JJ)
          LEXP2(JJ,IS34(I)) = LEXP2(JJ,IS34(I)) +
     1      MEXS34(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NS78
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX78(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX78(I),JJ))
          IF (IX78(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX78(I),JJ)
          IF (IY78(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY78(I),JJ))
          IF (IY78(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY78(I),JJ)
          LEXP2(JJ,IS78(I)) = LEXP2(JJ,IS78(I)) +
     1      MEXS78(JJ)*ZMUL
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPROCESSEA(SCALE,LEXP1,IEALL,NEALL,IXEALL,IYEALL,
     1  IE1357,NE1357,IX1357,IY1357,IE13,NE13,IX13,IY13,IE57,NE57,
     2  IX57,IY57,IE1,NE1,IX1,IY1,IE3,NE3,IX3,IY3,IE5,NE5,IX5,IY5,
     3  IE7,NE7,IX7,IY7,MEXEALL,MEXE1357,MEXE13,MEXE57,
     4  MEXE1,MEXE3,MEXE5,MEXE7,XS,YS,ZS,NEXPTOTP,MNEXPTOTP)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine processes the east interaction lists.
c
c  on input:
c
c     ieall(neall), ixeall(neall), iyeall(neall) are the boxes
c          receiving all child box data and the x and y offsets from
c          child 1, respectively.
c     the other east lists are similarly defined (see mkealist).
c
c     mexeall is the exponential expansion for all eight children, etc.
c     xs,ys,zs are the shift coefficients computed by subroutine
c          mkexps.
c
c  on output:
c
c     lexp1, which contains the local east expansion information for
c            all boxes, is incremented for each box in the
c            interaction lists.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NEXPTOTP,MNEXPTOTP
      INTEGER *4 NEALL,NE1357,NE13,NE57,NE1,NE3,NE5,NE7
      INTEGER *4 IEALL(NEALL),IXEALL(NEALL),IYEALL(NEALL)
      INTEGER *4 IE1357(NE1357),IX1357(NE1357),IY1357(NE1357)
      INTEGER *4 IE13(NE13),IX13(NE13),IY13(NE13)
      INTEGER *4 IE57(NE57),IX57(NE57),IY57(NE57)
      INTEGER *4 IE1(NE1),IX1(NE1),IY1(NE1)
      INTEGER *4 IE3(NE3),IX3(NE3),IY3(NE3)
      INTEGER *4 IE5(NE5),IX5(NE5),IY5(NE5)
      INTEGER *4 IE7(NE7),IX7(NE7),IY7(NE7)
c
      REAL *8 SCALE
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 LEXP1(MNEXPTOTP,*)
      COMPLEX *16 MEXEALL(NEXPTOTP)
      COMPLEX *16 MEXE1357(NEXPTOTP)
      COMPLEX *16 MEXE13(NEXPTOTP)
      COMPLEX *16 MEXE57(NEXPTOTP)
      COMPLEX *16 MEXE1(NEXPTOTP)
      COMPLEX *16 MEXE3(NEXPTOTP)
      COMPLEX *16 MEXE5(NEXPTOTP)
      COMPLEX *16 MEXE7(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 I, JJ
      COMPLEX *16 ZMUL
c
      DO I = 1,NEALL
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(3,JJ)
          IF (IXEALL(I).GT.0)
     1      ZMUL = ZMUL*XS(IXEALL(I),JJ)
          IF (IXEALL(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IXEALL(I),JJ))
          IF (IYEALL(I).GT.0)
     1      ZMUL = ZMUL*YS(IYEALL(I),JJ)
          IF (IYEALL(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IYEALL(I),JJ))
          LEXP1(JJ,IEALL(I)) = LEXP1(JJ,IEALL(I)) +
     1      MEXEALL(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE1357
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX1357(I).GT.0)
     1      ZMUL = ZMUL*XS(IX1357(I),JJ)
          IF (IX1357(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX1357(I),JJ))
          IF (IY1357(I).GT.0)
     1      ZMUL = ZMUL*YS(IY1357(I),JJ)
          IF (IY1357(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY1357(I),JJ))
          LEXP1(JJ,IE1357(I)) = LEXP1(JJ,IE1357(I)) +
     1      MEXE1357(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE13
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX13(I).GT.0)
     1      ZMUL = ZMUL*XS(IX13(I),JJ)
          IF (IX13(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX13(I),JJ))
          IF (IY13(I).GT.0)
     1      ZMUL = ZMUL*YS(IY13(I),JJ)
          IF (IY13(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY13(I),JJ))
          LEXP1(JJ,IE13(I)) = LEXP1(JJ,IE13(I)) +
     1      MEXE13(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE57
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX57(I).GT.0)
     1      ZMUL = ZMUL*XS(IX57(I),JJ)
          IF (IX57(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX57(I),JJ))
          IF (IY57(I).GT.0)
     1      ZMUL = ZMUL*YS(IY57(I),JJ)
          IF (IY57(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY57(I),JJ))
          LEXP1(JJ,IE57(I)) = LEXP1(JJ,IE57(I)) +
     1      MEXE57(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE1
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX1(I).GT.0)
     1      ZMUL = ZMUL*XS(IX1(I),JJ)
          IF (IX1(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX1(I),JJ))
          IF (IY1(I).GT.0)
     1      ZMUL = ZMUL*YS(IY1(I),JJ)
          IF (IY1(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY1(I),JJ))
          LEXP1(JJ,IE1(I)) = LEXP1(JJ,IE1(I)) +
     1      MEXE1(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE3
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX3(I).GT.0)
     1      ZMUL = ZMUL*XS(IX3(I),JJ)
          IF (IX3(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX3(I),JJ))
          IF (IY3(I).GT.0)
     1      ZMUL = ZMUL*YS(IY3(I),JJ)
          IF (IY3(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY3(I),JJ))
          LEXP1(JJ,IE3(I)) = LEXP1(JJ,IE3(I)) +
     1      MEXE3(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE5
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX5(I).GT.0)
     1      ZMUL = ZMUL*XS(IX5(I),JJ)
          IF (IX5(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX5(I),JJ))
          IF (IY5(I).GT.0)
     1      ZMUL = ZMUL*YS(IY5(I),JJ)
          IF (IY5(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY5(I),JJ))
          LEXP1(JJ,IE5(I)) = LEXP1(JJ,IE5(I)) +
     1      MEXE5(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NE7
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IX7(I).GT.0)
     1      ZMUL = ZMUL*XS(IX7(I),JJ)
          IF (IX7(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(XS(-IX7(I),JJ))
          IF (IY7(I).GT.0)
     1      ZMUL = ZMUL*YS(IY7(I),JJ)
          IF (IY7(I).LT.0)
     1      ZMUL = ZMUL*DCONJG(YS(-IY7(I),JJ))
          LEXP1(JJ,IE7(I)) = LEXP1(JJ,IE7(I)) +
     1      MEXE7(JJ)*ZMUL
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YPROCESSWE(SCALE,LEXP2,IWALL,NWALL,IXWALL,IYWALL,
     1  IW2468,NW2468,IX2468,IY2468,IW24,NW24,IX24,IY24,IW68,NW68,
     2  IX68,IY68,IW2,NW2,IX2,IY2,IW4,NW4,IX4,IY4,IW6,NW6,IX6,IY6,
     3  IW8,NW8,IX8,IY8,MEXWALL,MEXW2468,MEXW24,MEXW68,MEXW2,MEXW4,
     4  MEXW6,MEXW8,XS,YS,ZS,NEXPTOTP,MNEXPTOTP)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine processes the west interaction lists.
c
c  on input:
c
c     iwall(nwall), ixwall(nwall), iywall(nwall) are the boxes
c          receiving all child box data and the x and y offsets from
c          child 1, respectively.
c     the other west lists are similarly defined (see mkwelist).
c
c     mexeall is the exponential expansion for all eight children, etc.
c     xs,ys,zs are the shift coefficients computed by subroutine
c          mkexps.
c
c  on output:
c
c     lexp1, which contains the local west expansion information for
c            all boxes, is incremented for each box in the
c            interaction lists.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NEXPTOTP,MNEXPTOTP,NWALL,NW2468,NW24,NW68,NW2
      INTEGER *4 NW4,NW6,NW8
      INTEGER *4 IWALL(NWALL),IXWALL(NWALL),IYWALL(NWALL)
      INTEGER *4 IW2468(NW2468),IX2468(NW2468),IY2468(NW2468)
      INTEGER *4 IW24(NW24),IX24(NW24),IY24(NW24)
      INTEGER *4 IW68(NW68),IX68(NW68),IY68(NW68)
      INTEGER *4 IW2(NW2),IX2(NW2),IY2(NW2)
      INTEGER *4 IW4(NW4),IX4(NW4),IY4(NW4)
      INTEGER *4 IW6(NW6),IX6(NW6),IY6(NW6)
      INTEGER *4 IW8(NW8),IX8(NW8),IY8(NW8)
      INTEGER *4 IEXP(16),IEXP1(1)
c
      REAL *8 SCALE
      REAL *8 ZS(3,NEXPTOTP)
c
      COMPLEX *16 LEXP2(MNEXPTOTP,*)
      COMPLEX *16 MEXWALL(NEXPTOTP)
      COMPLEX *16 MEXW2468(NEXPTOTP)
      COMPLEX *16 MEXW24(NEXPTOTP)
      COMPLEX *16 MEXW68(NEXPTOTP)
      COMPLEX *16 MEXW2(NEXPTOTP)
      COMPLEX *16 MEXW4(NEXPTOTP)
      COMPLEX *16 MEXW6(NEXPTOTP)
      COMPLEX *16 MEXW8(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
c
c-----local variables.
c
      INTEGER *4 I,JJ
      COMPLEX *16 ZMUL
c
      DO I = 1,NWALL
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(2,JJ)
          IF (IXWALL(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IXWALL(I),JJ))
          IF (IXWALL(I).LT.0)
     1      ZMUL = ZMUL*XS(-IXWALL(I),JJ)
          IF (IYWALL(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IYWALL(I),JJ))
          IF (IYWALL(I).LT.0)
     1      ZMUL = ZMUL*YS(-IYWALL(I),JJ)
          LEXP2(JJ,IWALL(I)) = LEXP2(JJ,IWALL(I)) +
     1      MEXWALL(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW2468
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX2468(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX2468(I),JJ))
          IF (IX2468(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX2468(I),JJ)
          IF (IY2468(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY2468(I),JJ))
          IF (IY2468(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY2468(I),JJ)
          LEXP2(JJ,IW2468(I)) = LEXP2(JJ,IW2468(I)) +
     1      MEXW2468(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW24
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX24(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX24(I),JJ))
          IF (IX24(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX24(I),JJ)
          IF (IY24(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY24(I),JJ))
          IF (IY24(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY24(I),JJ)
          LEXP2(JJ,IW24(I)) = LEXP2(JJ,IW24(I)) +
     1      MEXW24(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW68
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX68(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX68(I),JJ))
          IF (IX68(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX68(I),JJ)
          IF (IY68(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY68(I),JJ))
          IF (IY68(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY68(I),JJ)
          LEXP2(JJ,IW68(I)) = LEXP2(JJ,IW68(I)) +
     1      MEXW68(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW2
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX2(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX2(I),JJ))
          IF (IX2(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX2(I),JJ)
          IF (IY2(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY2(I),JJ))
          IF (IY2(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY2(I),JJ)
          LEXP2(JJ,IW2(I)) = LEXP2(JJ,IW2(I)) +
     1      MEXW2(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW4
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX4(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX4(I),JJ))
          IF (IX4(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX4(I),JJ)
          IF (IY4(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY4(I),JJ))
          IF (IY4(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY4(I),JJ)
          LEXP2(JJ,IW4(I)) = LEXP2(JJ,IW4(I)) +
     1      MEXW4(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW6
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX6(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX6(I),JJ))
          IF (IX6(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX6(I),JJ)
          IF (IY6(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY6(I),JJ))
          IF (IY6(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY6(I),JJ)
          LEXP2(JJ,IW6(I)) = LEXP2(JJ,IW6(I)) +
     1      MEXW6(JJ)*ZMUL
        ENDDO
      ENDDO
c
      DO I = 1,NW8
        DO JJ = 1,NEXPTOTP
          ZMUL = ZS(1,JJ)
          IF (IX8(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(XS(IX8(I),JJ))
          IF (IX8(I).LT.0)
     1      ZMUL = ZMUL*XS(-IX8(I),JJ)
          IF (IY8(I).GT.0)
     1      ZMUL = ZMUL*DCONJG(YS(IY8(I),JJ))
          IF (IY8(I).LT.0)
     1      ZMUL = ZMUL*YS(-IY8(I),JJ)
          LEXP2(JJ,IW8(I)) = LEXP2(JJ,IW8(I)) +
     1      MEXW8(JJ)*ZMUL
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c     this following file contains all of the expansion creation
c     routines for a parent box from its four children.
c
c     mkudexp creates all up and down expansions centered at child 1
c
c     mknsexp creates all north and south expansions centered at child 1
c
c     mkewexp creates all east and west expansions centered at child 1
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMKUDEXP(IBOX,LADDER,NLEV,NTERMS,MPOLE,RLAMS,NLAMBS,
     1           NUMFOUR,NUMPHYS,NTHMAX,NEXPTOT,NEXPTOTP,MEXPUP,
     2           MEXPDOWN,MEXPUPHYS,MEXPDPHYS,MEXUALL,MEXU1234,
     3           MEXDALL,MEXD5678,XS,YS,ZS,
     4           FEXPE,FEXPO,RLSC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine creates the up (+z)  and down (-z) exponential
c     expansions for a parent box due to all eight children, if exist.
c
c  note:
c
c     some intelligence is used in the order of summation. thus
c     mexu1234 is computed first and used to initialize mexuall.
c     the contributions from boxes 5 6 7 8  are then added in
c     separately, etc.
c
c  on input:
c
c     ibox: current box number.
c     box: current box information.
c     nterms: number of terms in the multipole expansion.
c     mpole: the multipole expansion coefficients.
c     rlams: exponential expansion coefficients.
c     nlambs: number of terms in the exponential expansion.
c     numfour: number of fourier modes in the expansion.
c     numphys: number of modes in the plane wave expansion.
c     nthmax: max number of terms in the exponential expansion.
c     nexptot: total number of fourier modes in the expansion.
c     nextpotp: half of the fourier modes.
c
c  precomputed tables:
c
c    xs,ys,zs: stores the diagonal translation operators when shifting exponential
c        expansions.
c    fexpe,fexp0: how exponential expansions will be merged.
c    rlsc: stores p_n^m for different lambda_k.
c
c  on output:
c
c    mexuall: up expansion from all boxes.
c    mexu1234: up expansion from box 1-4.
c    mexdall: down expansion from all box.
c    mexd5678: down expansion from 5-8.
c
c  variables:
c    mexpup:
c    mexpdown
c    mexpuphys:
c    mexpdphys:
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 LADDER(*),IBOX,NLEV
      INTEGER *4 NTERMS
      INTEGER *4 NLAMBS,NUMFOUR(NLAMBS),NEXPTOT,NTHMAX
      INTEGER *4 NUMPHYS(NLAMBS),NEXPTOTP
c
      REAL *8 ZS(3,NEXPTOTP)
      REAL *8 RLAMS(NLAMBS)
      REAL *8 RLSC(0:NTERMS,0:NTERMS,NLAMBS)
c
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS,*)
      COMPLEX *16 MEXPUP(NEXPTOT),MEXPDOWN(NEXPTOT)
      COMPLEX *16 MEXPUPHYS(NEXPTOTP),MEXPDPHYS(NEXPTOTP)
      COMPLEX *16 MEXUALL(NEXPTOTP),MEXU1234(NEXPTOTP)
      COMPLEX *16 MEXDALL(NEXPTOTP),MEXD5678(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP)
      COMPLEX *16 YS(3,NEXPTOTP)
      COMPLEX *16 FEXPE(1),FEXPO(1)
c
c-----local varialbles:
c
      INTEGER *4 JJ,ISTART,MYBOX
      INTEGER *4 ICHILD(8)
      COMPLEX *16 ZTMP
c
      ISTART = LADDER(NLEV)
      MYBOX = ISTART + IBOX
      CALL MKCHILD(NLEV,LADDER,MYBOX,ICHILD)
c
c-----add contributions from child 1
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(1)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXU1234(JJ) = MEXPUPHYS(JJ)
        MEXDALL(JJ) =  MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 2
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(2)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXPUPHYS(JJ) = MEXPUPHYS(JJ)*DCONJG(XS(1,JJ))
        MEXU1234(JJ) = MEXU1234(JJ) + MEXPUPHYS(JJ)
        MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*XS(1,JJ)
        MEXDALL(JJ) = MEXDALL(JJ) + MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 3
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(3)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
         MEXPUPHYS(JJ) = MEXPUPHYS(JJ)*DCONJG(YS(1,JJ))
         MEXU1234(JJ) = MEXU1234(JJ) + MEXPUPHYS(JJ)
         MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*YS(1,JJ)
         MEXDALL(JJ) = MEXDALL(JJ) + MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 4
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(4)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = XS(1,JJ)*YS(1,JJ)
        MEXPUPHYS(JJ) = MEXPUPHYS(JJ)*DCONJG(ZTMP)
        MEXU1234(JJ) = MEXU1234(JJ) + MEXPUPHYS(JJ)
        MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*ZTMP
        MEXDALL(JJ) = MEXDALL(JJ) + MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 5
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(5)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXPUPHYS(JJ) = MEXPUPHYS(JJ)/ZS(1,JJ)
        MEXUALL(JJ) = MEXU1234(JJ) + MEXPUPHYS(JJ)
        MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*ZS(1,JJ)
        MEXD5678(JJ) = MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 6
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(6)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = XS(1,JJ)*ZS(1,JJ)
        MEXPUPHYS(JJ) = MEXPUPHYS(JJ)/ZTMP
        MEXUALL(JJ) = MEXUALL(JJ) + MEXPUPHYS(JJ)
        MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*ZTMP
        MEXD5678(JJ) = MEXD5678(JJ) + MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 7
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(7)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = YS(1,JJ)*ZS(1,JJ)
        MEXPUPHYS(JJ) = MEXPUPHYS(JJ)/ZTMP
        MEXUALL(JJ) = MEXUALL(JJ) + MEXPUPHYS(JJ)
        MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*ZTMP
        MEXD5678(JJ) = MEXD5678(JJ) + MEXPDPHYS(JJ)
      ENDDO
c
c-----add contributions from child 8
c
      CALL YMPOLETOEXP(MPOLE(0,0,ICHILD(8)),NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPUP,MEXPDOWN,RLSC)
      CALL YFTOPHYS(MEXPUP,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPUPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPDOWN,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPDPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = XS(1,JJ)*YS(1,JJ)*ZS(1,JJ)
        MEXPUPHYS(JJ) = MEXPUPHYS(JJ)/ZTMP
        MEXUALL(JJ) = MEXUALL(JJ) + MEXPUPHYS(JJ)
        MEXPDPHYS(JJ) = MEXPDPHYS(JJ)*ZTMP
        MEXD5678(JJ) = MEXD5678(JJ) + MEXPDPHYS(JJ)
        MEXDALL(JJ) = MEXDALL(JJ) + MEXD5678(JJ)
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMKNSEXP(IBOX,LADDER,NLEV,NTERMS,MPOLE,MROTATE,MWORK,
     1  RLAMS,NLAMBS,NUMFOUR,NUMPHYS,NTHMAX,NEXPTOT,NEXPTOTP,
     2  MEXPNOF,MEXPSOF,MEXPNPHYS,MEXPSPHYS,RDMINUS,
     3  MEXNALL,MEXN1256,MEXN12,MEXN56,MEXSALL,MEXS3478,MEXS34,
     4  MEXS78,XS,YS,ZS,FEXPE,FEXPO,RLSC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine creates the north (+y)  and south (-y) exponential
c     expansions for a parent box due to all eight children.
c
c  note:
c
c     some intelligence is used in the order of summation. thus
c     mexn12 and mexn56 are computed separately. mexn1256 is then
c     obtained by adding these two expansions together, etc.
c
c  on input:
c
c     ibox: current box number.
c     box: current box information.
c     nterms: number of terms in the multipole expansion.
c     mpole: the multipole expansion coefficients.
c     rlams: exponential expansion coefficients.
c     nlambs: number of terms in the exponential expansion.
c     numfour: number of fourier modes in the expansion.
c     numphys: number of modes in the plane wave expansion.
c     nthmax: max number of terms in the exponential expansion.
c     nexptot: total number of fourier modes in the expansion.
c     nextpotp: half of the fourier modes.
c
c  precomputed tables:
c    mrotate: rotation matrix so we shift along z-axis.
c    rdminus:
c    xs,ys,zs: stores the diagonal translation operators when shifting exponential
c        expansions.
c    fexpe,fexp0: how exponential expansions will be merged.
c    rlsc: stores p_n^m for different lambda_k.
c
c  on output:
c    mexnall: up expansion from all boxes.
c    mexn1256: up expansion from box 1-4.
c    mexn12:
c    mexn56:
c    mexsall:
c    mexs3478:
c    mexs34:
c    mexs78:
c
c  variables:
c    mwork:
c    mexpnof:
c    mexpsof:
c    mexpnphys:
c    mexpsphys:
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 LADDER(1),IBOX,NLEV
      INTEGER *4 NTERMS
      INTEGER *4 NLAMBS,NUMFOUR(NLAMBS),NEXPTOT,NTHMAX
      INTEGER *4 NUMPHYS(NLAMBS),NEXPTOTP
c
      REAL *8 ZS(3,NEXPTOTP)
      REAL *8 RDMINUS(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      REAL *8 RLAMS(NLAMBS)
      REAL *8 RLSC(0:NTERMS,0:NTERMS,NLAMBS)
c
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS,*)
      COMPLEX *16 MROTATE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MWORK(0:NTERMS,0:NTERMS)
      COMPLEX *16 MEXPNOF(NEXPTOT),MEXPSOF(NEXPTOT)
      COMPLEX *16 MEXPNPHYS(NEXPTOTP),MEXPSPHYS(NEXPTOTP)
      COMPLEX *16 MEXNALL(NEXPTOTP),MEXN1256(NEXPTOTP)
      COMPLEX *16 MEXN12(NEXPTOTP),MEXN56(NEXPTOTP)
      COMPLEX *16 MEXSALL(NEXPTOTP),MEXS3478(NEXPTOTP)
      COMPLEX *16 MEXS34(NEXPTOTP),MEXS78(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP),YS(3,NEXPTOTP)
      COMPLEX *16 FEXPE(1),FEXPO(1)
c
c-----local variables.
c
      INTEGER *4 JJ,ISTART,MYBOX,ICHILD(8)
      COMPLEX *16 ZTMP
c
      ISTART = LADDER(NLEV)
      MYBOX = ISTART + IBOX
      CALL MKCHILD(NLEV,LADDER,MYBOX,ICHILD)
c
c-----include contributions from child 1
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(1)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXN12(JJ) = MEXPNPHYS(JJ)
        MEXSALL(JJ) = MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 2
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(2)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)*DCONJG(YS(1,JJ))
        MEXN12(JJ) = MEXN12(JJ) + MEXPNPHYS(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*YS(1,JJ)
        MEXSALL(JJ) = MEXSALL(JJ) + MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 3
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(3)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)/ZS(1,JJ)
        MEXNALL(JJ) = MEXPNPHYS(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*ZS(1,JJ)
        MEXS34(JJ) = MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 4
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(4)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = YS(1,JJ)*ZS(1,JJ)
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)/ZTMP
        MEXNALL(JJ) = MEXNALL(JJ) + MEXPNPHYS(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*ZTMP
        MEXS34(JJ) = MEXS34(JJ) + MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 5
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(5)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)*DCONJG(XS(1,JJ))
        MEXN56(JJ) = MEXPNPHYS(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*XS(1,JJ)
        MEXSALL(JJ) = MEXSALL(JJ) + MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 6
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(6)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = YS(1,JJ)*XS(1,JJ)
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)*DCONJG(ZTMP)
        MEXN56(JJ) = MEXN56(JJ) + MEXPNPHYS(JJ)
        MEXN1256(JJ) = MEXN56(JJ) + MEXN12(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*ZTMP
        MEXSALL(JJ) = MEXSALL(JJ) + MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 7
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(7)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = ZS(1,JJ)*XS(1,JJ)
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)/ZTMP
        MEXNALL(JJ) = MEXNALL(JJ) + MEXPNPHYS(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*ZTMP
        MEXS78(JJ) = MEXPSPHYS(JJ)
      ENDDO
c
c-----include contributions from child 8
c
      CALL ROTZTOY(NTERMS,MPOLE(0,0,ICHILD(8)),MWORK,MROTATE,RDMINUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPNOF,MEXPSOF,RLSC)
      CALL YFTOPHYS(MEXPNOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPNPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPSOF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPSPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = ZS(1,JJ)*YS(1,JJ)*XS(1,JJ)
        MEXPNPHYS(JJ) = MEXPNPHYS(JJ)/ZTMP
        MEXNALL(JJ) = MEXNALL(JJ) + MEXPNPHYS(JJ) + MEXN1256(JJ)
        MEXPSPHYS(JJ) = MEXPSPHYS(JJ)*ZTMP
        MEXS78(JJ) = MEXS78(JJ) + MEXPSPHYS(JJ)
        MEXS3478(JJ) = MEXS78(JJ) + MEXS34(JJ)
        MEXSALL(JJ) = MEXSALL(JJ) + MEXS3478(JJ)
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YMKEWEXP(IBOX,LADDER,NLEV,NTERMS,MPOLE,MROTATE,
     1  RLAMS,NLAMBS,NUMFOUR,NUMPHYS,NTHMAX,NEXPTOT,NEXPTOTP,
     2  MEXPEF,MEXPWF,MEXPEPHYS,MEXPWPHYS,RDPLUS,
     3  MEXEALL,MEXE1357,MEXE13,MEXE57,MEXE1,MEXE3,MEXE5,
     4  MEXE7,MEXWALL,MEXW2468,MEXW24,MEXW68,MEXW2,MEXW4,
     5  MEXW6,MEXW8,XS,YS,ZS,FEXPE,FEXPO,RLSC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     this subroutine creates the east (+x)  and west (-x) exponential
c     expansions for a parent box due to all eight children.
c
c  note:
c
c     some intelligence is used in the order of summation. thus
c     mexe1, mexe3, mexe5, mexe7 are computed separately. mexe13, mexe57
c     mex1357 are then obtained by adding these two expansions
c     together, etc.
c
c  on input:
c     ibox: current box number.
c     box: current box information.
c     nterms: number of terms in the multipole expansion.
c     mpole: the multipole expansion coefficients.
c     rlams: exponential expansion coefficients.
c     nlambs: number of terms in the exponential expansion.
c     numfour: number of fourier modes in the expansion.
c     numphys: number of modes in the plane wave expansion.
c     nthmax: max number of terms in the exponential expansion.
c     nexptot: total number of fourier modes in the expansion.
c     nextpotp: half of the fourier modes.
c
c  precomputed tables:
c    mrotate: rotation matrix so we shift along z-axis.
c    rdminus:
c    xs,ys,zs: stores the diagonal translation operators when shifting exponential
c        expansions.
c    fexpe,fexp0: how exponential expansions will be merged.
c    rlsc: stores p_n^m for different lambda_k.
c
c  on output:
c    mexeall: east expansion from all boxes.
c    mexe1357: east expansion from box 1-4.
c    mexe13:
c    mexe57:
c    mexe1:
c    mexe3:
c    mexe5:
c    mexe7:
c    mexwall: west expansion from all boxes.
c    mexw2468: west expansion from box 1-4.
c    mexw24:
c    mexw68:
c    mexw2:
c    mexw4:
c    mexw6:
c    mexw8:
c
c  variables:
c    mwork:
c    mexpnof:
c    mexpsof:
c    mexpnphys:
c    mexpsphys:
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 LADDER(1),IBOX,NLEV
      INTEGER *4 NTERMS
      INTEGER *4 NLAMBS,NUMFOUR(NLAMBS),NEXPTOT,NTHMAX
      INTEGER *4 NUMPHYS(NLAMBS),NEXPTOTP
      INTEGER *4 IEXP(16)
c
      REAL *8 RDPLUS(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      REAL *8 ZS(3,NEXPTOTP)
      REAL *8 RLAMS(NLAMBS)
      REAL *8 RLSC(0:NTERMS,0:NTERMS,NLAMBS)
c
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS,*)
      COMPLEX *16 MROTATE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MEXPEF(NEXPTOT),MEXPWF(NEXPTOT)
      COMPLEX *16 MEXPEPHYS(NEXPTOTP),MEXPWPHYS(NEXPTOTP)
      COMPLEX *16 MEXEALL(NEXPTOTP),MEXE1357(NEXPTOTP)
      COMPLEX *16 MEXE13(NEXPTOTP),MEXE57(NEXPTOTP)
      COMPLEX *16 MEXE1(NEXPTOTP),MEXE3(NEXPTOTP)
      COMPLEX *16 MEXE5(NEXPTOTP),MEXE7(NEXPTOTP)
      COMPLEX *16 MEXWALL(NEXPTOTP),MEXW2468(NEXPTOTP)
      COMPLEX *16 MEXW24(NEXPTOTP),MEXW68(NEXPTOTP)
      COMPLEX *16 MEXW2(NEXPTOTP),MEXW4(NEXPTOTP)
      COMPLEX *16 MEXW6(NEXPTOTP),MEXW8(NEXPTOTP)
      COMPLEX *16 XS(3,NEXPTOTP),YS(3,NEXPTOTP)
      COMPLEX *16 FEXPE(1),FEXPO(1)
c
c-----local variables.
c
      INTEGER *4 JJ,ISTART,MYBOX
      INTEGER *4 ICHILD(8)
      COMPLEX *16 ZTMP
c
      ISTART = LADDER(NLEV)
      MYBOX = ISTART + IBOX
      CALL MKCHILD(NLEV,LADDER,MYBOX,ICHILD)
c
c-----include contributions from child 1
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(1)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXE1(JJ) = MEXPEPHYS(JJ)
        MEXWALL(JJ) = MEXPWPHYS(JJ)
      ENDDO
c
c-----include contributions from child 2
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(2)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXEALL(JJ) = MEXPEPHYS(JJ)/ZS(1,JJ)
        MEXW2(JJ) = MEXPWPHYS(JJ)*ZS(1,JJ)
      ENDDO
c
c-----include contributions from child 3
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(3)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXE3(JJ) = MEXPEPHYS(JJ)*DCONJG(YS(1,JJ))
        MEXE13(JJ) = MEXE1(JJ) + MEXE3(JJ)
        MEXPWPHYS(JJ) = MEXPWPHYS(JJ)*YS(1,JJ)
        MEXWALL(JJ) = MEXWALL(JJ) + MEXPWPHYS(JJ)
      ENDDO
c
c-----include contributions from child 4
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(4)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = YS(1,JJ)*ZS(1,JJ)
        MEXPEPHYS(JJ) = MEXPEPHYS(JJ)/ZTMP
        MEXEALL(JJ) = MEXEALL(JJ) + MEXPEPHYS(JJ)
        MEXW4(JJ) = MEXPWPHYS(JJ)*ZTMP
        MEXW24(JJ) = MEXW2(JJ) + MEXW4(JJ)
      ENDDO
c
c-----include contributions from child 5
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(5)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        MEXE5(JJ) = MEXPEPHYS(JJ)*XS(1,JJ)
        MEXPWPHYS(JJ) = MEXPWPHYS(JJ)*DCONJG(XS(1,JJ))
        MEXWALL(JJ) = MEXWALL(JJ) + MEXPWPHYS(JJ)
      ENDDO
c
c-----include contributions from child 6
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(6)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = XS(1,JJ)/ZS(1,JJ)
        MEXPEPHYS(JJ) = MEXPEPHYS(JJ)*ZTMP
        MEXEALL(JJ) = MEXEALL(JJ) + MEXPEPHYS(JJ)
        MEXW6(JJ) = MEXPWPHYS(JJ)/ZTMP
      ENDDO
c
c-----include contributions from child 7
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(7)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = XS(1,JJ)*DCONJG(YS(1,JJ))
        MEXE7(JJ) = MEXPEPHYS(JJ)*ZTMP
        MEXE57(JJ) = MEXE5(JJ) + MEXE7(JJ)
        MEXE1357(JJ) = MEXE13(JJ) + MEXE57(JJ)
        MEXPWPHYS(JJ) = MEXPWPHYS(JJ)*DCONJG(ZTMP)
        MEXWALL(JJ) = MEXWALL(JJ) + MEXPWPHYS(JJ)
      ENDDO
c
c-----include contributions from child 8
c
      CALL ROTZTOX(NTERMS,MPOLE(0,0,ICHILD(8)),MROTATE,RDPLUS)
      CALL YMPOLETOEXP(MROTATE,NTERMS,NLAMBS,
     1  NUMFOUR,NEXPTOT,MEXPEF,MEXPWF,RLSC)
      CALL YFTOPHYS(MEXPEF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPEPHYS,FEXPE,FEXPO)
      CALL YFTOPHYS(MEXPWF,NLAMBS,RLAMS,NUMFOUR,NUMPHYS,
     1  NTHMAX,MEXPWPHYS,FEXPE,FEXPO)
      DO JJ = 1,NEXPTOTP
        ZTMP = XS(1,JJ)*DCONJG(YS(1,JJ))/ZS(1,JJ)
        MEXPEPHYS(JJ) = MEXPEPHYS(JJ)*ZTMP
        MEXEALL(JJ) = MEXEALL(JJ) + MEXPEPHYS(JJ) + MEXE1357(JJ)
        MEXW8(JJ) = MEXPWPHYS(JJ)/ZTMP
        MEXW68(JJ) = MEXW8(JJ) + MEXW6(JJ)
        MEXW2468(JJ) = MEXW68(JJ) + MEXW24(JJ)
        MEXWALL(JJ) = MEXWALL(JJ) + MEXW2468(JJ)
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YHFSTRTN(NTERMS,D,SQC,THETA)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    implement the fast version of rotation matrices from
c      the recurrences formulas.
c
c  on input:
c    nterms: an integer indicates the dimension of d.
c    sqc: an array contains the square root of the
c       binormial coefficients.
c    theta:  the rotate angle about the y-axis.
c
c  on output:
c    d: an array which contains the rotation matrix.
c
c  note: only half of d are evaluated, the other
c    half can be obtained by using the symmetricity.
c
c  subroutine called : dabs(), dcos(), dsin(), dsqrt()
c
c  called from : rotgen()
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
c
      REAL *8 D(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      REAL *8 SQC(0:4*NTERMS, 0:4*NTERMS)
      REAL *8 THETA
c
c-----local variables.
c
      INTEGER *4 N,M,IJ,IM,IMP,MP,MPABS
      REAL *8 CTHETA, STHETA, HSTHTA
      REAL *8 CTHTAP, CTHTAN, PRECIS
      REAL *8 WW
      REAL *8 FACTS(0:200)
c
      DATA PRECIS/1.0D-19/
      DATA WW/0.7071067811865476D+00/
c
      CTHETA=DCOS(THETA)
      IF (DABS(CTHETA).LE.PRECIS) CTHETA=0.0D0
      STHETA=DSIN(-THETA)
      IF (DABS(STHETA).LE.PRECIS) STHETA=0.0D0
      HSTHTA=WW*STHETA
      CTHTAP=WW*(1.0D0+CTHETA)
      CTHTAN=-WW*(1.0D0-CTHETA)
c
c-----initial setup for some coefficient matrix.
c
      D(0,0,0)=1.0D0
c
      DO IJ=1,NTERMS
c
c-------compute the result for m'=0 case, use formula (1).
c
        DO IM=-IJ,-1
          D(IJ,0,IM)=-SQC(IJ-IM,2)*D(IJ-1,0,IM+1)
          IF (IM.GT.(1-IJ)) THEN
            D(IJ,0,IM)=D(IJ,0,IM)+SQC(IJ+IM,2)*D(IJ-1,0,IM-1)
          ENDIF
          D(IJ,0,IM)=D(IJ,0,IM)*HSTHTA
          IF (IM.GT.-IJ) THEN
            D(IJ,0,IM)=D(IJ,0,IM)+
     1        D(IJ-1,0,IM)*CTHETA*SQC(IJ+IM,1)*SQC(IJ-IM,1)
          ENDIF
          D(IJ,0,IM)=D(IJ,0,IM)/IJ
        ENDDO
c
        D(IJ,0,0)=D(IJ-1,0,0)*CTHETA
        IF (IJ.GT.1) THEN
          D(IJ,0,0)=D(IJ,0,0)+HSTHTA*SQC(IJ,2)*(D(IJ-1,0,-1)+
     1      D(IJ-1,0,1))/IJ
        ENDIF
c
        DO IM=1,IJ
          D(IJ,0,IM)=-SQC(IJ+IM,2)*D(IJ-1,0,IM-1)
          IF (IM.LT.(IJ-1)) THEN
            D(IJ,0,IM)=D(IJ,0,IM)+SQC(IJ-IM,2)*D(IJ-1,0,IM+1)
          ENDIF
          D(IJ,0,IM)=D(IJ,0,IM)*HSTHTA
          IF (IM.LT.IJ) THEN
            D(IJ,0,IM)=D(IJ,0,IM)+
     1        D(IJ-1,0,IM)*CTHETA*SQC(IJ+IM,1)*SQC(IJ-IM,1)
          ENDIF
          D(IJ,0,IM)=D(IJ,0,IM)/IJ
        ENDDO
c
c-------compute the result for 0<m'<=j case, use formula (2).
c
        DO IMP=1,IJ
          DO IM=-IJ,-1
            D(IJ,IMP,IM)=D(IJ-1,IMP-1,IM+1)*CTHTAN*SQC(IJ-IM,2)
            IF (IM.GT.(1-IJ)) THEN
              D(IJ,IMP,IM)=D(IJ,IMP,IM)-
     1          D(IJ-1,IMP-1,IM-1)*CTHTAP*SQC(IJ+IM,2)
            ENDIF
            IF (IM.GT.-IJ) THEN
              D(IJ,IMP,IM)=D(IJ,IMP,IM)+
     1          D(IJ-1,IMP-1,IM)*STHETA*SQC(IJ+IM,1)*SQC(IJ-IM,1)
            ENDIF
            D(IJ,IMP,IM)=D(IJ,IMP,IM)*WW/SQC(IJ+IMP,2)
          ENDDO
c
          D(IJ,IMP,0)=IJ*STHETA*D(IJ-1,IMP-1,0)
          IF (IJ.GT.1) THEN
            D(IJ,IMP,0)=D(IJ,IMP,0)-SQC(IJ,2)*(
     1        D(IJ-1,IMP-1,-1)*CTHTAP+D(IJ-1,IMP-1,1)*CTHTAN)
          ENDIF
          D(IJ,IMP,0)=D(IJ,IMP,0)*WW/SQC(IJ+IMP,2)
c
          DO IM=1,IJ
            D(IJ,IMP,IM)=D(IJ-1,IMP-1,IM-1)*CTHTAP*SQC(IJ+IM,2)
            IF (IM.LT.(IJ-1)) THEN
              D(IJ,IMP,IM)=D(IJ,IMP,IM)-
     1          D(IJ-1,IMP-1,IM+1)*CTHTAN*SQC(IJ-IM,2)
            ENDIF
            IF (IM.LT.IJ) THEN
              D(IJ,IMP,IM)=D(IJ,IMP,IM)+
     1          D(IJ-1,IMP-1,IM)*STHETA*SQC(IJ+IM,1)*SQC(IJ-IM,1)
            ENDIF
            D(IJ,IMP,IM)=D(IJ,IMP,IM)*WW/SQC(IJ+IMP,2)
          ENDDO
c
c---------note: the lower part of the matrix can be computed using
c           symmetricity, i.e. formula (3.80) in biedenharn & louck's
c           book.
c
        ENDDO
      ENDDO
c
c-----now scale the rotation matrix to avoid y_n^m
c       note : since in yukawa, i will use p_n^m instead of
c            y_n^m.
c
      FACTS(0)=1.0D0
      DO 1 N=1, 2*NTERMS
        FACTS(N)=FACTS(N-1)*DBLE(N)
1     CONTINUE
      DO N=0, NTERMS
        DO M=0, N
          DO MP=-N, N
            MPABS=IABS(MP)
            D(N,M,MP)=D(N,M,MP)*DSQRT(FACTS(N+M)/FACTS(N+MPABS)
     1        *FACTS(N-MPABS)/FACTS(N-M) )
          ENDDO
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ROTZTOY(NTERMS,MPOLE,MWORK,MROTATE,RDMINUS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    the rotation matrix used in the subroutine mknsexp()
c    so the north_south expansions are made the same way
c    as the up-down expansions.
c
c  on input :
c
c    nterms : number of terms in the multipole expansion.
c    mpole : the multipole expansion.
c    rdminus : the rotation matrix generated in subroutine
c      rotgen<- fstrtn()
c
c  output :
c    mrotate : the rotated multiple expansion coefficients.
c
c  working space : mwork().
c
c  called from :  mknsexp()
c
c  subroutine called : none.
c
c  note : this can be further simplified?
c
c     end result         z_new <- y_old
c                        y_new <- x_old
c                        x_new <- z_old
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 RDMINUS(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MWORK(0:NTERMS,0:NTERMS)
      COMPLEX *16 MROTATE(0:NTERMS,0:NTERMS)
c
c-----local varialbles.
c
      INTEGER *4 M,L,MP
      COMPLEX *16 EPHI(0:100)
c
c-----a rotation of -pi/2 radians about the z-axis in the
c       original coordinate system.
c
      EPHI(0) = 1.0D0
      DO M=1,NTERMS
        EPHI(M)=EPHI(M-1)*DCMPLX(0.0D0,-1.0D0)
      ENDDO
c
      DO L=0,NTERMS
        DO M=0,L
          MWORK(L,M)=EPHI(M)*MPOLE(L,M)
        ENDDO
      ENDDO
c
c-----a rotation of -pi/2 radians about the y'-axis in the
c       new coordinate system, bringing the +z-axis in line
c       with the original +y axis.
c
      DO L=0,NTERMS
        DO M=0,L
          MROTATE(L,M)=MWORK(L,0)*RDMINUS(L,0,M)
          DO MP=1,L
            MROTATE(L,M)=MROTATE(L,M)+MWORK(L,MP)*RDMINUS(L,MP,M)+
     1        DCONJG(MWORK(L,MP))*RDMINUS(L,MP,-M)
          ENDDO
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ROTYTOZ(NTERMS,MPOLE,MWORK,MROTATE,RDPLUS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    the rotation matrix used in the subroutine ladapfmm(), yadapfmm()
c      rotate the totated north_south expansions back to the
c      normal direction.
c
c  on input :
c    nterms : number of terms in the multipole expansion.
c    mpole : the local expansion.
c    rdplus : the rotation matrix generated in subroutine
c      rotgen<- fstrtn()
c
c  on output :
c    mrotate : the rotated local expansion coefficients.
c
c  working space : mwork().
c
c  called from :  ladapfmm(), yadapfmm()
c  subroutine called : none.
c
c  note : this can be further simplified?
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 RDPLUS(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MWORK(0:NTERMS,0:NTERMS)
      COMPLEX *16 MROTATE(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M,MP
      COMPLEX *16 EPHI(0:100)
c
c----a rotation of pi/2 radians about the y'-axis in the
c      new coordinate system, bringing the +z-axis in line
c      with the original +z axis.
c
      DO L=0,NTERMS
        DO M=0,L
          MWORK(L,M)=MPOLE(L,0)*RDPLUS(L,0,M)
          DO MP=1,L
            MWORK(L,M)=MWORK(L,M)+MPOLE(L,MP)*RDPLUS(L,MP,M)+
     1        DCONJG(MPOLE(L,MP))*RDPLUS(L,MP,-M)
          ENDDO
        ENDDO
      ENDDO
c
c-----a rotation of pi/2 radians about the z-axis in the
c       original coordinate system.
c
      EPHI(0) = 1.0D0
      DO M=1,NTERMS
        EPHI(M)=EPHI(M-1)*DCMPLX(0.0D0,1.0D0)
      ENDDO
c
      DO L=0,NTERMS
        DO M=0,L
          MROTATE(L,M)=EPHI(M)*MWORK(L,M)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ROTZTOX(NTERMS,MPOLE,MROTATE,RD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    the rotation matrix used in the subroutine ladapfmm(), yadapfmm(),
c    and mkewexp(). rotate the east-west expansions from or to the
c    normal direction, depending on rd.
c
c  on input :
c
c    nterms : number of terms in the multipole/local expansion.
c    mpole : the multipole/local expansion.
c    rd : the rotation matrix generated in subroutine
c         rotgen<- fstrtn(), it can either be rdplus or rdminus
c
c  on output :
c    mrotate : the rotated local expansion coefficients.
c
c  called from :  mkewexp(), ladapfmm(), yadapfmm()
c
c  subroutine called : none.
c
c  note : this can be further simplified?
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 RD(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
      COMPLEX *16 MROTATE(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M,MP
c
      DO L=0,NTERMS
        DO M=0,L
          MROTATE(L,M)=MPOLE(L,0)*RD(L,0,M)
          DO MP=1,L
            MROTATE(L,M)=MROTATE(L,M)+MPOLE(L,MP)*RD(L,MP,M)+
     1        DCONJG(MPOLE(L,MP))*RD(L,MP,-M)
          ENDDO
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YHFRMINI(NTERMS,C)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c----- initialization entry point - create factorial scaling factors
c
c     ...... not the most stable way of doing this ......
c
c     the vector fact() is only used here. and c(l,m) will used
c       in subroutine form_mp() and entry brtaev()
c
c     note :
c       the current program is changed a little bit from entry to
c       a subroutine, since the result c is needed at several places.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 C(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M
      REAL *8 FACT(0:120)
c
c-----functions called
c
      REAL *8 DBLE
c
      PRINT *, ' initializing'
      FACT(0) = 1.0D0
c
      DO 6000 L=1,2*NTERMS+1
        FACT(L) = FACT(L-1)*DBLE(L)
6000  CONTINUE
c
      DO 6200 L=0,NTERMS
        DO 6100 M = 0,L
          C(L,M) = FACT(L-M)/FACT(L+M)*DBLE(2*L+1)
6100    CONTINUE
6200  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YHROTGEN(NTERMS,CARRAY,RDPI2,RDMPI2,RDSQ3,RDMSQ3)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    precomputes the rotation matrix for
c     1. mp->mp
c     2. local->local
c     3. east-west expansion
c     4. north-south expansion
c
c  on input :
c    nterms : the number of terms in the multipole expansion.
c
c  on output :
c    rdpi2, rdmpi2 : the rotation matrix for 3 and 4.
c    rdsq3, rdmsq3 : the rotation matrix for 1 and 2.
c
c  workspace :
c    carray : the square root of the binomial numbers.
c      these numbers are only used here in this subroutine.
c
c  subroutine called :
c    bnlcft(), fstrtn(), datan, dacos,
c
c  called from : main()
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 CARRAY(0:4*NTERMS,0:4*NTERMS)
      REAL *8 RDPI2(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      REAL *8 RDMPI2(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      REAL *8 RDSQ3(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
      REAL *8 RDMSQ3(0:NTERMS,0:NTERMS,-NTERMS:NTERMS)
c
c-----local varialbles.
c
      REAL *8 PI,THETA
c
c-----call initialization routines
c
      PI = 4.0D0*DATAN(1.0D0)
c
      CALL BNLCFT(CARRAY,4*NTERMS)
      THETA = PI/2.0D0
      CALL YHFSTRTN(NTERMS,RDPI2,CARRAY,THETA)
      THETA = -PI/2.0D0
      CALL YHFSTRTN(NTERMS,RDMPI2,CARRAY,THETA)
      THETA = DACOS(DSQRT(3.0D0)/3.0D0)
      CALL YHFSTRTN(NTERMS,RDSQ3,CARRAY,THETA)
      THETA = DACOS(-DSQRT(3.0D0)/3.0D0)
      CALL YHFSTRTN(NTERMS,RDMSQ3,CARRAY,THETA)
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ADJ_MEM_PTR(MEM_PTR,ASSIGN_PTR,SIZE)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    computes the total required memory and the pointer for the
c    memory chunk.
c
c  on input:
c    mem_ptr: the previous pointer to the end of the memory chunk.
c    size: the size of the memory to be added.
c
c  on output:
c    assign_ptr: the pointer to the new memory request.
c    mem_ptr: the updated total memory requirement.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER*4 MEM_PTR,ASSIGN_PTR,SIZE
c
      ASSIGN_PTR = MEM_PTR
      MEM_PTR = MEM_PTR + SIZE
c
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ICOPY(N,DX,DY)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    copies an integer vector x of length n to another vector y of
c      the same length.
c
c  on input:
c    n: the number of elements to be copied.
c    dx: the vector to be copied.
c
c  on output:
c    dy: the duplicated vector.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 N
      INTEGER *4 DX(N),DY(N)
      INTEGER I
c
      IF(N.LE.0)RETURN
c
      DO 10 I = 1,N
        DY(I) = DX(I)
10    CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ADDEXP(B,A,NTERMS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose: adds one vector to another.
c
c  note: this can be replaced by future blast subroutines.
c
c  on input: expansions a and b of size (0:nterms,0:nterms)
c
c  on output: a is over written by (a+b).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4  NTERMS
c
      COMPLEX *16 B(0:NTERMS,0:NTERMS),A(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 I,J
c
      DO I = 0,NTERMS
        DO J = 0,NTERMS
          A(I,J) = A(I,J) + B(I,J)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE BNLCFT(C, NTERMS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    computes the binomial coefficients c_nterms^n, where n=0,1,2,...,nterms.
c
c  on input:
c
c    nterms: an integer indicates the number we are going to choose from.
c
c  on output:
c
c    c:    an array consists of the squre root of the
c                 binomial coefficients.
c
c  note : this is a different version from the laplace equation.
c    since the binomial coefficients are not needed, but
c    only the square root is needed, so we will not store
c    the binomial coefficients. and we will use that space
c    for some must-be computed coefficients at different levels.
c
c  subroutine called : dsqrt()
c
c  called from : rotgen()
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      REAL *8 C(0:NTERMS,0:NTERMS)
c
      INTEGER *4 N,M
c
c-----compute c(n,m)
c
      DO N=0,NTERMS
        C(N,0)=1.0D0
      ENDDO
c
      DO M=1,NTERMS
        C(M,M)=1.0D0
        DO N=M+1,NTERMS
          C(N,M)=C(N-1,M)+C(N-1,M-1)
        ENDDO
      ENDDO
c
c-----compute the square root of c(n,m)
c
      DO M=1,NTERMS
        DO N=M+1,NTERMS
          C(N,M)=DSQRT(C(N,M))
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE WRITETIMES(NATOMS,NLEV,NTERMS,TIME_UPPASS,TMKEXPS,
     1  TSHIFTS,TIME_EXPTOMP,TIME_TATA,TIME_LCEVAL)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose: output the cpu time information for the fmm.
c
c  note: this subroutine is only used by the uniform code.
c    therefore it is ok to put it here.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NATOMS,NLEV,NTERMS
      REAL *8 TIME_UPPASS,TMKEXPS,TSHIFTS,TIME_EXPTOMP,TIME_TATA
      REAL *8 TIME_LCEVAL
c
      WRITE(11,*)'n     nlev  nterms up pass mkexps  shifts  ',
     1              'exptomp tata    lceval'
      WRITE(11,1001)NATOMS,NLEV,NTERMS,TIME_UPPASS,TMKEXPS,
     1                TSHIFTS,TIME_EXPTOMP,TIME_TATA,TIME_LCEVAL
1001  FORMAT(I7,2X,I2,4X,I2,1X,6F8.2)
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  the files for the driver.
c    subroutines in this file are used for both adaptive and uniform
c    code.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE DUMMY1(NMOLS,ZAT,CHARGE)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose: sets the particle distributions, either random, or on a
c    sphere, or other types for adaptive fmm testing.
c
c  subroutine called : rand()
c
c  called from : main()
c
c  on input :
c
c    nmols : the number of random particles to be generated.
c
c  on output :
c
c    zat(3,:) : the location of the particles.
c    charge : the charge of the particles.
c
c  note: for adaptive code, particles can be located arbitrary.
c        for uniform, only unit box is allowed.
c
c  RE-written by Tsamouris Panayiotis for benchmarking purposes
c  added the functionality of binary dump input
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NMOLS
      REAL *8 VAL
      REAL *8 ZAT(3,NMOLS),CHARGE(NMOLS)
      REAL *8 X(NMOLS)
      REAL *8 Y(NMOLS)
      REAL *8 Z(NMOLS)
c
c-----local variables
c
      INTEGER *4 J
      REAL *4 RAND
      REAL *8 THETA, PHI
      

      OPEN (UNIT = 8 , FILE = "binary_x", FORM = 'unformatted',
     1      access = 'direct' , recl = 8)
      OPEN (UNIT = 7 , FILE = "binary_y", FORM = 'unformatted',
     1      access = 'direct' , recl = 8)
      OPEN (UNIT = 6 , FILE = "binary_z", FORM = 'unformatted',
     1      access = 'direct' , recl = 8)
      OPEN (UNIT = 5 , FILE = "binary_charge", FORM = 'unformatted',
     1      access = 'direct' , recl = 8)      

      OPEN (UNIT = 9, FILE = "test")

       do   J = 1, NMOLS
       read( 8 , rec = J ) ZAT(1,J)
       read( 7 , rec = J ) ZAT(2,J)
       read( 6 , rec = J ) ZAT(3,j)
       read( 5 , rec = J ) CHARGE(J)       
       WRITE (9,*) ZAT(1,J) ,ZAT(2,J) ,ZAT(3,J) ,CHARGE(J) 
	end do

c debugging write to file this works for text input
c     DO 1000 J = 1,NMOLS

c
c      READ (8,*) ZAT(1,J) ,ZAT(2,J) ,ZAT(3,J) ,CHARGE(J)
c      WRITE (9,*) ZAT(1,J) ,ZAT(2,J) ,ZAT(3,J) ,CHARGE(J)
c1000  CONTINUE
       
      CLOSE (9)
      CLOSE (5)
      CLOSE (6)
      CLOSE (7)
      CLOSE (8) 
c
      RETURN
      END
c

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE DUMMY(NMOLS,ZAT,CHARGE)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose: sets the particle distributions, either random, or on a
c    sphere, or other types for adaptive fmm testing.
c
c  subroutine called : rand()
c
c  called from : main()
c
c  on input :
c
c    nmols : the number of random particles to be generated.
c
c  on output :
c
c    zat(3,:) : the location of the particles.
c    charge : the charge of the particles.
c
c  note: for adaptive code, particles can be located arbitrary.
c        for uniform, only unit box is allowed.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NMOLS
      REAL *8 VAL
      REAL *8 ZAT(3,NMOLS),CHARGE(NMOLS)
c
c-----local variables
c
      INTEGER *4 J
      REAL *4 RAND
      REAL *8 THETA, PHI
c
c-----initial seed.
c
      VAL = RAND(1)
c
c-----distribute molecules randomly in box and assign random charges.
c
c     note: for uniform code, the particles are inside the unit box.
c       for nonuniform code, particles can be arbitrary.
c
      DO 1000 J = 1,NMOLS
c
c-------the following generates a distribution on a sphere.
c
      theta=rand(0)*3.1415926d0
      phi=rand(0)*2.0d0*3.1415926d0
      zat(1,j) =1.0d0*dsin(theta)*dsin(phi)
      zat(2,j) =1.0d0*dsin(theta)*dcos(phi)
      zat(3,j) =1.0d0*dcos(theta)
c
c-------the following generates a distribution inside a cylinder.
c
c        ZAT(1,J) =0.10D0*( RAND(0)-0.5D0)
c        ZAT(2,J) =( RAND(0)-0.5D0)
c        ZAT(3,J) =( RAND(0)-0.5D0)
c
c-------generate the charge distribution.
c
        CHARGE(J) = RAND(0) - 0.5D0
1000  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YDIRECI(BETA,NPARTS,ZPARTS,CHARGE,I,
     1                   RPOT,FIELD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    calculates directly the potential for particle i.
c
c  on input :
c
c    beta : the frequency. lap u-beta^2 u=f
c    nparts : the total number of particles.
c    zparts(3,:) : the location of the particles.
c    charge : the charge of the particle.
c    i : the location of the target particle.
c
c  on output :
c    rpot : the potential.
c
c  subroutine called : dsqrt, dexp,
c
c  called from : main()
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      IMPLICIT NONE
c
      INTEGER *4 NPARTS,I
      REAL *8 BETA,ZPARTS(3,NPARTS),CHARGE(NPARTS)
      REAL *8 RPOT,FIELD(3)
c
      INTEGER *4 J
      REAL *8 PI,RX,RY,RZ,RR,RDIS,EXPR,TERM1,TERM2
c
      PI = 4.0D0*DATAN(1.0D0)
      RPOT = 0.0D0
      FIELD(1)=0.0D0
      FIELD(2)=0.0D0
      FIELD(3)=0.0D0
c

      DO 1000 J = 1,NPARTS
        IF (J .EQ. I) GOTO 1000
        RX = ZPARTS(1,I) - ZPARTS(1,J)
        RY = ZPARTS(2,I) - ZPARTS(2,J)
        RZ = ZPARTS(3,I) - ZPARTS(3,J)
        RR = RX*RX + RY*RY + RZ*RZ
        RDIS = DSQRT(RR)*BETA
        EXPR= DEXP(-RDIS)
        TERM1 = CHARGE(J)/RDIS*EXPR*PI/2.0D0
        RPOT = RPOT + TERM1
c
c-------the field.
c
        TERM2 = -TERM1*(1.0D0+RDIS)/RR
        FIELD(1)=FIELD(1)+RX*TERM2
        FIELD(2)=FIELD(2)+RY*TERM2
        FIELD(3)=FIELD(3)+RZ*TERM2

1000  CONTINUE

c
c
      RETURN
      END
c
