SUBROUTINE ESELF(CMETHD,CXZP,NX,NY,NZ,AKD,CXZC,CXZW,CXZE)
C Arguments:
      CHARACTER*6 CMETHD
      REAL AKD
      INTEGER NX,NY,NZ
      COMPLEX CXZC(NX+1,NY+1,NZ+1,6),CXZE(NX,NY,NZ,3),
     &        CXZP(NX,NY,NZ,3),CXZW(2*NX,2*NY,2*NZ,*)
C Local scalars:
      COMPLEX CXEX,CXEY,CXEZ,CXFAC,CXPHAS,CXUNIT,CXXX,CXXY,CXXZ,
     &        CXYY,CXYZ,CXZERO,CXZZ
      REAL FAC,R,R2,R3,WOLD,XX
      INTEGER I,IER,II,IR,ISGN,J,JJ,JR,JSGN,K,KR,KSGN,M,NGRID
C Local arrays:
      REAL X(3)
      INTEGER ISYM(3),NN(3)
C
      EXTERNAL C3DFFT,EXTND,FOURX,PAD,TRIM
      INTRINSIC EXP,MIN,NINT,SIGN,SQRT
      SAVE CXUNIT,CXZERO,NN,NGRID,WOLD
***********************************************************************
C
C       Given the dipole moments, CXZP, at all points on
C       a rectangular grid,oscillating at frequency AKD,
C       compute the electric field amplitude, CXZE,
C       at each point produced by all the other dipoles except the one
C       at that point.
C       The relationship between the dipoles and the field
C       values at the grid points can be expressed as a convolution,
C       and the convolution is efficiently evaluated using 3D FFTs.
C
C       3-dimensional FFT is performed using either
C       a) FOURX code by  N. Brenner + slight PJF mods
C          Works with either MXMEM=0 or 1 in DDSCAT.f
C          Select CMETHD='BRENNR' (in the DDSCAT parameter file)
C       b) CXFFT3 interface to Temperton's CFFT99F code.
C          This option is more suitable for vector machines
C          (CRAY) but takes (8*nx*ny*nz) more memory in comparison
C          with FOURX.  It also requires that NX,NY,NZ each be of form
C          (2**I)*(3**J)*(5**K), for integer I,J,K.  Subroutine EXTEND
C          takes care of choosing suitable NX,NY,NZ.
C          *** Requires that DDSCAT be compiled with MXMEM=1 ***
C          Select CMETHD='TMPRTN' (in the DDSCAT parameter file).
C       c) Native 3d FFT routine in Convex vector library.
C          This option is available only on Convex hardware, and
C          linking must be done with -lveclib option (see the Makefile).
C          This also requires that NX,NY,NZ be of form
C          (2**I)*(3**J)*(5**K); subroutine EXTEND takes care of choosing
C          suitable NX,NY,NZ.
C          Works with either MXMEM=0 or 1 in DDSCAT.f
C          Select CMETHD='CONVEX' (in the DDSCAT parameter file).
C       d) CXFFT3N interface to GPFA code of Temperton. 
C          Good points: 
C             -On CRAY the code is on average 30-40% faster in comparison
C              to TMPRTN (and 10-15 faster in comparison to BRENNR)
C             -On scalar machines the code is 2-8 faster in comparison
C              to BRENNR or TMPRTN
C             -Doesn't require additional storage
C          Problems:
C              -Requires that NX,NY,NZ be of form (2**I)*(3**J)*(5**K); 
C               subroutine EXTEND takes care of choosing suitable NX,NY,NZ.
C              -The code is NOT ON PUBLIC DOMAIN and we do not distribute it.
C              -The choice of "lvr" variable in gpfa2f, gpfa3f, gpfa5f
C               depends on machine. WARNING: on C90 use lvr=128, on
C               all other CRAY's use lvr=64; wrong lvr will produce WRONG
C               RESULTS. On scalar machines lvr depends on cache length
C               but the bad choice degrades performance but still produces
C               correct results.
C        
C       INPUT:
C       CXZP(I,J,K,L)   Lth cartesian component of the dipole
C                       moment at the grid point (I,J,K);
C                       the DIMENSIONed length of CXZP in
C                       the calling routine is CXZP(NX,NY,NZ,3)
C                       [or CXZP(NX*NY*NZ,3) or CXZP(3*NX*NY*NZ)]
C
C       NX,NY,NZ        Size of grid in x,y,z directions (INTEGER).

C       AKD             Frequency of oscillation of dipoles and electric
C                       field; also absolute value of wave vector of
C                       incident wave [c=1]  (REAL*4).

C       CXZC            (NX+1) by (NY+1) by (NZ+1) by 6 array of Green
C                       function coefficients used
C                       internally by ESELF and
C                       *************NOT TO BE OVERWRITTEN***********
C                       between calls, because these coefficients are
C                       recomputed only if W has changed since the last
C                       call to ESELF.

C       CXZW            Complex, scratch-space vector of length:
C                       2*NX*2*NY*2*NY*3   --- for  FOURX FFT
C                       2*NX*2*NY*2*NZ*4   --- for CXFFT3 solver
C                       See comment about FFT usage and CMETHD
C                       flag.
C                       Can be overwritten between calls to ESELF

C       OUTPUT:

C       CXZE(I,J,K,L)   Lth component of dipole-induced electric field
C                       at grid point (I,J,K);
C                       the DECLARED length of ZE in the calling
C                       program is CXZE(NX,NY,NZ,3)
C                       [or CXZE(NX*NY*NZ,3) or CXZE(3*NX*NY*NZ)]
C
C Written by Jeremy J. Goodman, Princeton Univ. Observatory, 90.09.22
C History:
C 90.11.29 (BTD): Modified to set untransformed ZC(1,1,1,1-6)=0.
C 90.11.29 (PJF): Modified to use FOURX and CXFFT99
C 90.12.05 (BTD): Modified for new ordering of elements of polarization
C                 and electric field vectors in calling program.  Modified
C                 ESELF and PAD to remove distinction between NX,NY,NZ
C                 and dimensions CXZE and CXZP, since our new ordering
C                 always assumes this.
C 90.12.13 (BTD): Modified to include CONVEX option.
C 92.04.20 (BTD): removed ISYM from argument list of TRIM (was not used)
C 94.06.20 (PJF): modified to call CXFFT3N when CMETHD='NEWTMP'
C 96.10.18 (BTD): changed NEWTMP to GPFAFT
C 99.04.26 (BTD): changed notation CXY -> CXZE for clarity
C
C Copyright (C) 1993,1994,1996,1999 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      DATA WOLD/-999e-33/,CXUNIT/ (0.,1.)/,CXZERO/(0.,0.)/
C
      IF(AKD.EQ.WOLD)GOTO 70
C
C               AKD.NE.WOLD :
C
C               We have to recompute the
C               Green-function coefficients giving
C               components of field strength at a each grid point R
C               produced by unit-valued component of dipole moment at
C               point R', and then Fourier transform these components.
      WOLD=AKD
      NN(1)=2*NX
      NN(2)=2*NY
      NN(3)=2*NZ
      NGRID=8*NX*NY*NZ
C
C               Compute the actual coefficients:
C
      DO 60 K=1,NZ
         DO 50 J=1,NY
            DO 40 I=1,NX
               IF(I.EQ.1.AND.J.EQ.1.AND.K.EQ.1)GOTO 40
C Note: to use lattice with other than unit spacing in each direction,
C simply multiply RHS of each of next three lines by dx(1),dx(2),dx(3),
C declare real dx(3), and add dx to argument list
C This flexibility has been removed from present code for simplicity.
C 90.11.04 (BTD)
               X(1)=I-1
               X(2)=J-1
               X(3)=K-1
               R2=X(1)**2+X(2)**2+X(3)**2
               R=SQRT(R2)
               R3=R*R2
               CXPHAS=EXP(CXUNIT*AKD*R)
               M=0
               DO 30 II=1,3
                  M=M+1
                  CXFAC=(1.-CXUNIT*AKD*R)/R2
                  CXZC(I,J,K,M)=-CXPHAS/R3*
     &                          (AKD**2*(X(II)**2-R2)+
     &                          CXFAC*(R2-3*X(II)**2))
                  DO 20 JJ=II+1,3
                     M=M+1
                     XX=X(II)*X(JJ)
                     CXZC(I,J,K,M)=-CXPHAS/R3*
     &                             (AKD**2*XX-CXFAC*(3*XX))
   20             CONTINUE
   30          CONTINUE
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
C
C Set zero-lag elements to zero (they were skipped above)
      DO 65 M=1,6
         CXZC(1,1,1,M)=CXZERO
   65 CONTINUE
C
C          At this point, zc(i,j,k,*) contains the upper triangular
C          part of the symmetric 3 x 3 matrix giving the electric
C          field at grid point (i,j,k) caused by a dipole at (1,1,1)

C          Fill out zc to twice the size in each grid dimension
C          to accomodate negative lags [periodicity in each dimension is
C          assumed, so (e.g.)
C          nx < i <= 2*nx is equivalent to -nx < i <= 0],
C          exploiting symmetries, and Fourier transform.
C          Store only the first octant of the transform,
C          since the others can be obtained by symmetry.
C
      ISYM(1)=1
      ISYM(2)=1
      ISYM(3)=1
      CALL EXTND(CXZC(1,1,1,1),NX,NY,NZ,ISYM,CXZW(1,1,1,1))
      IF(CMETHD.EQ.'BRENNR')THEN
         CALL FOURX(CXZW(1,1,1,1),NN,3,+1)
      ELSEIF(CMETHD.EQ.'TMPRTN')THEN
         CALL CXFFT3(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
      ELSEIF(CMETHD.EQ.'GPFAFT')THEN
         CALL CXFFT3N(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1)
      ELSEIF(CMETHD.EQ.'CONVEX')THEN
         CALL C3DFFT(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
      ENDIF
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,ISYM,CXZC(1,1,1,1))
      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,CXZC(1,1,1,1))
C
      ISYM(1)=-1
      ISYM(2)=-1
      ISYM(3)=1
      CALL EXTND(CXZC(1,1,1,2),NX,NY,NZ,ISYM,CXZW(1,1,1,1))
      IF(CMETHD.EQ.'BRENNR')THEN
         CALL FOURX(CXZW(1,1,1,1),NN,3,+1)
      ELSEIF(CMETHD.EQ.'TMPRTN')THEN
         CALL CXFFT3(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
      ELSEIF(CMETHD.EQ.'GPFAFT')THEN
         CALL CXFFT3N(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1)
      ELSEIF(CMETHD.EQ.'CONVEX')THEN
         CALL C3DFFT(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
      ENDIF
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,ISYM,CXZC(1,1,1,2))
      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,CXZC(1,1,1,2))
C
      ISYM(1)=-1
      ISYM(2)=1
      ISYM(3)=-1
      CALL EXTND(CXZC(1,1,1,3),NX,NY,NZ,ISYM,CXZW(1,1,1,1))
      IF(CMETHD.EQ.'BRENNR')THEN
         CALL FOURX(CXZW(1,1,1,1),NN,3,+1)
      ELSEIF(CMETHD.EQ.'TMPRTN')THEN
         CALL CXFFT3(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
      ELSEIF(CMETHD.EQ.'GPFAFT')THEN
         CALL CXFFT3N(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1)
      ELSEIF(CMETHD.EQ.'CONVEX')THEN
         CALL C3DFFT(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
      ENDIF
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,ISYM,CXZC(1,1,1,3))
      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,CXZC(1,1,1,3))
C
      ISYM(1)=1
      ISYM(2)=1
      ISYM(3)=1
      CALL EXTND(CXZC(1,1,1,4),NX,NY,NZ,ISYM,CXZW(1,1,1,1))
      IF(CMETHD.EQ.'BRENNR')THEN
         CALL FOURX(CXZW(1,1,1,1),NN,3,+1)
      ELSEIF(CMETHD.EQ.'TMPRTN')THEN
         CALL CXFFT3(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
      ELSEIF(CMETHD.EQ.'GPFAFT')THEN
         CALL CXFFT3N(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1)
      ELSEIF(CMETHD.EQ.'CONVEX')THEN
         CALL C3DFFT(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
      ENDIF
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,ISYM,CXZC(1,1,1,4))
      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,CXZC(1,1,1,4))
C
      ISYM(1)=1
      ISYM(2)=-1
      ISYM(3)=-1
      CALL EXTND(CXZC(1,1,1,5),NX,NY,NZ,ISYM,CXZW(1,1,1,1))
      IF(CMETHD.EQ.'BRENNR')THEN
         CALL FOURX(CXZW(1,1,1,1),NN,3,+1)
      ELSEIF(CMETHD.EQ.'TMPRTN')THEN
         CALL CXFFT3(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
      ELSEIF(CMETHD.EQ.'GPFAFT')THEN
         CALL CXFFT3N(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1)
      ELSEIF(CMETHD.EQ.'CONVEX')THEN
         CALL C3DFFT(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
      ENDIF
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,ISYM,CXZC(1,1,1,5))
      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,CXZC(1,1,1,5))
C
      ISYM(1)=1
      ISYM(2)=1
      ISYM(3)=1
      CALL EXTND(CXZC(1,1,1,6),NX,NY,NZ,ISYM,CXZW(1,1,1,1))
      IF(CMETHD.EQ.'BRENNR')THEN
         CALL FOURX(CXZW(1,1,1,1),NN,3,+1)
      ELSEIF(CMETHD.EQ.'TMPRTN')THEN
         CALL CXFFT3(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
      ELSEIF(CMETHD.EQ.'GPFAFT')THEN
         CALL CXFFT3N(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,+1)
      ELSEIF(CMETHD.EQ.'CONVEX')THEN
         CALL C3DFFT(CXZW(1,1,1,1),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
      ENDIF
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,ISYM,CXZC(1,1,1,6))
      CALL TRIM(CXZW(1,1,1,1),NX,NY,NZ,CXZC(1,1,1,6))
C
C          End of recomputation of Green-function coefficients
C
   70 CONTINUE
C
C       Fourier transform the polarizations:
C
      DO 80 M=1,3
         CALL PAD(CXZP(1,1,1,M),NX,NY,NZ,CXZW(1,1,1,M))
         IF(CMETHD.EQ.'BRENNR')THEN
            CALL FOURX(CXZW(1,1,1,M),NN,3,+1)
         ELSEIF(CMETHD.EQ.'TMPRTN')THEN
            CALL CXFFT3(CXZW(1,1,1,M),2*NX,2*NY,2*NZ,+1,CXZW(1,1,1,4))
        ELSEIF(CMETHD.EQ.'GPFAFT')THEN
           CALL CXFFT3N(CXZW(1,1,1,M),2*NX,2*NY,2*NZ,+1)
         ELSEIF(CMETHD.EQ.'CONVEX')THEN
            CALL C3DFFT(CXZW(1,1,1,M),2*NX,2*NY,2*NZ,2*NX,2*NY,+1,IER)
         ENDIF
C
C***********************************************************************
   80 CONTINUE
C
C       Multiply by F.T. of Green function:
C
      DO 110 K=1,2*NZ
         KSGN=NINT(SIGN(1.,NZ+1.5-K))
         KR=MIN(K,2*NZ+2-K)
         DO 100 J=1,2*NY
            JSGN=NINT(SIGN(1.,NY+1.5-J))
            JR=MIN(J,2*NY+2-J)
            DO 90 I=1,2*NX
               ISGN=NINT(SIGN(1.,NX+1.5-I))
               IR=MIN(I,2*NX+2-I)
               CXXX=CXZC(IR,JR,KR,1)
               CXXY=CXZC(IR,JR,KR,2)*(ISGN*JSGN)
               CXXZ=CXZC(IR,JR,KR,3)*(ISGN*KSGN)
               CXYY=CXZC(IR,JR,KR,4)
               CXYZ=CXZC(IR,JR,KR,5)*(JSGN*KSGN)
               CXZZ=CXZC(IR,JR,KR,6)
               CXEX=CXXX*CXZW(I,J,K,1)+CXXY*CXZW(I,J,K,2)+
     &              CXXZ*CXZW(I,J,K,3)
               CXEY=CXXY*CXZW(I,J,K,1)+CXYY*CXZW(I,J,K,2)+
     &              CXYZ*CXZW(I,J,K,3)
               CXEZ=CXXZ*CXZW(I,J,K,1)+CXYZ*CXZW(I,J,K,2)+
     &              CXZZ*CXZW(I,J,K,3)
               CXZW(I,J,K,1)=CXEX
               CXZW(I,J,K,2)=CXEY
               CXZW(I,J,K,3)=CXEZ
   90       CONTINUE
  100    CONTINUE
  110 CONTINUE
C
C       Inverse Fourier transform to obtain electric field:
C
      DO 150 M=1,3
         IF(CMETHD.EQ.'BRENNR')THEN
            CALL FOURX(CXZW(1,1,1,M),NN,3,-1)
         ELSEIF(CMETHD.EQ.'TMPRTN')THEN
            CALL CXFFT3(CXZW(1,1,1,M),2*NX,2*NY,2*NZ,-1,CXZW(1,1,1,4))
         ELSEIF(CMETHD.EQ.'GPFAFT')THEN
           CALL CXFFT3N(CXZW(1,1,1,M),2*NX,2*NY,2*NZ,-1)
         ELSEIF(CMETHD.EQ.'CONVEX')THEN
            CALL C3DFFT(CXZW(1,1,1,M),2*NX,2*NY,2*NZ,2*NX,2*NY,-1,IER)
         ENDIF
C
C***********************************************************************
C Note: the Convex FFT routine already normalizes result.
C       For other FFT routines need to divide result by NGRID
C
         IF(CMETHD.EQ.'CONVEX')THEN
            DO 118 K=1,NZ
               DO 116 J=1,NY
                  DO 114 I=1,NX
                     CXZE(I,J,K,M)=CXZW(I,J,K,M)
  114             CONTINUE
  116          CONTINUE
  118       CONTINUE
         ELSE
            FAC=1./NGRID
            DO 140 K=1,NZ
               DO 130 J=1,NY
                  DO 120 I=1,NX
                     CXZE(I,J,K,M)=FAC*CXZW(I,J,K,M)
  120             CONTINUE
  130          CONTINUE
  140       CONTINUE
         ENDIF
  150 CONTINUE
C
      RETURN
      END
C***********************************************************************
C
      SUBROUTINE EXTND(CXA,NX,NY,NZ,ISYM,CXB)
C
C       Using symmetries, extend first octant of coefficient matrix to
C       other octants.
C
C Originally written by J.J. Goodman
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      INTEGER NX,NY,NZ
      COMPLEX CXA(NX+1,NY+1,NZ+1),CXB(2*NX,2*NY,2*NZ)
      INTEGER ISYM(3)
C Local variables:
      INTEGER I,J,K
      COMPLEX CXZERO
C SAVE statement:
      SAVE CXZERO
C DATA statement:
      DATA CXZERO/(0.,0.)/
C***********************************************************************
      DO 30 K=1,NZ
         DO 20 J=1,NY
            DO 10 I=1,NX
               CXB(I,J,K)=CXA(I,J,K)
   10       CONTINUE
   20    CONTINUE
   30 CONTINUE
C
      DO 60 K=1,NZ
         DO 50 J=1,NY
            CXB(NX+1,J,K)=CXZERO
            DO 40 I=NX+2,2*NX
               CXB(I,J,K)=CXA(2*NX+2-I,J,K)*ISYM(1)
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
C
      DO 90 K=1,NZ
         DO 80 I=1,2*NX
            CXB(I,NY+1,K)=CXZERO
            DO 70 J=NY+2,2*NY
               CXB(I,J,K)=CXB(I,2*NY+2-J,K)*ISYM(2)
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
C
      DO 120 J=1,2*NY
         DO 110 I=1,2*NX
            CXB(I,J,NZ+1)=CXZERO
            DO 100 K=NZ+2,2*NZ
               CXB(I,J,K)=CXB(I,J,2*NZ+2-K)*ISYM(3)
  100       CONTINUE
  110    CONTINUE
  120 CONTINUE
C
      RETURN
      END
C***********************************************************************
C
      SUBROUTINE TRIM(CXB,NX,NY,NZ,CXA)
C
C Copy the first octant of CXB into CXA
C
C Originally written by J.J. Goodman
C History:
C 92.04.20 (BTD) removed ISYM from argument list of TRIM:
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      INTEGER NX,NY,NZ
      COMPLEX CXA(NX+1,NY+1,NZ+1),CXB(2*NX,2*NY,2*NZ)
C Local variables:
      INTEGER I,J,K
C***********************************************************************
      DO 30 K=1,NZ+1
         DO 20 J=1,NY+1
            DO 10 I=1,NX+1
               CXA(I,J,K)=CXB(I,J,K)
   10       CONTINUE
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
C***********************************************************************
C
      SUBROUTINE PAD(CXA,NX,NY,NZ,CXB)
C
C Pad the array CXA with zeros out to twice its length in each
C dimension, and put the result in CXB.
C
C Originally written by J.J. Goodman
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      INTEGER NX,NY,NZ
      COMPLEX CXA(NX,NY,NZ),CXB(2*NX,2*NY,2*NZ)
C Local variables:
      COMPLEX CXZERO
      INTEGER I,J,K
      SAVE CXZERO
      DATA CXZERO/(0.,0.)/
C***********************************************************************
      DO 30 K=1,2*NZ
         DO 20 J=1,2*NY
            DO 10 I=1,2*NX
               CXB(I,J,K)=CXZERO
   10       CONTINUE
   20    CONTINUE
   30 CONTINUE
C
      DO 60 K=1,NZ
         DO 50 J=1,NY
            DO 40 I=1,NX
               CXB(I,J,K)=CXA(I,J,K)
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
C
      RETURN
      END

      SUBROUTINE EVALA(CXADIA,CXALPH,MXN3,NAT)
C
C*********************************************************************
C Purpose: to evaluate diagonal elements of matrix A for use in
C discrete-dipole scattering calculation.
C
C Given:   CXALPH(1-NAT3)=complex polarizabilities of dipoles
C          NAT=number of dipoles
C          MXN3,MXNAT=dimensioning information
C Returns: CXADIAG(1-3*NAT)=diagonal elements of matrix A
C
C B. T. Draine, Princeton Univ. Obs., 87/1/7
C History:
C 88/ 5/ 5 (BTD): Modifications...
C 90/11/ 6 (BTD): Modified to pass array dimensions
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C*** Arguments:
      INTEGER MXN3,NAT
      COMPLEX CXADIA(MXN3),CXALPH(MXN3)
C
C*** Local variables:
      INTEGER IAX
C
C Compute "diagonal" elements of A: those having I=J
      DO 100 IAX=1,3*NAT
          CXADIA(IAX)=1./CXALPH(IAX)
  100 CONTINUE
      RETURN
C
      END
      SUBROUTINE EVALE(CXE00,AKD,IXYZ0,MXNAT,MXN3,NAT,NAT0,NX,NY,NZ,CXE)
C*** Arguments:
      INTEGER MXN3,MXNAT,NAT,NAT0,NX,NY,NZ
      INTEGER*2 IXYZ0(MXNAT,3)
      REAL AKD(3)
C Note: CXE should be dimensioned to CXE(NAT,3) in this routine
C       so that first 3*NAT elements of CXE are employed.
C       XYZ0 should be dimensioned to 
      COMPLEX CXE(NAT,3),CXE00(3)
C
C***  Local variables:
      COMPLEX CXFAC,CXI
      REAL GPI,X,X1,X2
      INTEGER IA,IX,IY,IZ,M
C
C*** Intrinsic functions:
      INTRINSIC ATAN,EXP,REAL
C
C*** Data statements:
      DATA CXI/(0.,1.)/
C***********************************************************************
C
C Given:   CXE00(1-3)=Incident E field at origin (complex) at t=0
C          AKD(1-3)=(kx,ky,kz)*d for incident wave (d=lattice spacing)
C          IXYZ0(1-NAT0,3)=x/d,y/d,z/d for each of NAT0 physical dipoles
C          MXNAT,MXN3=dimensioning information
C          NAT0=number of dipoles in physical target
C          NAT=number of locations at which to calculate CXE
C
C Returns: CXE(1-NAT,3)=incident E field at NAT locations at t=0
C
C B.T.Draine, Princeton Univ. Obs., 88.05.09
C History:
C 90.11.06 (BTD): Modified to pass array dimension.
C 90.11.29 (BTD): Modified to allow option for either
C                   physical locations only (NAT=NAT0), or
C                   extended dipole array (NAT>NAT0)
C 90.11.30 (BTD): Corrected error for case NAT>NAT0
C 90.11.30 (BTD): Corrected another error for case NAT>NAT0
C 90.12.03 (BTD): Change ordering of XYZ0 and CXE
C 90.12.05 (BTD): Corrected error in dimensioning of CXE
C 90.12.10 (BTD): Remove XYZ0, replace with IXYZ0
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
      GPI=4.*ATAN(1.)
C
C Evaluate electric field vector at each dipole location.
C
C If NAT=NAT0, then evaluate E only at occupied sites.
C If NAT>NAT0, then evaluate E at all sites.
C
      IF(NAT.EQ.NAT0)THEN
         DO 10 IA=1,NAT0
            X=0.
            DO 5 M=1,3
               X=X+AKD(M)*IXYZ0(IA,M)
    5       CONTINUE
            CXFAC=EXP(CXI*X)
            DO 7 M=1,3
               CXE(IA,M)=CXE00(M)*CXFAC
    7       CONTINUE
   10    CONTINUE
      ELSE
         IA=0
         DO 20 IZ=1,NZ
            X1=AKD(3)*REAL(IZ)
            DO 19 IY=1,NY
               X2=X1+AKD(2)*REAL(IY)
               DO 18 IX=1,NX
                  IA=IA+1
                  X=X2+AKD(1)*REAL(IX)
                  CXFAC=EXP(CXI*X)
                  DO 17 M=1,3
                     CXE(IA,M)=CXE00(M)*CXFAC
   17             CONTINUE
   18          CONTINUE
   19       CONTINUE
   20    CONTINUE
      ENDIF
      RETURN
C
      END
      SUBROUTINE EVALQ(CXALPH,AK,NAT3,E02,CXE,CXP,CABS,CEXT,CPHA,MXN3,
     &                 IMETHD)
C*** Arguments:
      REAL CABS,CEXT,CPHA
      INTEGER IMETHD,MXN3,NAT3
      COMPLEX CXE(MXN3),
     &        CXP(MXN3),
     &        CXALPH(MXN3)
      REAL AK(3)
C
C*** Local variables:
      COMPLEX CXA,CXI
      REAL AK1,AK2,AK3,ALPHA2,ALPHSQ,E02,PI,RABS
      INTEGER J
C
C*** Intrinsic functions:
      INTRINSIC AIMAG,CONJG,REAL,SQRT
C
C*** SAVE statements:
      SAVE CXI
C
C*** Data statements:
C
      DATA CXI/(0.,1.)/
C***********************************************************************
C
C Given: CXALPH(1-NAT3) = polarizability vector
C        AK(1-3) = (k_x,k_y,k_z)*d (d=lattice spacing)
C        NAT3 = 3*number of dipoles
C        E02 = |E_0|^2 , where E_0 = incident complex E field.
C        CXE(1-NAT3) = components of E field at each dipole
C        CXP(1-NAT3) = components of polarization vector at each dipole
C        IMETHD = 0 or 1
C Finds: CEXT = extinction cross section /d**2
C  and, if IMETHD=1, also computes
C        CPHA = phase-lag cross section /d**2
C        CABS = absorption cross section /d**2
C        
C Note: present definition of CPHA is 1/2 of Martin's definition
C
C B.T.Draine, Princeton Univ. Obs., 87/1/4
C
C History:
C 88.04.28 (BTD): modifications
C 90.11.02 (BTD): modified to allow use of vacuum sites (now pass E02
C                 from calling routine instead of evaluating it here)
C 90.12.13 (BTD): modified to use IMETHD flag, to allow "fast" calls
C                 in which only CEXT is computed.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Initialization:
C Compute magnitude of kd
C
      CEXT=0.
      CABS=0.
      CPHA=0.
      PI=4.*ATAN(1.)
      AK2=0.
      DO 10 J=1,3
         AK2=AK2+AK(J)*AK(J)
   10 CONTINUE
      AK1=SQRT(AK2)
      AK3=AK1*AK2
C
C Initialization complete.
      IF(IMETHD.EQ.0)THEN
C*** Compute CEXT:
         DO 20 J=1,NAT3
            CEXT=CEXT+AIMAG(CXP(J))*REAL(CXE(J))
     &               -REAL(CXP(J))*AIMAG(CXE(J))
   20    CONTINUE
         CEXT=4.*PI*AK1*CEXT/E02
      ELSEIF(IMETHD.EQ.1)THEN
C Compute CABS
         DO 40 J=1,NAT3
            ALPHA2=AIMAG(CXALPH(J))
            ALPHSQ=CONJG(CXALPH(J))*CXALPH(J)
c*** diagnostic
c      if(alphsq.le.0.)then
c         write(0,*)'in evalq, having trouble, nat3=',nat3
c         write(0,*)'j=',j,' cxalpha(j)=',cxalph(j)
c      endif
c***
            RABS=ALPHA2/ALPHSQ-AK3/1.5
            CABS=CABS+CONJG(CXP(J))*CXP(J)*RABS
   40    CONTINUE
         CABS=4.*PI*AK1*CABS/E02
C Compute CEXT and CPHA
         CXA=(0.,0.)
         DO 50 J=1,NAT3
            CXA=CXA+CXP(J)*CONJG(CXE(J))
   50    CONTINUE
         CEXT=4.*PI*AK1*AIMAG(CXA)/E02
         CPHA=2.*PI*AK1*REAL(CXA)/E02
      ENDIF
      RETURN
C
      END
      SUBROUTINE EXTEND(CMETHD,ICOMP,ICOMP2,IDVERR,IDVOUT,IOCC,IX,IY,IZ,
     &                  MXNX,MXNY,MXNZ,MXNAT,MXN3,NAT,NAT0,NAT3,
     &                  NX,NY,NZ)
C Arguments:
      CHARACTER*6 CMETHD
      INTEGER
     &   IDVERR,IDVOUT,MXNAT,MXN3,MXNX,MXNY,MXNZ,
     &   NAT,NAT0,NAT3,NX,NY,NZ
      INTEGER*2
     &   ICOMP(MXN3),ICOMP2(MXN3),IOCC(MXNAT),
     &   IX(MXNAT),IY(MXNAT),IZ(MXNAT)
C Local parameters:
      INTEGER MX235
      PARAMETER(MX235=136)
C Local variables:
      INTEGER IXA,IXMAX,IXMIN,IYA,IYMAX,IYMIN,
     &        IZA,IZMAX,IZMIN,JA,JA2,NXY
      INTEGER NF235(MX235)
      SAVE NF235
C***********************************************************************
C Given:
C    NAT0 = number of dipoles in real target
C    IX(1-NAT0) = x-index of dipoles in desired target
C    IY(1-NAT0) = y-index of dipoles in desired target
C    IZ(1-NAT0) = z-index of dipoles in desired target
C    ICOMP(1-NAT0,1-3) = composition info. for dipoles in desired target
C                       [with storage locations according to
C                        dimensioning ICOMP(MXNAT,3) ]
C    MXNX,MXNY,MXNZ = maximum allowed target size in x,y,z directions
C    MXNAT,MXN3 = dimensioning information
C
C Returns:
C    NAT = number of dipoles in extended target = NX*NY*NZ
C    NAT3 = 3*NAT
C    IX(1-NAT0) = x-index of dipoles in reordered physical target
C    IY(1-NAT0) = y-index of dipoles in reordered physical target
C    IZ(1-NAT0) = z-index of dipoles in reordered physical target
C    ICOMP(1-NAT,1-3) = composition inf. for dipoles (ICOMP=0 if vacant)
C                       [with storage locations according to
C                        dimensioning ICOMP(NAT,3) ]
C    NX = range of IX values (1 to NX)
C    NY = range of IY values (1 to NY)
C    NZ = range of IZ values (1 to NZ)
C
C Note: it is necessary to reorder IX,IY,IZ because DDSCAT calls routine
C    REDUCE to take list of polarizations for "extended" target and
C    collapse it to list of polarizations for "physical" target.
C    Similarly, REDUCE is also used to collapse vector of incident E
C    and vector of polarizabilities.
C
C B.T.Draine, Princeton Univ. Obs., 90.10.30
C History:
C 90.11. 8 (BTD): Changed DO...ENDDO to DO #...# CONTINUE
C 90.11.30 (BTD): Cosmetic changes.
C 90.11.30 (BTD): Modify to choose NX,NY,NZ from NF235.
C 90.12.05 (BTD): Change argument list; correct bug.
C 90.12.07 (BTD): Previous version was stupidly written. Rewritten.
C 90.12.13 (BTD): Corrected data list for NF235. Modified to allow
C                 CMETHD='CONVEX' option.
C 92.06.01 (BTD+PJF): Corrected major flaw in EXTEND: previous version
C                 neglected to reorder elements of IX,IY,IZ to
C                 correspond to possible reordering of occupied sites
C                 when extended target is constructed.  This error was
C                 previously overlooked because ordering was fortuitously
C                 unchanged by extend for targets generated by TARREC,
C                 TARELL, TARCYL, TARHEX.  Note that differential
C                 scattering cross sections computed previously using
C                 DDSCAT ver4a and using TARTET are not correct.
C 92.09.09 (BTD): extended values of NF235 up to 1024, and added
C                 warning if calling with IXMAX,IYMAX,or IZMAX > 1024
C 93.01.20 (BTD): added MXNX,MXNY,MXNZ to argument list and
C                 added tests to ensure that target is consistent with
C                 limits MXNX,MXNY,MXNZ
C 93.03.12 (BTD): extended NF235 up to 4096
C 94.06.20 (PJF): added code for NEWTMP 3D FFT
C 96.10.18 (BTD): changed NEWTMP to GPFAFT (Generalized Prime Factor
C                 Algorithm for Fourier Transform).
C end history
C Copyright (C) 1993,1996 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Elements of NF235 are numbers.le.4096  which can be factored by 2,3,5:
      DATA NF235/2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,
     &     45,48,50,54,60,64,72,75,80,81,90,96,100,108,120,125,128,135,
     &     144,150,160,162,180,192,200,216,225,240,243,250,256,270,288,
     &     300,320,324,360,375,384,400,405,432,450,480,486,500,512,540,
     &     576,600,625,640,648,675,720,729,750,768,800,810,864,900,960,
     &     972,1000,1024,1080,1125,1152,1200,1215,1250,1280,1296,1350,
     &     1440,1458,1500,1536,1600,1620,1728,1800,1875,1920,1944,2000,
     &     2025,2048,2160,2187,2250,2304,2400,2430,2500,2560,2592,2700,
     &     2880,2916,3000,3072,3125,3200,3240,3375,3456,3600,3645,3750,
     &     3840,3888,4000,4050,4096/
C*** First determine max., min values of input IX,IY,IZ:
      IXMIN=IX(1)
      IXMAX=IX(1)
      IYMIN=IY(1)
      IYMAX=IY(1)
      IZMIN=IZ(1)
      IZMAX=IZ(1)
      DO 100 JA=2,NAT0
         IF(IX(JA).LT.IXMIN)IXMIN=IX(JA)
         IF(IX(JA).GT.IXMAX)IXMAX=IX(JA)
         IF(IY(JA).LT.IYMIN)IYMIN=IY(JA)
         IF(IY(JA).GT.IYMAX)IYMAX=IY(JA)
         IF(IZ(JA).LT.IZMIN)IZMIN=IZ(JA)
         IF(IZ(JA).GT.IZMAX)IZMAX=IZ(JA)
  100 CONTINUE
C
C*** Now shift so that IX runs from 1 to IXMAX-IXMIN+1, etc.
      DO 110 JA=1,NAT0
         IX(JA)=IX(JA)+1-IXMIN
         IY(JA)=IY(JA)+1-IYMIN
         IZ(JA)=IZ(JA)+1-IZMIN
  110 CONTINUE
      IXMAX=IXMAX-IXMIN+1
      IYMAX=IYMAX-IYMIN+1
      IZMAX=IZMAX-IZMIN+1
C
      IF(IXMAX.LE.MXNX.AND.IYMAX.LE.MXNY.AND.IZMAX.LE.MXNZ)THEN
         WRITE(IDVOUT,6100)IXMAX,IYMAX,IZMAX
      ELSE
         WRITE(IDVOUT,6200)IXMAX,IYMAX,IZMAX,MXNX,MXNY,MXNZ
         CALL ERRMSG('FATAL','EXTEND',
     &   ' need to change MXNX,MXNY,MXNZ in DDSCAT and recompile')
      ENDIF
C
C*** Now determine NX,NY,NZ
      IF(CMETHD.EQ.'BRENNR')THEN
C Using FOURX: no restriction on NX,NY,NZ:
         NX=IXMAX
         NY=IYMAX
         NZ=IZMAX
      ELSEIF(CMETHD.EQ.'TMPRTN'.OR.CMETHD.EQ.'CONVEX'
     $   .OR.CMETHD.EQ.'GPFAFT')THEN
C Using Temperton's routine CFFT99F or Convex veclib routine C3DFFT or
C Temperton's new GPFA routine:
C     NX,NY,NZ must be factorable by 2,3,5  !!
C We require that NX,NY,NZ be selected from list NF235
C Check that list is long enough:
         IF(IXMAX.GT.NF235(MX235))THEN
            WRITE(IDVERR,*)'EXTEND> Error: IXMAX=',IXMAX,
     &                     '.gt.NF235(MX235)=',NF235(MX235)
            CALL ERRMSG('FATAL','EXTEND',
     &                  ' need to add to NF235 in EXTEND')
         ELSEIF(IYMAX.GT.NF235(MX235))THEN
            WRITE(IDVERR,*)'EXTEND> Error: IYMAX=',IYMAX,
     &                     '.gt.NF235(MX235)=',NF235(MX235)
            CALL ERRMSG('FATAL','EXTEND',
     &                  ' need to add to NF235 in EXTEND')
         ELSEIF(IZMAX.GT.NF235(MX235))THEN
            WRITE(IDVERR,*)'EXTEND> Error: IZMAX=',IZMAX,
     &                     '.gt.NF235(MX235)=',NF235(MX235)
            CALL ERRMSG('FATAL','EXTEND',
     &                  ' need to add to NF235 in EXTEND')
         ENDIF
         NX=0
         NY=0
         NZ=0
         DO 115 JA=1,MX235
            IF(IXMAX.LE.NF235(JA).AND.NX.EQ.0)NX=NF235(JA)
            IF(IYMAX.LE.NF235(JA).AND.NY.EQ.0)NY=NF235(JA)
            IF(IZMAX.LE.NF235(JA).AND.NZ.EQ.0)NZ=NF235(JA)
  115    CONTINUE
      ENDIF
      IF(NX.GT.MXNX.OR.NY.GT.MXNY.OR.NZ.GT.MXNZ)THEN
         WRITE(IDVOUT,6300)NX,NY,NZ,MXNX,MXNY,MXNZ
         CALL ERRMSG('FATAL','EXTEND',
     &   ' need to change MXNX,MXNY,MXNZ in DDSCAT and recompile')
      ENDIF
      NAT=NX*NY*NZ
      NAT3=3*NAT
C
C*** Now extend target: generate vectors of dipole properties
C        ICOMP(jx,jy,jz,1-3)=ICOMP(ja,1-3)=ICOMP(ja3)
C        IOCC(jx,jy,jz)=IOCC(ja)
C    with indices jx,jy,jz running from 1-NX,1-NY,1-NZ
C              or ja          "     "   1-(NX*NY*NZ)
C              or ja3         "     "   1-(3*NX*NY*NZ)
C
C
C First initialize ICOMP and IOCC to zero:
      DO 200 JA2=1,NAT
         IOCC(JA2)=0
         ICOMP2(JA2)=0
         ICOMP2(JA2+NAT)=0
         ICOMP2(JA2+2*NAT)=0
  200 CONTINUE
C Now reset values at occupied sites
      NXY=NX*NY
      DO 300 JA=1,NAT0
         JA2=IX(JA)+NX*(IY(JA)-1)+NXY*(IZ(JA)-1)
         IOCC(JA2)=1
         ICOMP2(JA2)=ICOMP(JA)
         ICOMP2(JA2+NAT)=ICOMP(JA+MXNAT)
         ICOMP2(JA2+2*NAT)=ICOMP(JA+2*MXNAT)
  300 CONTINUE
C
C Now overwrite original ICOMP
      DO 400 JA=1,NAT3
         ICOMP(JA)=ICOMP2(JA)
  400 CONTINUE
C
C Now reorder IX,IY,IZ to correspond to new ordering of occupied sites
      JA=0
      DO 500 JA2=1,NAT
         IF(IOCC(JA2).EQ.1)THEN
            JA=JA+1
            IXA=1+MOD(JA2-1,NX)
            IYA=1+MOD(JA2-IXA,NXY)/NX
            IZA=1+(JA2-IXA-NX*(IYA-1))/NXY
            IX(JA)=IXA
            IY(JA)=IYA
            IZ(JA)=IZA
         ENDIF
  500 CONTINUE
      RETURN
 6100 FORMAT(' >EXTEND: target extent in x,y,z directions =',I3,',',I3,
     &       ',',I3,' (in Target Frame)')
 6200 FORMAT(' >EXTEND: Fatal error!:',/,
     &       '          target extent in x,y,z directions =',I3,',',I3,
     &       ',',I3,/,
     &       '          exceeds dimension limits MXNX,MXNY,MXNZ=',I3,
     &       ',',I3,',',I3)
 6300 FORMAT(' >EXTEND: Fatal error after extending target!:',/,
     &       '          Extended target extent in x,y,z directions =',
     &       I3,',',I3,',',I3,/,
     &       '          exceeds dimension limits MXNX,MXNY,MXNZ=',I3,
     &       ',',I3,',',I3)
      END