cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  this file contains the translation operators for both uniform and
c   adaptive fmm yukawa solvers.
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
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
ccc        theta=rand(0)*3.1415926d0
ccc        phi=rand(0)*2.0d0*3.1415926d0
ccc        zat(1,j) =1.0d0*dsin(theta)*dsin(phi)
ccc        zat(2,j) =1.0d0*dsin(theta)*dcos(phi)
ccc        zat(3,j) =1.0d0*dcos(theta)
c
c-------the following generates a distribution inside a cylinder.
c
        ZAT(1,J) =0.10D0*( RAND(0)-0.5D0)
        ZAT(2,J) =( RAND(0)-0.5D0)
        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 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
c      REAL *4 RAND
c      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 1000  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) 
1000   CONTINUE

      CLOSE (8)
      CLOSE (7)
      CLOSE (6)
      CLOSE (5)
      CLOSE (9) 
      
c      OPEN (UNIT = 8, FILE = "fortran.dump")
c      OPEN (UNIT = 9, FILE = "test")

c
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
      
c      CLOSE (8)
c      CLOSE (9) 
c
      
      RETURN
      END
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE DNDUMMY(NMOLS,ZAT,CHARGE,DN)
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    dn(3,:) : the normal direction.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NMOLS
      REAL *8 VAL
      REAL *8 ZAT(3,NMOLS),CHARGE(NMOLS),DN(3,NMOLS)
c
c-----local variables
c
      INTEGER *4 J
      REAL *4 RAND
      REAL *8 THETA,PHI,XX,YY,ZZ,RR
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.1D0*( RAND(0)-0.5D0)
c        ZAT(2,J) =0.1D0*( 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
c
c-------generate the normal derivative vector.
c
        XX =( RAND(0)-0.5D0)
        YY =( RAND(0)-0.5D0)
        ZZ =( RAND(0)-0.5D0)
        RR=DSQRT(XX*XX+YY*YY+ZZ*ZZ)
        DN(1,J)=XX/RR
        DN(2,J)=YY/RR
        DN(3,J)=ZZ/RR
ccc        dn(1,j)=0.0d0
ccc        dn(2,j)=0.0d0
ccc        dn(3,j)=zz/rr
c
1000  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE DUMMYST(NSOU,ZSOU,CHARGE,NTAR,ZTAR)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose: sets the source and target particle distributions, either
c    random, or on a sphere, or other types for adaptive fmm testing.
c
c  subroutine called : rand()
c
c  called from : main()
c
c  on input :
c
c    nsou : the number of sourc3 particles to be generated.
c    ntar : the number of target particles to be generated.
c
c  on output :
c    zsou(3,:) : the location of source particles.
c    charge : the charge of the source particles.
c    ztar(3,:) : the location of source 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 NSOU,NTAR
      REAL *8 VAL
      REAL *8 ZSOU(3,NSOU),CHARGE(NSOU),ZTAR(3,NTAR)
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,NSOU
c
c-------the following generates a distribution on a sphere.
c
ccc        theta=rand(0)*3.1415926d0
ccc        phi=rand(0)*2.0d0*3.1415926d0
ccc        zat(1,j) =1.0d0*dsin(theta)*dsin(phi)
ccc        zat(2,j) =1.0d0*dsin(theta)*dcos(phi)
ccc        zat(3,j) =1.0d0*dcos(theta)
c
c-------the following generates a distribution inside a cylinder.
c
        ZSOU(1,J) =1D0*( RAND(0)-0.5D0)
        ZSOU(2,J) =( RAND(0)-0.5D0)
        ZSOU(3,J) =( RAND(0)-0.5D0)
c
c-------generate the charge distribution.
c
        CHARGE(J) = RAND(0) - 0.5D0
1000  CONTINUE
c
      DO 2000 J = 1,NTAR
c
c-------the following generates a distribution on a sphere.
c
ccc        theta=rand(0)*3.1415926d0
ccc        phi=rand(0)*2.0d0*3.1415926d0
ccc        zat(1,j) =1.0d0*dsin(theta)*dsin(phi)
ccc        zat(2,j) =1.0d0*dsin(theta)*dcos(phi)
ccc        zat(3,j) =1.0d0*dcos(theta)
c
c-------the following generates a distribution inside a cylinder.
c
        ZTAR(1,J) =0.1D0*( RAND(0)-0.5D0)
        ZTAR(2,J) =( RAND(0)-0.5D0)
        ZTAR(3,J) =( RAND(0)-0.5D0)
2000  CONTINUE
cccc---debug
c       CALL PRINF('nsou is *', NSOU,1)
c       CALL PRINF('ntar is *', NTAR,1)
c       ZSOU(1,1)=0D0
c       ZSOU(2,1)=0D0
c       ZSOU(3,1)=0D0
c       ZTAR(1,1)=120D0
c       ZTAR(2,1)=0D0
c       ZTAR(3,1)=0D0
c       CHARGE(1)=50D0
cccc----debug
c
      RETURN
      END
c

      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 PRINM(MPOLE,NTERMS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  print out coefficients of multipole expansion
c
c  note: this subroutine is mostly used for debugging purposes.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NTERMS
      COMPLEX *16 MPOLE(0:NTERMS,0:NTERMS)
c
      INTEGER *4 L,M
c
      DO 100 L = 0,NTERMS
         WRITE(6,1000)(MPOLE(L,M),M=0,L)
         WRITE(13,1000)(MPOLE(L,M),M=0,L)
         WRITE(6,1001)
         WRITE(13,1001)
100   CONTINUE
1000  FORMAT(6D12.5)
1001  FORMAT(/)
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
      SUBROUTINE FSTRTN(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:2*NTERMS,0:2*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
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---------use symmetricity, i.e. formula (3.80) in biedenharn & louck's
c           book, to compute the lower part of the matrix
c
        ENDDO
      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 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 LGNDR(NMAX, X, Y)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine computes the lengre polynomial expansion using
c    a recursive expansion.
c
c  on input:
c    nmax: the max number of terms in the expansion.
c    x: where we want to evaluate the expansion.
c
c  on output:
c    y: the function value at x.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NMAX
      REAL *8 X,Y(0:NMAX,0:NMAX)
c
c-----local variables.
c
      INTEGER *4 M,N
      REAL *8 U,DSQRT,DBLE
c
      U=-DSQRT(1.0D0-X*X)
      Y(0,0)=1.0D0
      DO 10 M=0, NMAX
        IF (M.GT.0)  Y(M,M)=Y(M-1,M-1)*U*DBLE(2*M-1)
        IF (M.LT.NMAX)  Y(M+1,M)=DBLE(2*M+1)*X*Y(M,M)
        DO 20 N=M+2, NMAX
          Y(N,M)=((2.0D0*DBLE(N)-1.0D0)*X*Y(N-1,M)-DBLE(N+M-1)
     1      *Y(N-2,M)) / DBLE(N-M)
 20     CONTINUE
 10   CONTINUE
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
c
c  the following files are for the case when both negative and positive
c    terms are stored for the expansion.
c
c  this is important for helmholtz kernels because half the terms are
c    used for the multipole expansion, but all are used for the local
c    expansion.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
         SUBROUTINE LOCADD(B,A,NTERMS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose: adds one vector to another.
c
c  note: this subroutine is used by the helmholtz equation solver.
c        as the m in the local expansion is now from -nterms to nterms
c
c  input :  local expansions a and b of size (0:nterms,-nterms:nterms)
c
c  output:     a is over written by (a+b).
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4  NTERMS
      COMPLEX *16 B(0:NTERMS,-NTERMS:NTERMS),A(0:NTERMS,-NTERMS:NTERMS)
c
      INTEGER *4 I,J
c
      DO I = 0,NTERMS
        DO J = -NTERMS,NTERMS
          A(I,J) = A(I,J) + B(I,J)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE LOCROTYTOZ(NTERMS,MPOLE,MWORK,MROTATE,RDPLUS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    the rotation matrix used in the subroutine hbrfrc()
c      rotate the totated north_south expansions back to the
c      normal direction.
c
c  input :
c    nterms : number of terms in the multipole expansion.
c    mpole : the local expansion.
c    rdplus : the rotation matrix generated in subroutine
c             hrotgen<- hfstrtn()
c
c  output :
c    mrotate : the rotated local expansion coefficients.
c
c  working space : mwork().
c
c  called from :  hbrfrc()
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,-NTERMS:NTERMS)
      COMPLEX *16 MWORK(0:NTERMS,-NTERMS:NTERMS)
      COMPLEX *16 MROTATE(0:NTERMS,-NTERMS:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M,MP
      COMPLEX *16 EPHI(-100: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=-L,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        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)
        EPHI(-M)=DCONJG(EPHI(M))
      ENDDO
c
      DO L=0,NTERMS
        DO M=-L,L
          MROTATE(L,M)=EPHI(M)*MWORK(L,M)
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE LOCROTZTOX(NTERMS,MPOLE,MROTATE,RD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    the rotation matrix used in the subroutine hbrfrc(), and hmkewexp()
c      rotate the east-west expansions from or to the
c      normal direction, depending on rd.
c
c  input :
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         hrotgen<- hfstrtn(), it can either be rdplus or rdminus
c
c  output :
c    mrotate : the rotated local expansion coefficients.
c
c  called from :  hmkewexp(), hbrfrc()
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,-NTERMS:NTERMS)
      COMPLEX *16 MROTATE(0:NTERMS,-NTERMS:NTERMS)
c
c-----local variables.
c
      INTEGER *4 L,M,MP
c
      DO L=0,NTERMS
        DO M=-NTERMS,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        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 FRMINI(C,CS,CSINV)
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       csinv and cs will be used in brtaev only.
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
      REAL *8 C(0:60,0:60),CS(0:60,0:60),CSINV(0:60,0:60)
c
c-----local variables.
c
      INTEGER *4 L,M
      REAL *8 D,FACT(0:120)
c
      D = 1.0D0
      FACT(0) = D
      DO 6000 L=1,120
        D=D*DSQRT(L+0.0D0)
        FACT(L) = D
6000  CONTINUE
c
      CS(0,0) = 1.0D0
      CSINV(0,0) = 1.0D0
      DO 6200 L=1,60
        DO 6100 M = 0,L
          C(L,M) = FACT(L-M)/FACT(L+M)
          CSINV(L,M) = FACT(L-M)*FACT(L+M)
          CS(L,M) = 1.0D0/CSINV(L,M)
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 ROTGEN(NTERMS,CARRAY,RDPI2,RDMPI2,RDSQ3,RDMSQ3,DC)
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:2*NTERMS,0:2*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)
      REAL *8 DC(0:2*NTERMS,0:2*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(DC,2*NTERMS)
c
      THETA = PI/2.0D0
      CALL FSTRTN(NTERMS,RDPI2,DC,THETA)
      THETA = -PI/2.0D0
      CALL FSTRTN(NTERMS,RDMPI2,DC,THETA)
      THETA = DACOS(DSQRT(3.0D0)/3.0D0)
      CALL FSTRTN(NTERMS,RDSQ3,DC,THETA)
      THETA = DACOS(-DSQRT(3.0D0)/3.0D0)
      CALL FSTRTN(NTERMS,RDMSQ3,DC,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

      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
        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^{2\pi} e^{i\lambda(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^{2\pi} e^{-i\lambda(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 IN(SCAL, X, NB, B, NCALC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    calculates the i_n(z)=sqrt(pi/2/z)*i_(n+1/2)(z).
c
c  on input:
c    scal: the scaling factor to avoid underflow.
c    x: the parameter for i_n(x)
c    nb: number of terms for the subindex n, n=0:nb.
c
c  on output :
c    b: contains the i_n(x), n=0:nb.
c    ncalc: an error index. this is from the subroutine
c           for the calculation of bessel function.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NB,NCALC
      REAL *8 X,B(0:NB),SCAL
c
c-----local variables.
c
      INTEGER *4 I,IZE
      REAL *8 CONST,ENSIG,ENMTEN,HALFPI
      REAL *8 XSCAL,TERM1,TERM2,ALPHA
      DATA ENSIG, ENMTEN/1.0D-4, 1.0D-300/
c
c-----function called
c
      REAL *8 DABS
c
      IF (X.LT. 0.0D0) THEN
        PRINT *,' error in input for in.f'
      ENDIF
c
c-----if x.le. 10e-4, then use the 2 terms taylor expansion.
c
      IF (X.LE.ENSIG) THEN
        XSCAL=X/SCAL
        TERM1=1.0D0
        TERM2=0.5D0*X*X
        B(0)=TERM1*(1.0D0+TERM2/3.0D0)
        DO 1 I=1, NB
          TERM1=TERM1*XSCAL/DBLE(2*I+1)
          IF (TERM1.LE. ENMTEN) THEN
            TERM1=0.0D0
          ENDIF
          B(I)=TERM1*(1.0D0+TERM2/DBLE(2*I+3) )
1       CONTINUE
        NCALC=NB+1
c
c-------usually, scal should be less than one, however, in the
c         calculation of h_n, scal will be greater than one,
c         in this case, we should modify the program and replace every
c         small item with zero.
c
      ELSEIF (X.GT.1.0D+2) THEN
        DO I=0,NB
          B(I)=0.0D0
        ENDDO
        NCALC=NB+1
      ELSE
c
c-------otherwise, we will call besselj() and then scale the
c         jn by the scaling factor.
c
        HALFPI=DATAN(1.0D0)*2.0D0
        CONST=DSQRT(HALFPI/X)
c
c-------calculate i_(n+1/2).
c
        ALPHA=0.5D0
        IZE=1
        CALL RIBESL(X, ALPHA, NB+1, IZE, B, NCALC)
c
c-------now scale to i_n*scal**(-n).
c
        DO 2 I=0,NB
          B(I)=B(I)*CONST
          CONST=CONST/SCAL
          IF (DABS(B(I)) .LE. ENMTEN ) THEN
            CONST=0.0D0
          ENDIF
2       CONTINUE
      ENDIF
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE KN(SCAL,X,NB,BY,NCALC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c purpose:
c  this routine calculates scaled sphecial bessel functions y_n (x)
c  for non-negative argument x, and n=0,...,nb.
c
c on input:
c   scal: the scaling factor.
c   x: working precision non-negative real argument for which
c      y's are to be calculated.
c   nb: integer number of functions to be calculated, nb .gt. 0.
c       the first function calculated is of order 0, and the
c       last is of order nb.
c   by: working precision output vector of length nb+1.  if the
c       routine terminates normally (ncalc=nb+1), the vector by
c       contains the functions y(0,x), ... , y(nb,x),
c       if (0 .lt. ncalc .lt. nb), by(i) contains correct function
c       values for i .le. ncalc, and contains the ratios
c       y(i-1,x)/y(i-2,x) for the rest of the array.
c   ncalc: integer output variable indicating possible errors.
c       before using the vector by, the user should check that
c       ncalc=nb, i.e., all orders have been calculated to
c       the desired accuracy.  see error returns below.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c explanation of machine-dependent constants
c
c   beta   = radix for the floating-point system
c   p      = number of significant base-beta digits in the
c            significand of a floating-point number
c   minexp = smallest representable power of beta
c   maxexp = smallest power of beta that overflows
c   eps    = beta ** (-p)
c   del    = machine number below which sin(x)/x = 1; approximately
c            sqrt(eps).
c   xmin   = smallest acceptable argument for rbesy; approximately
c            max(2*beta**minexp,2/xinf), rounded up
c   xinf   = largest positive machine number; approximately
c            beta**maxexp
c   thresh = lower bound for use of the asymptotic form; approximately
c            aint(-log10(eps/2.0))+1.0
c   xlarge = upper bound on x; approximately 1/del, because the sine
c            and cosine functions have lost about half of their
c            precision at that point.
c
c
c     approximate values for some important machines are:
c
c                        beta    p     minexp      maxexp      eps
c
c  cray-1        (s.p.)    2    48     -8193        8191    3.55e-15
c  cyber 180/185
c    under nos   (s.p.)    2    48      -975        1070    3.55e-15
c  ieee (ibm/xt,
c    sun, etc.)  (s.p.)    2    24      -126         128    5.96e-8
c  ieee (ibm/xt,
c    sun, etc.)  (d.p.)    2    53     -1022        1024    1.11d-16
c  ibm 3033      (d.p.)   16    14       -65          63    1.39d-17
c  vax           (s.p.)    2    24      -128         127    5.96e-8
c  vax d-format  (d.p.)    2    56      -128         127    1.39d-17
c  vax g-format  (d.p.)    2    53     -1024        1023    1.11d-16
c
c
c                         del      xmin      xinf     thresh  xlarge
c
c cray-1        (s.p.)  5.0e-8  3.67e-2466 5.45e+2465  15.0e0  2.0e7
c cyber 180/855
c   under nos   (s.p.)  5.0e-8  6.28e-294  1.26e+322   15.0e0  2.0e7
c ieee (ibm/xt,
c   sun, etc.)  (s.p.)  1.0e-4  2.36e-38   3.40e+38     8.0e0  1.0e4
c ieee (ibm/xt,
c   sun, etc.)  (d.p.)  1.0d-8  4.46d-308  1.79d+308   16.0d0  1.0d8
c ibm 3033      (d.p.)  1.0d-8  2.77d-76   7.23d+75    17.0d0  1.0d8
c vax           (s.p.)  1.0e-4  1.18e-38   1.70e+38     8.0e0  1.0e4
c vax d-format  (d.p.)  1.0d-9  1.18d-38   1.70d+38    17.0d0  1.0d9
c vax g-format  (d.p.)  1.0d-8  2.23d-308  8.98d+307   16.0d0  1.0d8
c
c*******************************************************************
c
c error returns
c
c  in case of an error, ncalc .ne. nb, and not all y's are
c  calculated to the desired accuracy.
c
c  ncalc .lt. -1:  an argument is out of range. for example,
c       nb .le. 0, ize is not 1 or 2, or ize=1 and abs(x) .ge.
c       xmax.  in this case, by(1) = 0.0, the remainder of the
c       by-vector is not calculated, and ncalc is set to
c       min0(nb,0)-2  so that ncalc .ne. nb.
c  ncalc = -1:  y(alpha,x) .ge. xinf.  the requested function
c       values are set to 0.0.
c  1 .lt. ncalc .lt. nb: not all requested function values could
c       be calculated accurately.  by(i) contains correct function
c       values for i .le. ncalc, and and the remaining nb-ncalc
c       array elements contain 0.0.
c
c
c intrinsic functions required are:
c
c     dble, exp, int, max, min, real, sqrt
c
c
c acknowledgement
c
c  this program draws heavily on temme's algol program for y(a,x)
c  and y(a+1,x) and on campbell's programs for y_nu(x).  temme's
c  scheme is used for  x < thresh, and campbell's scheme is used
c  in the asymptotic region.  segments of code from both sources
c  have been translated into fortran 77, merged, and heavily modified.
c  modifications include parameterization of machine dependencies,
c  use of a new approximation for ln(gamma(x)), and built-in
c  protection against over/underflow.
c
c references: "bessel functions j_nu(x) and y_nu(x) of real
c              order and real argument," campbell, j. b.,
c              comp. phy. comm. 18, 1979, pp. 133-142.
c
c             "on the numerical evaluation of the ordinary
c              bessel function of the second kind," temme,
c              n. m., j. comput. phys. 21, 1976, pp. 343-350.
c
c  latest modification: march 19, 1990
c
c  modified by: w. j. cody
c               applied mathematics division
c               argonne national laboratory
c               argonne, il  60439
c
c----------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 I,NB,NCALC
      REAL *8 EX, BY(0:NB), P, SCAL, XINF, XMIN, XLARGE
      REAL *8 U1,U2
c----------------------------------------------------------------------
c  machine-dependent constants
c----------------------------------------------------------------------
      DATA XMIN,XINF/4.46D-308,1.79D308/
      DATA XLARGE/1.0D8/
c----------------------------------------------------------------------
      HALFPI=DATAN(1.0D0)*2.0D0
      EX = X
      IF ((NB.GE.0) .AND. (X.GE.XMIN) .AND. (EX.LT.XLARGE) )  THEN
c
c----------------------------------------------------------------------
c  now have first one or two y's
c----------------------------------------------------------------------
         P=DEXP(-EX)/EX*HALFPI
         BY(0) = P
         BY(1) = P*SCAL*(1.0D0+1.0D0/EX)
         NCALC = 1
         U1=SCAL/EX
         U2=SCAL*SCAL
         DO 400 I = 2, NB
           IF (DABS(BY(I-1))*U1 .GE. XINF/DBLE(2*I-1) ) GO TO 450
           BY(I) =  DBLE(2*I-1)*U1*BY(I-1)+U2*BY(I-2)
           NCALC = NCALC + 1
  400   CONTINUE
  450   DO 460 I = NCALC+1, NB
                BY(I) = ZERO
  460   CONTINUE
      ELSE
        BY(0)=0.0D0
        NCALC=MIN(NB+1,0)-1
      ENDIF
c
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE RIBESL(X,ALPHA,NB,IZE,B,NCALC)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c    this routine calculates bessel functions i sub(n+alpha) (x)
c    for non-negative argument x, and non-negative order n+alpha,
c    with or without exponential scaling.
c
c  on input:
c    x: working precision non-negative real argument for which
c       i's or exponentially scaled i's (i*exp(-x))
c       are to be calculated.  if i's are to be calculated,
c       x must be less than exparg (see below).
c    alpha: working precision fractional part of order for which
c           i's or exponentially scaled i's (i*exp(-x)) are
c           to be calculated.  0 .le. alpha .lt. 1.0.
c    nb: integer number of functions to be calculated, nb .gt. 0.
c        the first function calculated is of order alpha, and the
c        last is of order (nb - 1 + alpha).
c    ize: integer type.  ize = 1 if unscaled i's are to calculated,
c         and 2 if exponentially scaled i's are to be calculated.
c    b: working precision output vector of length nb.  if the routine
c       terminates normally (ncalc=nb), the vector b contains the
c       functions i(alpha,x) through i(nb-1+alpha,x), or the
c       corresponding exponentially scaled functions.
c    ncalc: integer output variable indicating possible errors.
c           before using the vector b, the user should check that
c           ncalc=nb, i.e., all orders have been calculated to
c           the desired accuracy.  see error returns below.
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c explanation of machine-dependent constants
c
c   beta   = radix for the floating-point system
c   minexp = smallest representable power of beta
c   maxexp = smallest power of beta that overflows
c   it     = number of bits in the mantissa of a working precision
c            variable
c   nsig   = decimal significance desired.  should be set to
c            int(log10(2)*it+1).  setting nsig lower will result
c            in decreased accuracy while setting nsig higher will
c            increase cpu time without increasing accuracy.  the
c            truncation error is limited to a relative error of
c            t=.5*10**(-nsig).
c   enten  = 10.0 ** k, where k is the largest integer such that
c            enten is machine-representable in working precision
c   ensig  = 10.0 ** nsig
c   rtnsig = 10.0 ** (-k) for the smallest integer k such that
c            k .ge. nsig/4
c   enmten = smallest abs(x) such that x/4 does not underflow
c   xlarge = upper limit on the magnitude of x when ize=2.  bear
c            in mind that if abs(x)=n, then at least n iterations
c            of the backward recursion will be executed.  the value
c            of 10.0 ** 4 is used on every machine.
c   exparg = largest working precision argument that the library
c            exp routine can handle and upper limit on the
c            magnitude of x when ize=1; approximately
c            log(beta**maxexp)
c
c
c     approximate values for some important machines are:
c
c                        beta       minexp      maxexp       it
c
c  cray-1        (s.p.)    2        -8193        8191        48
c  cyber 180/855
c    under nos   (s.p.)    2         -975        1070        48
c  ieee (ibm/xt,
c    sun, etc.)  (s.p.)    2         -126         128        24
c  ieee (ibm/xt,
c    sun, etc.)  (d.p.)    2        -1022        1024        53
c  ibm 3033      (d.p.)   16          -65          63        14
c  vax           (s.p.)    2         -128         127        24
c  vax d-format  (d.p.)    2         -128         127        56
c  vax g-format  (d.p.)    2        -1024        1023        53
c
c
c                        nsig       enten       ensig      rtnsig
c
c cray-1        (s.p.)    15       1.0e+2465   1.0e+15     1.0e-4
c cyber 180/855
c   under nos   (s.p.)    15       1.0e+322    1.0e+15     1.0e-4
c ieee (ibm/xt,
c   sun, etc.)  (s.p.)     8       1.0e+38     1.0e+8      1.0e-2
c ieee (ibm/xt,
c   sun, etc.)  (d.p.)    16       1.0d+308    1.0d+16     1.0d-4
c ibm 3033      (d.p.)     5       1.0d+75     1.0d+5      1.0d-2
c vax           (s.p.)     8       1.0e+38     1.0e+8      1.0e-2
c vax d-format  (d.p.)    17       1.0d+38     1.0d+17     1.0d-5
c vax g-format  (d.p.)    16       1.0d+307    1.0d+16     1.0d-4
c
c
c                         enmten      xlarge   exparg
c
c cray-1        (s.p.)   1.84e-2466   1.0e+4    5677
c cyber 180/855
c   under nos   (s.p.)   1.25e-293    1.0e+4     741
c ieee (ibm/xt,
c   sun, etc.)  (s.p.)   4.70e-38     1.0e+4      88
c ieee (ibm/xt,
c   sun, etc.)  (d.p.)   8.90d-308    1.0d+4     709
c ibm 3033      (d.p.)   2.16d-78     1.0d+4     174
c vax           (s.p.)   1.17e-38     1.0e+4      88
c vax d-format  (d.p.)   1.17d-38     1.0d+4      88
c vax g-format  (d.p.)   2.22d-308    1.0d+4     709
c
c*******************************************************************
c*******************************************************************
c
c error returns
c
c  in case of an error,  ncalc .ne. nb, and not all i's are
c  calculated to the desired accuracy.
c
c  ncalc .lt. 0:  an argument is out of range. for example,
c     nb .le. 0, ize is not 1 or 2, or ize=1 and abs(x) .ge. exparg.
c     in this case, the b-vector is not calculated, and ncalc is
c     set to min0(nb,0)-1 so that ncalc .ne. nb.
c
c  nb .gt. ncalc .gt. 0: not all requested function values could
c     be calculated accurately.  this usually occurs because nb is
c     much larger than abs(x).  in this case, b(n) is calculated
c     to the desired accuracy for n .le. ncalc, but precision
c     is lost for ncalc .lt. n .le. nb.  if b(n) does not vanish
c     for n .gt. ncalc (because it is too small to be represented),
c     and b(n)/b(ncalc) = 10**(-k), then only the first nsig-k
c     significant figures of b(n) can be trusted.
c
c
c intrinsic functions required are:
c
c     dble, exp, dgamma, gamma, int, max, min, real, sqrt
c
c
c acknowledgement
c
c  this program is based on a program written by david j.
c  sookne (2) that computes values of the bessel functions j or
c  i of real argument and integer order.  modifications include
c  the restriction of the computation to the i bessel function
c  of non-negative real argument, the extension of the computation
c  to arbitrary positive order, the inclusion of optional
c  exponential scaling, and the elimination of most underflow.
c  an earlier version was published in (3).
c
c references: "a note on backward recurrence algorithms," olver,
c              f. w. j., and sookne, d. j., math. comp. 26, 1972,
c              pp 941-947.
c
c             "bessel functions of real argument and integer order,"
c              sookne, d. j., nbs jour. of res. b. 77b, 1973, pp
c              125-132.
c
c             "algorithm 597, sequence of modified bessel functions
c              of the first kind," cody, w. j., trans. math. soft.,
c              1983, pp. 242-245.
c
c  latest modification: may 30, 1989
c
c  modified by: w. j. cody and l. stoltz
c               applied mathematics division
c               argonne national laboratory
c               argonne, il  60439
c
c-------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 IZE,K,L,MAGX,N,NB,NBMX,NCALC,NEND,NSIG,NSTART
      DOUBLE PRECISION DGAMMA,
     1 ALPHA,B,CONST,CONV,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG,
     2 ENTEN,EXPARG,FUNC,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL,
     3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO
      DIMENSION B(NB)
c-------------------------------------------------------------------
c  mathematical constants
c-------------------------------------------------------------------
      DATA ONE,TWO,ZERO,HALF,CONST/1.0D0,2.0D0,0.0D0,0.5D0,1.585D0/
c-------------------------------------------------------------------
c  machine-dependent parameters
c-------------------------------------------------------------------
      DATA NSIG,XLARGE,EXPARG /16,1.0D4,709.0D0/
      DATA ENTEN,ENSIG,RTNSIG/1.0D308,1.0D16,1.0D-4/
      DATA ENMTEN/8.9D-308/
c-------------------------------------------------------------------
c  statement functions for conversion
c-------------------------------------------------------------------
      CONV(N) = DBLE(N)
      FUNC(X) = DGAMMA(X)
c-------------------------------------------------------------------
c check for x, nb, or ize out of range.
c-------------------------------------------------------------------
      IF ((NB.GT.0) .AND. (X .GE. ZERO) .AND.
     1    (ALPHA .GE. ZERO) .AND. (ALPHA .LT. ONE) .AND.
     2    (((IZE .EQ. 1) .AND. (X .LE. EXPARG)) .OR.
     3     ((IZE .EQ. 2) .AND. (X .LE. XLARGE)))) THEN
c-------------------------------------------------------------------
c use 2-term ascending series for small x
c-------------------------------------------------------------------
            NCALC = NB
            MAGX = INT(X)
            IF (X .GE. RTNSIG) THEN
c-------------------------------------------------------------------
c initialize the forward sweep, the p-sequence of olver
c-------------------------------------------------------------------
                  NBMX = NB-MAGX
                  N = MAGX+1
                  EN = CONV(N+N) + (ALPHA+ALPHA)
                  PLAST = ONE
                  P = EN / X
c-------------------------------------------------------------------
c calculate general significance test
c-------------------------------------------------------------------
                  TEST = ENSIG + ENSIG
                  IF (2*MAGX .GT. 5*NSIG) THEN
                        TEST = DSQRT(TEST*P)
                     ELSE
                        TEST = TEST / CONST**MAGX
                  END IF
                  IF (NBMX .GE. 3) THEN
c-------------------------------------------------------------------
c calculate p-sequence until n = nb-1.  check for possible overflow.
c-------------------------------------------------------------------
                     TOVER = ENTEN / ENSIG
                     NSTART = MAGX+2
                     NEND = NB - 1
                     DO 100 K = NSTART, NEND
                        N = K
                        EN = EN + TWO
                        POLD = PLAST
                        PLAST = P
                        P = EN * PLAST/X + POLD
                        IF (P .GT. TOVER) THEN
c-------------------------------------------------------------------
c to avoid overflow, divide p-sequence by tover.  calculate
c p-sequence until abs(p) .gt. 1.
c-------------------------------------------------------------------
                           TOVER = ENTEN
                           P = P / TOVER
                           PLAST = PLAST / TOVER
                           PSAVE = P
                           PSAVEL = PLAST
                           NSTART = N + 1
   60                      N = N + 1
                              EN = EN + TWO
                              POLD = PLAST
                              PLAST = P
                              P = EN * PLAST/X + POLD
                           IF (P .LE. ONE) GO TO 60
                           TEMPB = EN / X
c-------------------------------------------------------------------
c calculate backward test, and find ncalc, the highest n
c such that the test is passed.
c-------------------------------------------------------------------
                           TEST = POLD*PLAST / ENSIG
                           TEST = TEST*(HALF-HALF/(TEMPB*TEMPB))
                           P = PLAST * TOVER
                           N = N - 1
                           EN = EN - TWO
                           NEND = MIN0(NB,N)
                           DO 80 L = NSTART, NEND
                              NCALC = L
                              POLD = PSAVEL
                              PSAVEL = PSAVE
                              PSAVE = EN * PSAVEL/X + POLD
                              IF (PSAVE*PSAVEL .GT. TEST) GO TO 90
   80                      CONTINUE
                           NCALC = NEND + 1
   90                      NCALC = NCALC - 1
                           GO TO 120
                        END IF
  100                CONTINUE
                     N = NEND
                     EN = CONV(N+N) + (ALPHA+ALPHA)
c-------------------------------------------------------------------
c calculate special significance test for nbmx .gt. 2.
c-------------------------------------------------------------------
                     TEST = MAX(TEST,DSQRT(PLAST*ENSIG)*DSQRT(P+P))
                  END IF
c-------------------------------------------------------------------
c calculate p-sequence until significance test passed.
c-------------------------------------------------------------------
  110             N = N + 1
                     EN = EN + TWO
                     POLD = PLAST
                     PLAST = P
                     P = EN * PLAST/X + POLD
                  IF (P .LT. TEST) GO TO 110
c-------------------------------------------------------------------
c initialize the backward recursion and the normalization sum.
c-------------------------------------------------------------------
  120             N = N + 1
                  EN = EN + TWO
                  TEMPB = ZERO
                  TEMPA = ONE / P
                  EM = CONV(N) - ONE
                  EMPAL = EM + ALPHA
                  EMP2AL = (EM - ONE) + (ALPHA + ALPHA)
                  SUM = TEMPA * EMPAL * EMP2AL / EM
                  NEND = N - NB
                  IF (NEND .LT. 0) THEN
c-------------------------------------------------------------------
c n .lt. nb, so store b(n) and set higher orders to zero.
c-------------------------------------------------------------------
                        B(N) = TEMPA
                        NEND = -NEND
                        DO 130 L = 1, NEND
  130                      B(N+L) = ZERO
                     ELSE
                        IF (NEND .GT. 0) THEN
c-------------------------------------------------------------------
c recur backward via difference equation, calculating (but
c not storing) b(n), until n = nb.
c-------------------------------------------------------------------
                           DO 140 L = 1, NEND
                              N = N - 1
                              EN = EN - TWO
                              TEMPC = TEMPB
                              TEMPB = TEMPA
                              TEMPA = (EN*TEMPB) / X + TEMPC
                              EM = EM - ONE
                              EMP2AL = EMP2AL - ONE
                              IF (N .EQ. 1) GO TO 150
                              IF (N .EQ. 2) EMP2AL = ONE
                              EMPAL = EMPAL - ONE
                              SUM = (SUM + TEMPA*EMPAL) * EMP2AL / EM
  140                      CONTINUE
                        END IF
c-------------------------------------------------------------------
c store b(nb)
c-------------------------------------------------------------------
  150                   B(N) = TEMPA
                        IF (NB .LE. 1) THEN
                           SUM = (SUM + SUM) + TEMPA
                           GO TO 230
                        END IF
c-------------------------------------------------------------------
c calculate and store b(nb-1)
c-------------------------------------------------------------------
                        N = N - 1
                        EN = EN - TWO
                        B(N)  = (EN*TEMPA) / X + TEMPB
                        IF (N .EQ. 1) GO TO 220
                        EM = EM - ONE
                        EMP2AL = EMP2AL - ONE
                        IF (N .EQ. 2) EMP2AL = ONE
                        EMPAL = EMPAL - ONE
                        SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM
                  END IF
                  NEND = N - 2
                  IF (NEND .GT. 0) THEN
c-------------------------------------------------------------------
c calculate via difference equation and store b(n), until n = 2.
c-------------------------------------------------------------------
                     DO 200 L = 1, NEND
                        N = N - 1
                        EN = EN - TWO
                        B(N) = (EN*B(N+1)) / X +B(N+2)
                        EM = EM - ONE
                        EMP2AL = EMP2AL - ONE
                        IF (N .EQ. 2) EMP2AL = ONE
                        EMPAL = EMPAL - ONE
                        SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM
  200                CONTINUE
                  END IF
c-------------------------------------------------------------------
c calculate b(1)
c-------------------------------------------------------------------
                  B(1) = TWO*EMPAL*B(2) / X + B(3)
  220             SUM = (SUM + SUM) + B(1)
c-------------------------------------------------------------------
c normalize.  divide all b(n) by sum.
c-------------------------------------------------------------------
  230             IF (ALPHA .NE. ZERO)
     1               SUM = SUM * FUNC(ONE+ALPHA) * (X*HALF)**(-ALPHA)
                  IF (IZE .EQ. 1) SUM = SUM * DEXP(-X)
                  TEMPA = ENMTEN
                  IF (SUM .GT. ONE) TEMPA = TEMPA * SUM
                  DO 260 N = 1, NB
                     IF (B(N) .LT. TEMPA) B(N) = ZERO
                     B(N) = B(N) / SUM
  260             CONTINUE
                  RETURN
c-------------------------------------------------------------------
c two-term ascending series for small x.
c-------------------------------------------------------------------
               ELSE
                  TEMPA = ONE
                  EMPAL = ONE + ALPHA
                  HALFX = ZERO
                  IF (X .GT. ENMTEN) HALFX = HALF * X
                  IF (ALPHA .NE. ZERO) TEMPA = HALFX**ALPHA /FUNC(EMPAL)
                  IF (IZE .EQ. 2) TEMPA = TEMPA * DEXP(-X)
                  TEMPB = ZERO
                  IF ((X+ONE) .GT. ONE) TEMPB = HALFX * HALFX
                  B(1) = TEMPA + TEMPA*TEMPB / EMPAL
                  IF ((X .NE. ZERO) .AND. (B(1) .EQ. ZERO)) NCALC = 0
                  IF (NB .GT. 1) THEN
                     IF (X .EQ. ZERO) THEN
                           DO 310 N = 2, NB
                              B(N) = ZERO
  310                      CONTINUE
                        ELSE
c-------------------------------------------------------------------
c calculate higher-order functions.
c-------------------------------------------------------------------
                           TEMPC = HALFX
                           TOVER = (ENMTEN + ENMTEN) / X
                           IF (TEMPB .NE. ZERO) TOVER = ENMTEN / TEMPB
                           DO 340 N = 2, NB
                              TEMPA = TEMPA / EMPAL
                              EMPAL = EMPAL + ONE
                              TEMPA = TEMPA * TEMPC
                              IF (TEMPA .LE. TOVER*EMPAL) TEMPA = ZERO
                              B(N) = TEMPA + TEMPA*TEMPB / EMPAL
                              IF ((B(N) .EQ. ZERO) .AND. (NCALC .GT. N))
     1                             NCALC = N-1
  340                      CONTINUE
                     END IF
                  END IF
            END IF
         ELSE
            NCALC = MIN0(NB,0)-1
      END IF
      RETURN
c---------- last line of ribesl ----------
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      DOUBLE PRECISION FUNCTION DGAMMA(X)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c this routine calculates the gamma function for a real argument x.
c   computation is based on an algorithm outlined in reference 1.
c   the program uses rational functions that approximate the gamma
c   function to at least 20 significant decimal digits.  coefficients
c   for the approximation over the interval (1,2) are unpublished.
c   those for the approximation for x .ge. 12 are from reference 2.
c   the accuracy achieved depends on the arithmetic system, the
c   compiler, the intrinsic functions, and proper selection of the
c   machine-dependent constants.
c
c
c*******************************************************************
c*******************************************************************
c
c explanation of machine-dependent constants
c
c beta   - radix for the floating-point representation
c maxexp - the smallest positive power of beta that overflows
c xbig   - the largest argument for which gamma(x) is representable
c          in the machine, i.e., the solution to the equation
c                  gamma(xbig) = beta**maxexp
c xinf   - the largest machine representable floating-point number;
c          approximately beta**maxexp
c eps    - the smallest positive floating-point number such that
c          1.0+eps .gt. 1.0
c xminin - the smallest positive floating-point number such that
c          1/xminin is machine representable
c
c     approximate values for some important machines are:
c
c                            beta       maxexp        xbig
c
c cray-1         (s.p.)        2         8191        966.961
c cyber 180/855
c   under nos    (s.p.)        2         1070        177.803
c ieee (ibm/xt,
c   sun, etc.)   (s.p.)        2          128        35.040
c ieee (ibm/xt,
c   sun, etc.)   (d.p.)        2         1024        171.624
c ibm 3033       (d.p.)       16           63        57.574
c vax d-format   (d.p.)        2          127        34.844
c vax g-format   (d.p.)        2         1023        171.489
c
c                            xinf         eps        xminin
c
c cray-1         (s.p.)   5.45e+2465   7.11e-15    1.84e-2466
c cyber 180/855
c   under nos    (s.p.)   1.26e+322    3.55e-15    3.14e-294
c ieee (ibm/xt,
c   sun, etc.)   (s.p.)   3.40e+38     1.19e-7     1.18e-38
c ieee (ibm/xt,
c   sun, etc.)   (d.p.)   1.79d+308    2.22d-16    2.23d-308
c ibm 3033       (d.p.)   7.23d+75     2.22d-16    1.39d-76
c vax d-format   (d.p.)   1.70d+38     1.39d-17    5.88d-39
c vax g-format   (d.p.)   8.98d+307    1.11d-16    1.12d-308
c
c*******************************************************************
c*******************************************************************
c
c error returns
c
c  the program returns the value xinf for singularities or
c     when overflow would occur.  the computation is believed
c     to be free of underflow and overflow.
c
c
c  intrinsic functions required are:
c
c     int, dble, exp, log, real, sin
c
c
c references: "an overview of software development for special
c              functions", w. j. cody, lecture notes in mathematics,
c              506, numerical analysis dundee, 1975, g. a. watson
c              (ed.), springer verlag, berlin, 1976.
c
c              computer approximations, hart, et. al., wiley and
c              sons, new york, 1968.
c
c  latest modification: october 12, 1989
c
c  authors: w. j. cody and l. stoltz
c           applied mathematics division
c           argonne national laboratory
c           argonne, il 60439
c
c----------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 I,N
      LOGICAL PARITY
      DOUBLE PRECISION
     1    C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE,
     2    TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
      DIMENSION C(7),P(8),Q(8)
c----------------------------------------------------------------------
c  mathematical constants
c----------------------------------------------------------------------
      DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
     1     SQRTPI/0.9189385332046727417803297D0/,
     2     PI/3.1415926535897932384626434D0/
c----------------------------------------------------------------------
c  machine dependent parameters
c----------------------------------------------------------------------
      DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
     1     XINF/1.79D308/
c----------------------------------------------------------------------
c  numerator and denominator coefficients for rational minimax
c     approximation over (1,2).
c----------------------------------------------------------------------
      DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
     1       -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
     2       8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
     3       -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
      DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
     1      -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
     2        2.25381184209801510330112D+4,4.75584627752788110767815D+3,
     3      -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
c----------------------------------------------------------------------
c  coefficients for minimax approximation over (12, inf).
c----------------------------------------------------------------------
      DATA C/-1.910444077728D-03,8.4171387781295D-04,
     1     -5.952379913043012D-04,7.93650793500350248D-04,
     2     -2.777777777777681622553D-03,8.333333333333333331554247D-02,
     3      5.7083835261D-03/
c----------------------------------------------------------------------
c  statement functions for conversion between integer and float
c----------------------------------------------------------------------
      CONV(I) = DBLE(I)
      PARITY = .FALSE.
      FACT = ONE
      N = 0
      Y = X
      IF (Y .LE. ZERO) THEN
c----------------------------------------------------------------------
c  argument is negative
c----------------------------------------------------------------------
            Y = -X
            Y1 = DINT(Y)
            RES = Y - Y1
            IF (RES .NE. ZERO) THEN
                  IF (Y1 .NE. DINT(Y1*HALF)*TWO) PARITY = .TRUE.
                  FACT = -PI / SIN(PI*RES)
                  Y = Y + ONE
               ELSE
                  RES = XINF
                  GO TO 900
            END IF
      END IF
c----------------------------------------------------------------------
c  argument is positive
c----------------------------------------------------------------------
      IF (Y .LT. EPS) THEN
c----------------------------------------------------------------------
c  argument .lt. eps
c----------------------------------------------------------------------
            IF (Y .GE. XMININ) THEN
                  RES = ONE / Y
               ELSE
                  RES = XINF
                  GO TO 900
            END IF
         ELSE IF (Y .LT. TWELVE) THEN
            Y1 = Y
            IF (Y .LT. ONE) THEN
c----------------------------------------------------------------------
c  0.0 .lt. argument .lt. 1.0
c----------------------------------------------------------------------
                  Z = Y
                  Y = Y + ONE
               ELSE
c----------------------------------------------------------------------
c  1.0 .lt. argument .lt. 12.0, reduce argument if necessary
c----------------------------------------------------------------------
                  N = INT(Y) - 1
                  Y = Y - CONV(N)
                  Z = Y - ONE
            END IF
c----------------------------------------------------------------------
c  evaluate approximation for 1.0 .lt. argument .lt. 2.0
c----------------------------------------------------------------------
            XNUM = ZERO
            XDEN = ONE
            DO 260 I = 1, 8
               XNUM = (XNUM + P(I)) * Z
               XDEN = XDEN * Z + Q(I)
  260       CONTINUE
            RES = XNUM / XDEN + ONE
            IF (Y1 .LT. Y) THEN
c----------------------------------------------------------------------
c  adjust result for case  0.0 .lt. argument .lt. 1.0
c----------------------------------------------------------------------
                  RES = RES / Y1
               ELSE IF (Y1 .GT. Y) THEN
c----------------------------------------------------------------------
c  adjust result for case  2.0 .lt. argument .lt. 12.0
c----------------------------------------------------------------------
                  DO 290 I = 1, N
                     RES = RES * Y
                     Y = Y + ONE
  290             CONTINUE
            END IF
         ELSE
c----------------------------------------------------------------------
c  evaluate for argument .ge. 12.0,
c----------------------------------------------------------------------
            IF (Y .LE. XBIG) THEN
                  YSQ = Y * Y
                  SUM = C(7)
                  DO 350 I = 1, 6
                     SUM = SUM / YSQ + C(I)
  350             CONTINUE
                  SUM = SUM/Y - Y + SQRTPI
                  SUM = SUM + (Y-HALF)*DLOG(Y)
                  RES = DEXP(SUM)
               ELSE
                  RES = XINF
                  GO TO 900
            END IF
      END IF
c----------------------------------------------------------------------
c  final adjustments and return
c----------------------------------------------------------------------
      IF (PARITY) RES = -RES
      IF (FACT .NE. ONE) RES = FACT / RES
  900 DGAMMA = RES
      RETURN
c ---------- last line of gamma ----------
      END
c
c**********************************************************************
      SUBROUTINE LGNDRGT1(SCAL,NMAX, X, Y)
c**********************************************************************
c
c-----this subroutine calculates the legendre function for x >1.
c     the scaled version. scal^n* p_n^m
c     scal should be less than 1.
c
c     input :
c       scal : the scaling factor.
c       nmax : the number of terms we want to compute.
c       x : the value at which we want to compute.
c
c     output :
c       y() : the scal^n* p_n^m.
c
c***********************************************************************
      IMPLICIT REAL *8 (A-H,O-Z)
      INTEGER *4 NMAX, M
      REAL *8  SCAL, X, Y(0:NMAX,0:NMAX), U, V, W
c
      V=SCAL*X
      W=SCAL*SCAL
      U=DSQRT(X*X-1)*SCAL
      Y(0,0)=1.0D0
      DO 10 M=0, NMAX
         IF (M.GT.0)  Y(M,M)=Y(M-1,M-1)*U*DBLE(2*M-1)
         IF (M.LT.NMAX)  Y(M+1,M)=DBLE(2*M+1)*V*Y(M,M)
         DO 20 N=M+2, NMAX
            Y(N,M)=((2.0D0*DBLE(N)-1.0D0)*V*Y(N-1,M)-DBLE(N+M-1)
     1         *W*Y(N-2,M)) / DBLE(N-M)
 20      CONTINUE
 10   CONTINUE
      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 YDNFORMMP(BETA,X0Y0Z0,ZPARTS,CHARGE,DN,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. this subroutine is for the calculation of
c      double layer potential.
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    dn(3,nparts): the normal direction.
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 DN(3,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,CPZ2
      COMPLEX *16 CDN1,CDN2
      DATA IMAG/(0.0D0,1.0D0)/
      DATA PRECIS/1.0D-14/
c
c-----first form the multipole expansion.
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
        CDN1=DN(1,I)-DN(2,I)*IMAG
        CDN2=DCONJG(CDN1)
c
c-------contribution from (0,0) term.
c
        CP = CHARGE(I)*BI(0)*BETA/SCALE
c
        MPOLE(1,1)=MPOLE(1,1)+(CP/2.0D0)*CDN1
        MPOLE(1,0)=MPOLE(1,0)-CP*DN(3,I)
c
c-------compute legendre polynomials of argument cos(theta) = ctheta
c         and add contributions from legendre polynomials
c
        CALL LGNDR(NTERMS,CTHETA,P)
c
c-------contribution from (1,0) term.
c
        CP=CHARGE(I)*P(1,0)*BI(1)*C(1,0)*BETA/3.0D0
c
        MPOLE(2,1)=MPOLE(2,1)+(CP/2.0D0/SCALE)*CDN1
        MPOLE(0,0)=MPOLE(0,0)-CP*SCALE*DN(3,I)
        MPOLE(2,0)=MPOLE(2,0)-CP/SCALE*2D0*DN(3,I)
c
c-------contribution from (l,0) terms, nterms>l>1.
c
        DO 4300 L = 2,NTERMS-1
          CP = CHARGE(I)*P(L,0)*BI(L)*C(L,0)*BETA/DBLE(2*L+1)
          CPZ= (CP/2.0D0)*CDN1
          CP=CP*DN(3,I)
c
          MPOLE(L-1,1)=MPOLE(L-1,1)-CPZ*SCALE
          MPOLE(L+1,1)=MPOLE(L+1,1)+CPZ/SCALE
          MPOLE(L-1,0)=MPOLE(L-1,0)-CP*SCALE*DBLE(L)
          MPOLE(L+1,0)=MPOLE(L+1,0)-CP/SCALE*DBLE(L+1)
4300    CONTINUE
c
c-------contribution from (nterms,0).
c
        L=NTERMS
        CP = CHARGE(I)*P(L,0)*BI(L)*C(L,0)*BETA/DBLE(2*L+1)*SCALE
c
        MPOLE(L-1,1)=MPOLE(L-1,1)-(CP/2.0D0)*CDN1
        MPOLE(L-1,0)=MPOLE(L-1,0)-CP*DBLE(L)*DN(3,I)
c
c-------contribution from (1,1) and (1,-1) terms.
c
        CP= CHARGE(I)*BI(1)*C(1,1)*P(1,1)*BETA/3.0D0
        CPZ= CP*EPHI(1)
        CP=2D0*DREAL(CPZ*CDN2)
        MPOLE(0,0)=MPOLE(0,0)+CP*SCALE
        MPOLE(2,0)=MPOLE(2,0)-CP/SCALE
        CPZ=CPZ/SCALE
        MPOLE(2,2)=MPOLE(2,2)+(CPZ/2.0D0)*CDN1
        MPOLE(2,1)=MPOLE(2,1)-CPZ*DN(3,I)
c
        DO 4500 L = 2,NTERMS-1
          DO 4400 M=1,L-2
            CP = CHARGE(I)*BI(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)
            CPZ = CP*EPHI(M)
c
c-----------contributions from l=2:nterms-1, m=1,l-2
c
            CPZ2= CPZ*CDN2
            IF (M.EQ.1) THEN
              CP=DBLE(L*(L+1))*DREAL(CPZ2)
              MPOLE(L-1,0)=MPOLE(L-1,0)+CP*SCALE
              MPOLE(L+1,0)=MPOLE(L+1,0)-CP/SCALE
            ELSE
              CPZ2=CPZ2/2.0D0
              MPOLE(L-1,M-1)=MPOLE(L-1,M-1)+CPZ2*
     1          (DBLE((L+M-1)*(L+M))*SCALE)
              MPOLE(L+1,M-1)=MPOLE(L+1,M-1)-CPZ2*
     1          (DBLE((L-M+1)*(L-M+2))/SCALE)
            ENDIF
c
            CPZ2=CPZ*CDN1/2.0D0
            MPOLE(L-1,M+1)=MPOLE(L-1,M+1)-CPZ2*SCALE
            MPOLE(L+1,M+1)=MPOLE(L+1,M+1)+CPZ2/SCALE
            CPZ2=CPZ*DN(3,I)
            MPOLE(L-1,M)=MPOLE(L-1,M)-CPZ2*(DBLE(L+M)*SCALE)
            MPOLE(L+1,M)=MPOLE(L+1,M)-CPZ2*(DBLE(L-M+1)/SCALE)
4400      CONTINUE
c
c---------m=l-1
c
          M=L-1
          CP = CHARGE(I)*BI(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)
          CPZ = CP*EPHI(M)
c
          CPZ2=CPZ*CDN2
          IF (M.EQ.1) THEN
            CP=DREAL(CPZ2)
            MPOLE(L-1,0)=MPOLE(L-1,0)+CP*DBLE((L)*(L+1))*SCALE
            MPOLE(L+1,0)=MPOLE(L+1,0)-CP*6.0D0/SCALE
          ELSE
            CPZ2=CPZ2/2D0
            MPOLE(L-1,M-1)=MPOLE(L-1,M-1)+
     1        CPZ2*(DBLE((L+M-1)*(L+M))*SCALE)
            MPOLE(L+1,M-1)=MPOLE(L+1,M-1)-
     1        CPZ2*(6.0D0/SCALE)
          ENDIF
c
          CPZ2=CPZ*CDN1/2D0
          MPOLE(L+1,M+1)=MPOLE(L+1,M+1)+CPZ2/SCALE
          CPZ2=CPZ*DN(3,I)
          MPOLE(L-1,M)=MPOLE(L-1,M)-CPZ2*(DBLE(L+M)*SCALE)
          MPOLE(L+1,M)=MPOLE(L+1,M)-CPZ2*(2D0/SCALE)
c
c---------m=l
c
          M=L
          CP = CHARGE(I)*BI(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)
          CPZ = CP*EPHI(M)
c
          CPZ2=CPZ/2.0D0*CDN2
          MPOLE(L-1,M-1)=MPOLE(L-1,M-1)+
     1      CPZ2*(DBLE((L+M-1)*(L+M))*SCALE)
          MPOLE(L+1,M-1)=MPOLE(L+1,M-1)-CPZ2*(2.0D0/SCALE)
          MPOLE(L+1,M+1)=MPOLE(L+1,M+1)+CPZ/(2.0D0*SCALE)*CDN1
          MPOLE(L+1,M)=MPOLE(L+1,M)-CPZ*(DN(3,I)/SCALE)
4500    CONTINUE
c
c-------contribution from l=nterms, m=1,l-1
c
        L=NTERMS
c
c-------contribution from (nterms, 1)
c
        DO 5400 M=1,L-2
          CP=CHARGE(I)*BI(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)*SCALE
          CPZ = CP*EPHI(M)
c
          CPZ2=CPZ*CDN2
          IF (M.EQ.1) THEN
            CP=DREAL(CPZ2)
            MPOLE(L-1,0)=MPOLE(L-1,0)+CP*DBLE(L*(L+1))
          ELSE
            MPOLE(L-1,M-1)=MPOLE(L-1,M-1)+CPZ2
     1        *(DBLE((L+M-1)*(L+M))/2.0D0)
          ENDIF
c
          MPOLE(L-1,M+1)=MPOLE(L-1,M+1)-CPZ/2.0D0*CDN1
          MPOLE(L-1,M)=MPOLE(L-1,M)-CPZ*(DBLE(L+M)*DN(3,I))
5400    CONTINUE
c
c-------m=l-1
c
        M=L-1
        CP=CHARGE(I)*BI(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)*SCALE
        CPZ=CP*EPHI(M)
c
        MPOLE(L-1,M-1)=MPOLE(L-1,M-1)+CPZ
     1    *(DBLE((L+M-1)*(L+M))/2.0D0)*CDN2
        MPOLE(L-1,M)=MPOLE(L-1,M)-CPZ*(DBLE(L+M)*DN(3,I))
c
c-------m=l
c
        M=L
        CP = CHARGE(I)*BI(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)*SCALE
        CPZ = CP*EPHI(M)
c
        MPOLE(L-1,M-1)=MPOLE(L-1,M-1)+CPZ
     1    *(DBLE((L+M-1)*(L+M))/2.0D0)*CDN2
c
4900  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YDNLC(BETA,X0Y0Z0,ZPARTS,CHARGE,DN,NPARTS,
     1  LOCAL,NTERMS,SCALE,P,C)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine forms the local expansion cause by the nparts
c      particles in the box. this subroutine is for the calculation of
c      double layer potential.
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    dn(3,nparts): the normal direction.
c    nparts: the total number of particles.
c    nterms: order of desired expansion
c    scale: the scaling factor.
c
c  on output:
c
c    local: 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 DN(3,NPARTS)
      REAL *8 P(0:NTERMS,0:NTERMS),C(0:NTERMS,0:NTERMS)
      COMPLEX *16 LOCAL(0:NTERMS,0:NTERMS)
c
c-----local variables.
c
      INTEGER *4 I,L,M,NCALC
      REAL *8 RK,BK(0:60)
      REAL *8 PRECIS,D,CP,PROJ
      REAL *8 RX,RY,RZ,RR,CTHETA
      COMPLEX *16 IMAG,EPHI(1:60),CPZ,CPZ2
      DATA IMAG/(0.0D0,1.0D0)/
      DATA PRECIS/1.0D-14/
c
c-----first form the local expansion.
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 KN(SCALE,RK,NTERMS,BK,NCALC)
        IF (NCALC .NE. NTERMS) 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 expansion.'
          STOP
        ENDIF
c
c-------contribution from (0,0) term.
c
        CP = CHARGE(I)*BK(0)*BETA*SCALE
        LOCAL(1,1)=LOCAL(1,1)-CP/2.0D0*(DN(1,I)-DN(2,I)*IMAG)
        LOCAL(1,0)=LOCAL(1,0)+CP*DN(3,I)
c
c-------compute legendre polynomials of argument cos(theta) = ctheta
c         and add contributions from legendre polynomials
c
        CALL LGNDR(NTERMS,CTHETA,P)
c
c-------contribution from (1,0) term.
c
        CP=CHARGE(I)*P(1,0)*BK(1)*C(1,0)*BETA/3.0D0
c
        LOCAL(2,1)=LOCAL(2,1)-CP/2.0D0*SCALE*
     1    (DN(1,I)-DN(2,I)*IMAG)
        LOCAL(0,0)=LOCAL(0,0)+CP/SCALE*DN(3,I)
        LOCAL(2,0)=LOCAL(2,0)+CP*SCALE*2D0*DN(3,I)
c
c-------contribution from (l,0) terms, nterms>l>1.
c
        DO 4300 L = 2,NTERMS-1
          CP = CHARGE(I)*P(L,0)*BK(L)*C(L,0)*BETA/DBLE(2*L+1)
          CPZ= CP/2D0*(DN(1,I)-DN(2,I)*IMAG)
          CP = CP*DN(3,I)
c
          LOCAL(L-1,1)=LOCAL(L-1,1)+CPZ/SCALE
          LOCAL(L+1,1)=LOCAL(L+1,1)-CPZ*SCALE
          LOCAL(L-1,0)=LOCAL(L-1,0)+CP/SCALE*DBLE(L)
          LOCAL(L+1,0)=LOCAL(L+1,0)+CP*SCALE*DBLE(L+1)
4300    CONTINUE
c
c-------contribution from (nterms,0).
c
        L=NTERMS
        CP=CHARGE(I)*P(L,0)*BK(L)*C(L,0)*BETA/DBLE(2*L+1)/SCALE
c
        LOCAL(L-1,1)=LOCAL(L-1,1)+CP/2D0*(DN(1,I)-DN(2,I)*IMAG)
        LOCAL(L-1,0)=LOCAL(L-1,0)+CP*DBLE(L)*DN(3,I)
c
c-------contribution from (1,1) term.
c
        CP=CHARGE(I)*BK(1)*C(1,1)*P(1,1)*BETA/3D0
        CPZ=CP*EPHI(1)
        CP=2D0*DREAL(CPZ*(DN(1,I)+DN(2,I)*IMAG))
        LOCAL(0,0)=LOCAL(0,0)-CP/SCALE
        LOCAL(2,0)=LOCAL(2,0)+CP*SCALE
        CPZ=CPZ*SCALE
        LOCAL(2,2)=LOCAL(2,2)-CPZ/2D0*(DN(1,I)-DN(2,I)*IMAG)
        LOCAL(2,1)=LOCAL(2,1)+CPZ*DN(3,I)
c
        DO 4500 L = 2,NTERMS-1
          DO 4400 M=1,L-2
            CP=CHARGE(I)*BK(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)
            CPZ=CP*EPHI(M)
c
c-----------contributions from l=2:nterms-1, m=1,l-2
c
            CPZ2=CPZ*(DN(1,I)+IMAG*DN(2,I))
            IF (M.EQ.1) THEN
              CP=DBLE(L*(L+1))*DREAL(CPZ2)
              LOCAL(L-1,0)=LOCAL(L-1,0)-CP/SCALE
              LOCAL(L+1,0)=LOCAL(L+1,0)+CP*SCALE
            ELSE
              CPZ2=CPZ2/2D0
              LOCAL(L-1,M-1)=LOCAL(L-1,M-1)-CPZ2*
     1          DBLE((L+M-1)*(L+M))/SCALE
              LOCAL(L+1,M-1)=LOCAL(L+1,M-1)+CPZ2*
     1          DBLE((L-M+1)*(L-M+2))*SCALE
            ENDIF
c
            CPZ2=CPZ*(DN(1,I)-IMAG*DN(2,I))/2.0D0
            LOCAL(L-1,M+1)=LOCAL(L-1,M+1)+CPZ2/SCALE
            LOCAL(L+1,M+1)=LOCAL(L+1,M+1)-CPZ2*SCALE
            CPZ2=CPZ*DN(3,I)
            LOCAL(L-1,M)=LOCAL(L-1,M)+CPZ2*DBLE(L+M)/SCALE
            LOCAL(L+1,M)=LOCAL(L+1,M)+CPZ2*DBLE(L-M+1)*SCALE
4400      CONTINUE
c
c---------m=l-1
c
          M=L-1
          CP=CHARGE(I)*BK(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)
          CPZ=CP*EPHI(M)
c
          CPZ2=CPZ*(DN(1,I)+IMAG*DN(2,I))
          IF (M.EQ.1) THEN
            CP=DREAL(CPZ2)
            LOCAL(L-1,0)=LOCAL(L-1,0)-CP*DBLE((L)*(L+1))/SCALE
            LOCAL(L+1,0)=LOCAL(L+1,0)+CP*6.0D0*SCALE
          ELSE
            CPZ2=CPZ2/2D0
            LOCAL(L-1,M-1)=LOCAL(L-1,M-1)-CPZ2*
     1        DBLE((L+M-1)*(L+M))/SCALE
            LOCAL(L+1,M-1)=LOCAL(L+1,M-1)+CPZ2*
     1        6.0D0*SCALE
          ENDIF
c
          CPZ2=CPZ*(DN(1,I)-IMAG*DN(2,I))/2D0
          LOCAL(L+1,M+1)=LOCAL(L+1,M+1)-CPZ2*SCALE
          CPZ2=CPZ*DN(3,I)
          LOCAL(L-1,M)=LOCAL(L-1,M)+CPZ2*DBLE(L+M)/SCALE
          LOCAL(L+1,M)=LOCAL(L+1,M)+CPZ2*2.0D0*SCALE
c
c---------m=l
c
          M=L
          CP=CHARGE(I)*BK(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)
          CPZ=CP*EPHI(M)
c
          CPZ2=CPZ/2.0D0*(DN(1,I)+IMAG*DN(2,I))
          LOCAL(L-1,M-1)=LOCAL(L-1,M-1)-CPZ2*DBLE((L+M-1)*(L+M))/SCALE
          LOCAL(L+1,M-1)=LOCAL(L+1,M-1)+CPZ2*2D0*SCALE
          LOCAL(L+1,M+1)=LOCAL(L+1,M+1)-CPZ/2.0D0
     1      *SCALE*(DN(1,I)-IMAG*DN(2,I))
          LOCAL(L+1,M)=LOCAL(L+1,M)+CPZ*SCALE*DN(3,I)
4500    CONTINUE
c
c-------contribution from l=nterms, m=1,l-1
c
        L=NTERMS
c
c-------contribution from (nterms, 1)
c
        DO 5400 M=1,L-2
          CP=CHARGE(I)*BK(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)/SCALE
          CPZ=CP*EPHI(M)
c
          CPZ2=CPZ*(DN(1,I)+IMAG*DN(2,I))
          IF (M.EQ.1) THEN
            CP=DREAL(CPZ2)
            LOCAL(L-1,0)=LOCAL(L-1,0)-CP*DBLE(L*(L+1))
          ELSE
            LOCAL(L-1,M-1)=LOCAL(L-1,M-1)-CPZ2/2.0D0
     1        *DBLE((L+M-1)*(L+M))
          ENDIF
c
          LOCAL(L-1,M+1)=LOCAL(L-1,M+1)+CPZ/2.0D0
     1      *(DN(1,I)-IMAG*DN(2,I))
          LOCAL(L-1,M)=LOCAL(L-1,M)+CPZ*DBLE(L+M)*DN(3,I)
5400    CONTINUE
c
c-------m=l-1
c
        M=L-1
        CP=CHARGE(I)*BK(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)/SCALE
        CPZ=CP*EPHI(M)
c
        LOCAL(L-1,M-1)=LOCAL(L-1,M-1)-CPZ/2.0D0
     1    *DBLE((L+M-1)*(L+M))*(DN(1,I)+IMAG*DN(2,I))
        LOCAL(L-1,M)=LOCAL(L-1,M)+CPZ*DBLE(L+M)*DN(3,I)
c
c-------m=l
c
        M=L
        CP=CHARGE(I)*BK(L)*C(L,M)*P(L,M)*BETA/DBLE(2*L+1)/SCALE
        CPZ=CP*EPHI(M)
c
        LOCAL(L-1,M-1)=LOCAL(L-1,M-1)-CPZ/2.0D0
     1    *DBLE((L+M-1)*(L+M))*(DN(1,I)+IMAG*DN(2,I))
c
4900  CONTINUE
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 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
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YDNDIRECI(BETA,NPARTS,ZPARTS,CHARGE,DN,I,
     1                   RPOT,FIELD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    calculates directly the potential for particle i. this is for
c      the double layer potential, so the normal direction for each
c      particle is specified in dn().
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    dn(3,:) : the normal direction for each 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),DN(3,NPARTS)
      REAL *8 RPOT,FIELD(3)
c
      INTEGER *4 J
      REAL *8 PI,RX,RY,RZ,RR,RDIS,EXPR,TERM1,TERM2,TERM3,TERM4
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
        TERM2 = -TERM1*(1.0D0+RDIS)/RR
        TERM3 = RX*DN(1,J)+RY*DN(2,J)+RZ*DN(3,J)
        RPOT = RPOT + TERM2*TERM3
c
c-------the field. see the mathematica file "yukawaforce" for
c         formula details.
c
        TERM4=TERM3*(3.0D0+3.0D0*RDIS+RDIS*RDIS)/RR/RR
        FIELD(1)=FIELD(1)+TERM2*DN(1,J)+TERM1*TERM4*RX
        FIELD(2)=FIELD(2)+TERM2*DN(2,J)+TERM1*TERM4*RY
        FIELD(3)=FIELD(3)+TERM2*DN(3,J)+TERM1*TERM4*RZ
1000  CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE YDIRECIST(BETA,NSOU,ZSOU,CHARGE,I,ZTAR,
     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    zsou(3,:) : the location of the source particles.
c    charge : the charge of the particle.
c    i : the index of the target particle.
c    ztar: the location of the target particles.
c
c  on output :
c    rpot : the potential.
c    field: the field at target points.
c
c  subroutine called : dsqrt, dexp,
c
c  called from : main()
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      IMPLICIT NONE
c
      INTEGER *4 NSOU,I
      REAL *8 BETA,ZSOU(3,NSOU),CHARGE(NSOU),ZTAR(3,1)
      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,NSOU
        RX = ZTAR(1,I) - ZSOU(1,J)
        RY = ZTAR(2,I) - ZSOU(2,J)
        RZ = ZTAR(3,I) - ZSOU(3,J)
        RR = RX*RX + RY*RY + RZ*RZ
c
c-------stop if two particles are too close.
c
        IF (RR.LE.1D-60) THEN
          PRINT *, 'source and target toooo close, stop'
          STOP
        ENDIF
c
        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
      RETURN
      END
c
