      SUBROUTINE SCAT(AK,AKS,EM1,EM2,E02,
     &    CBKSCA,CSCA,CSCAG,CXE00,CXF1L,CXF2L,CXP,CXSCR1,CXSCR2,CXSCR3,
     &    ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR,
     &    PHIN,SCRRS1,THETAN,IXYZ)
C
C Arguments:
      REAL CBKSCA,CSCA,CSCAG,E02
      INTEGER ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR
      INTEGER*2 IXYZ(MXNAT,3)
      COMPLEX CXE00(3),CXF1L(MXSCA),CXF2L(MXSCA),CXP(NAT,3),
     &        CXSCR1(MXN3),CXSCR2(MXN3),CXSCR3(MXN3)
      REAL AK(3),AKS(3,MXSCA),EM1(3,MXSCA),EM2(3,MXSCA),
     $     PHIN(MXSCA),SCRRS1(MXN3),THETAN(MXSCA)
C
C Local variables:
      COMPLEX CXI
      REAL AK2,AK3,AKK,CONST,COSPHI,COSTH,DCOSTH,DPHI,EINC,ES2,
     &     PI,PHI,RRR,SINPHI,SINTH,TERM,W,WTH
      INTEGER ICOSTH,IPHI,IPHM,J,K,ND
      COMPLEX CXES(3)
      REAL AKS0(3),AKS1(3),AKS2(3)
C
C Intrinsic functions:
      INTRINSIC COS,EXP,MIN,SIN,SQRT,CONJG,REAL
C
C Data statements:
      DATA CXI/(0.,1.)/
C***********************************************************************
C
C SCAT computes energy radiated by array of oscillating dipoles
C and corresponding scattering properties if oscillations are in
C response to incident E field.
C
C Given:
C     AK(1-3) = KX,KY,KZ = components of incident k-vector
C     AKS(3,MXSCA)=scattered k vectors
C     EM1(3,MXSCA)=pol.vector 1 for each scattering direction
C     EM2(3,MXSCA)=pol.vector 2 for each scattering direction
C     E02 = |E_0|^2 , where E_0 = incident complex E field
C     NAT = number of dipoles
C     NAT3 = 3*NAT
C     XYZ(NAT,1-3) = X,Y,Z for dipoles 1-NAT
C     CXE00(1-3) = (complex) reference polarization vector
C     CXP(1-NAT,1-3) = components of polarization of each dipole
C     ICTHM = number of THETA values for calculation of CSCA and G
C             (ICTHM must be odd; recommend ICTHM > 5*(2*pi*a/lambda)
C              for accuracy)
C     IPHIM = number of PHI values for calculation of CSCA and G
C     NDIR = number of directions at which scattering matrix elements
C            are to be computed
C     THETAN(1-NDIR) : scattering angle (radians);
C                      cos(THETAN) = n0 dot n
C     PHIN(1-NDIR)   : scattering angle (radians);
C                      PHIN is azimuthal angle of n from plane defined
C                      by n0 and Re(CXE00)
C
C Returns:
C     CBKSCA=differential scattering cross section for theta=pi
C            (lattice units)
C     CSCA = scattering cross section (lattice units)
C     G = <cos(theta)> , where theta=scattering angle
C     F1L(NDIR) = f_{1L} for direction NDIR
C     F2L(NDIR) = f_{2L} for direction NDIR
C                where f_{1L} connects incident polarization L to
C                outgoing theta polarization (E in scattering plane)
C                f_{2L} connects incident polarization L to outgoing phi
C                polarization (E perpendicular to scattering plane)
C                This is same f_{ml} notation as used by Draine (1988;
C                Ap.J., 333, 848).
C
C Notes:
C Scratch vectors introduced to permit vectorization:
C     SCRRS1 = real vector of length.GE.NAT
C     SCRCS1,SCRCS2,SCRCS3 = complex vectors of length.GE.NAT
C
C B.T.Draine, Princeton Univ. Obs., 87/1/4
C History:
C 88.06.29 (BTD): modified
C 89.11.28 (BTD): modified
C 90.11.02 (BTD): modified to enable use with vacuum sites.
C 90.11.20 (BTD): corrected error in calculation of f_ml (missing
C                 factor of k).
C 90.12.03 (BTD): changed ordering of XYZ0 and CXP
C                 eliminated cxe0 from argument list
C 90.12.10 (BTD): replaced XYZ0 with IXYZ0 to reduce memory use
C 91.08.14 (BTD): add CBKSCA to argument list for SCAT
C                 add code to eliminate sum over phi when theta=0 or pi
C                 add code to calculate CBKSCA
C 94.04.04 (BTD): removed superfluous CMPLX( ) operation from 2 lines
C                 in DO 200 loop. Also added new variable TERM to move
C                 some computation out of last DO loop.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      PI=4.*ATAN(1.)
C
C ICTHM,IPHIM determine angular grid used for calculation of total
C scattering cross section and G=<cos(theta)>
C ICTHM=number of theta values for which scattering is evaluated
C       (0 to pi)
C IPHIM=number of phi values for which scattering is evaluated
C       (0 to 2*pi)
C Note: ICTHM must be odd as Simpson's rule integration is employed.
C Check whether ICTHM is odd, and increase by 1 if not.
C
      IF(ICTHM.EQ.2*(ICTHM/2))ICTHM=ICTHM+1
C
C Incident k vector AK defines one axis (for spherical integration)
C Use (real part of) reference polarization vector CXE00 to define
C second coordinate axis.
C Construct third coordinate axis by taking cross product of first
C two axes.
C AKS1, AKS2 = vectors of length equal to AK lying along second and
C              third coordinate axes.
C
      AK2=AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3)
      AKK=SQRT(AK2)
      AK3=AKK*AK2
      AKS1(1)=REAL(CXE00(1))
      AKS1(2)=REAL(CXE00(2))
      AKS1(3)=REAL(CXE00(3))
C EINC = magnitude of (real part of) reference electric field
      EINC=SQRT(AKS1(1)**2+AKS1(2)**2+AKS1(3)**2)
      IF(EINC.LE.0.)CALL ERRMSG('FATAL','SCAT',
     &                          'CXE00 IS PURE IMAGINARY!')
      RRR=AKK/EINC
      AKS1(1)=RRR*AKS1(1)
      AKS1(2)=RRR*AKS1(2)
      AKS1(3)=RRR*AKS1(3)
      AKS2(1)=(AK(2)*AKS1(3)-AK(3)*AKS1(2))/AKK
      AKS2(2)=(AK(3)*AKS1(1)-AK(1)*AKS1(3))/AKK
      AKS2(3)=(AK(1)*AKS1(2)-AK(2)*AKS1(1))/AKK
C
      DCOSTH=2./REAL(ICTHM-1)
      DPHI=2.*PI/REAL(IPHIM)
      CSCA=0.
      CSCAG=0.
C
C CONST = multiplicative factor entering weights
C factor 2/(ICTHM-1) = d(cos(theta))
C factor 2*pi/IPHIM = d(phi)
C factor 3 in denom from Simpson's rule (later mult. by 1,4,2..)
C Use Simpson's rule integration over d(cos(theta))
      CONST=4.*PI/ (3.*REAL((ICTHM-1)*IPHIM))
C Note that for cos(theta)=-1 or 1 we do not need to sum over phi,
C so treat these two directions separately
      DO 80 ICOSTH=2,ICTHM-1
         WTH=4.
         IF(ICOSTH.GT.2*(ICOSTH/2)) WTH=2.
C         IF(ICOSTH.EQ.1) WTH=1.
C         IF(ICOSTH.EQ.ICTHM) WTH=1.
         COSTH=1.-REAL(ICOSTH-1)*DCOSTH
         SINTH=SQRT(1.-MIN(COSTH*COSTH,1.))
         IPHM=IPHIM
C Now evaluate total weight factor W
         W=CONST*WTH
         DO 70 IPHI=1,IPHM
            PHI=DPHI*(REAL(IPHI)-0.5)
            SINPHI=SIN(PHI)
            COSPHI=COS(PHI)
C Now compute scattered k vector
C Initialize CXES(1-3) = scattered electric field
            DO 10 J=1,3
               AKS0(J)=COSTH*AK(J)+SINTH*
     &                 (COSPHI*AKS1(J)+SINPHI*AKS2(J))
               CXES(J)=(0.,0.)
   10       CONTINUE
C
C Now compute CXES(1-3) = (r/k*k)*exp(-ikr+iwt)*E(ks,r,t)
C      where E(ks,r,t) is complex electric field vector propagating
C      with wave vector ks, at distance r from origin
C
            DO 50 J=1,NAT
               CXSCR1(J)=(0.,0.)
C Compute SCRCS1(J)=ks dot p(j) / ks**2
               DO 20 K=1,3
                  CXSCR1(J)=CXSCR1(J)+AKS0(K)*CXP(J,K)
   20          CONTINUE
               CXSCR1(J)=CXSCR1(J)/AK2
               SCRRS1(J)=0.
C SCRRS1(J) is ks dot r(j)
               DO 30 K=1,3
                  SCRRS1(J)=SCRRS1(J)+IXYZ(J,K)*AKS0(K)
   30          CONTINUE
C CXSCR2(J) is exp(-i*ks*r(j)) factor for atom j
               CXSCR2(J)=EXP(-CXI*SCRRS1(J))
               DO 40 K=1,3
                  CXES(K)=CXES(K)+(CXP(J,K)-AKS0(K)*CXSCR1(J))*
     &                    CXSCR2(J)
   40          CONTINUE
   50       CONTINUE
C
C Have completed computation of vector ES(1-3) for scattering direction
C (THETA,PHI).
C
            ES2=0.0
            DO 60 K=1, 3
               ES2=ES2+CXES(K)*CONJG(CXES(K))
   60       CONTINUE
            CSCA=CSCA+W*ES2
            CSCAG=CSCAG+W*COSTH*ES2
   70    CONTINUE
   80 CONTINUE
C***
C Now do special case theta=0 and theta=pi
      W=CONST*IPHM
      DO 180 ICOSTH=-1,1,2
C Compute scattered k vector
C Initialize CXES(1-3) = scattered electric field
         DO 110 J=1,3
            AKS0(J)=ICOSTH*AK(J)
            CXES(J)=(0.,0.)
  110    CONTINUE
C
C Now compute CXES(1-3) = (r/k*k)*exp(-ikr+iwt)*E(ks,r,t)
C      where E(ks,r,t) is complex electric field vector propagating
C      with wave vector ks, at distance r from origin
C
         DO 150 J=1,NAT
            CXSCR1(J)=(0.,0.)
C Compute SCRCS1(J)=ks dot p(j) / ks**2
            DO 120 K=1,3
               CXSCR1(J)=CXSCR1(J)+AKS0(K)*CXP(J,K)
  120       CONTINUE
            CXSCR1(J)=CXSCR1(J)/AK2
            SCRRS1(J)=0.
C SCRRS1(J) is ks dot r(j)
            DO 130 K=1,3
               SCRRS1(J)=SCRRS1(J)+IXYZ(J,K)*AKS0(K)
  130       CONTINUE
C CXSCR2(J) is exp(-i*ks*r(j)) factor for atom j
            CXSCR2(J)=EXP(-CXI*SCRRS1(J))
            DO 140 K=1,3
               CXES(K)=CXES(K)+(CXP(J,K)-AKS0(K)*CXSCR1(J))*CXSCR2(J)
  140       CONTINUE
  150    CONTINUE
C
C Have completed computation of vector ES(1-3) for scattering direction
C theta=-1 or 1
C
         ES2=0.0
         DO 160 K=1, 3
            ES2=ES2+CXES(K)*CONJG(CXES(K))
  160    CONTINUE
         CSCA=CSCA+W*ES2
         CSCAG=CSCAG+W*ICOSTH*ES2
C Save backscattering cross section
         IF(ICOSTH.EQ.-1)CBKSCA=AK2*AK2*ES2/E02
  180 CONTINUE
C***
      CSCA=AK2*AK2*CSCA/E02
      CSCAG=AK2*AK2*CSCAG/E02
      IF(NDIR.LE.0)RETURN
C
C*** Calculate the 2 matrix elements FM1 for selected directions
C F1L(NDIR) is to the polarization state e1 (E in scattering plane)
C F2L(NDIR) is to the polarization state e2 (E perp. to scatt. plane)
C IMPORTANT NOTE: The following calculation assumes that the incident
C wave is linearly polarized, with E00 having no imaginary component.
C
      TERM=AK3/SQRT(E02)
      DO 230 ND=1,NDIR
         CXF1L(ND)=(0.0,0.0)
         CXF2L(ND)=(0.0,0.0)
         DO 210 J=1,NAT
            CXSCR2(J)=(0.0,0.0)
            CXSCR3(J)=(0.0,0.0)
            DO 200 K=1, 3
               CXSCR2(J)=CXSCR2(J)+CXP(J,K)*EM1(K,ND)
               CXSCR3(J)=CXSCR3(J)+CXP(J,K)*EM2(K,ND)
  200       CONTINUE
  210    CONTINUE
         DO 220 J=1, NAT
            CXSCR1(J)=EXP(-CXI*(AKS(1,ND)*IXYZ(J,1)+
     &                          AKS(2,ND)*IXYZ(J,2)+
     &                          AKS(3,ND)*IXYZ(J,3)))
            CXF1L(ND)=CXF1L(ND)+CXSCR2(J)*CXSCR1(J)
            CXF2L(ND)=CXF2L(ND)+CXSCR3(J)*CXSCR1(J)
  220    CONTINUE
C
C   Units: Dipole polarizations are in units of [E]*d**3, where d=lattice
C          spacing, and [E]=dimensions of E field
C          Therefore, up to this point CXF1L, CXF2L are in units
C          of [E]*d**3
C          Now multiply by AK3/SQRT(E02) (units of 1/([E]d**3)) to get
C          dimensionless result.
C
         CXF1L(ND)=CXF1L(ND)*TERM
         CXF2L(ND)=CXF2L(ND)*TERM
  230 CONTINUE
      RETURN
C
      END
      SUBROUTINE SCAVEC(MXSCA,NSCAT,THETAN,PHIN,CXE01,ENSC,EM1,EM2)
C***********************************************************************
C Given:
C      MXSCA=dimension information for ENSC,PHIN,THETAN
C      NSCAT=number of scattering directions
C      THETAN(1-NSCAT)=scattering angles theta
C      PHIN(1-NSCAT)=scattering angles phi
C      CXE01(1-3)=complex polarization vector 1 (phi=0 direction
C                 is defined by Real component of CXE01)
C Returns:
C      ENSC(1-3,1-NSCAT)=scattering vectors in Lab Frame
C      EM1(1-3,1-NSCAT)=scattered pol vectors parallel to scat. plane
C                       in Lab Frame
C      EM2(1-3,1-NSCAT)=scattered pol vectors perp. to scat. plane
C                       in Lab Frame
C It is assumed that incident propagation vector is in x-direction
C in Lab Frame
C
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
C Arguments:
      INTEGER MXSCA,NSCAT
      COMPLEX CXE01(3)
      REAL EM1(3,MXSCA), EM2(3,MXSCA), ENSC(3,MXSCA), PHIN(MXSCA),
     &     THETAN(MXSCA)
C
C Local variables:
      INTEGER I
      REAL COSPHI,COSTHE,DENOM,PHI0,SINPHI,SINTHE
C
C***********************************************************************
C*** First determine orientation of phi=0 direction
      DENOM=SQRT((REAL(CXE01(1)))**2 + (REAL(CXE01(2)))**2 + 
     &           (REAL(CXE01(3)))**2 )
      PHI0=ACOS(REAL(CXE01(2))/DENOM)
      DO 1000 I=1,NSCAT
         COSPHI=COS(PHIN(I)+PHI0)
         SINPHI=SIN(PHIN(I)+PHI0)
         COSTHE=COS(THETAN(I))
         SINTHE=SIN(THETAN(I))
         ENSC(1,I)=COSTHE
         ENSC(2,I)=SINTHE*COSPHI
         ENSC(3,I)=SINTHE*SINPHI
         EM1(1,I)=-SINTHE
         EM1(2,I)=COSTHE*COSPHI
         EM1(3,I)=COSTHE*SINPHI
         EM2(1,I)=0.
         EM2(2,I)=-SINPHI
         EM2(3,I)=COSPHI
 1000 CONTINUE
      RETURN
      END
      SUBROUTINE TARCYL(A1,A2,A,B,AXIS,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
C***********************************************************************
C Purpose: to construct cylindrical target from "atoms"
C Input:
C        AXIS = 1.,2.,3. for symmetry axis along x,y,z axis
C        A = cylinder length (in units of lattice spacing d)
C        B = cylinder diameter (in units of lattice spacing d)
C        IOSHP=device number for "target.out" file
C             =-1 to suppress printing of "target.out"
C        MXNAT = dimensioning information
C Returns:
C        A1(1-3) = unit vector along cylinder axis
C        A2(1-3) = unit vector perpendicular to cylinder axis
C        CDESCR = string describing target (up to 67 chars.)
C        NAT = number of atoms in target
C        IX(1-NAT)=x/d for atoms in target
C        IY(1-NAT)=y/d
C        IZ(1-NAT)=z/d
C
C B.T.Draine, Princeton Univ. Obs., 87/4/3
C History:
C 89/12/21 (BTD): Modified for use by ddscat 89/12/21
C 90/12/03 (BTD): Modified for conformity with other target routines
C 91/01/05 (BTD): Changed I4 -> I7 when printing NAT
C 91/04/22 (BTD): Added XV,YV,ZV and A1,A2 to WRITE(IOSHP,FMT=9020)
C 93/03/12 (BTD): Changed CDESCR*(*) -> CDESCR*67
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 Scalar arguments:
      CHARACTER CDESCR*67
      INTEGER MXNAT,NAT
      REAL A,AXIS,B
C Array arguments:
      INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT)
      REAL A1(3),A2(3)
C Local scalars:
      INTEGER IAXIS,IOSHP,JX,JY,JZ,
     &    NB,NBMAX,NBMIN,NFAC,NLAY,
     &    NX1,NX2,NY1,NY2,NZ1,NZ2
      LOGICAL IN
      REAL ASPR,PI,REFF2,ROFF,X,XOFF,Y,YOFF,Z,ZOFF
C
      PI=4.*ATAN(1.)
      IAXIS=AXIS
C
C Determine limits for testing x,y,z values
C In case of disk axis, run from I=1 to I=INT(A+0.5)
C In radial directions, place atoms as follows:
C     If B is close to even number
C          x = i + 0.5
C          with i running from -int((b+.5)/2) to int((b+.5)/2)-1
C     If B is close to odd number
C          x = i
C          with i running from -int((b+.5)/2) to int((b+.5)/2)
      NB=INT(0.5*(B+0.5))
      IF(NB.EQ.2*(NB/2))THEN
         ROFF=0.5
         NBMIN=-NB
         NBMAX=NB-1
      ELSE
         ROFF=0.
         NBMIN=-NB
         NBMAX=NB
      ENDIF
C
      DO 0100 JX=1,3
         A1(JX)=0.
         A2(JX)=0.
 0100 CONTINUE
         NX1=NBMIN
         NX2=NBMAX
         NY1=NBMIN
         NY2=NBMAX
         NZ1=NBMIN
         NZ2=NBMAX
         XOFF=ROFF
         YOFF=ROFF
         ZOFF=ROFF
      IF(IAXIS.EQ.1)THEN
         A1(1)=1.
         A2(2)=1.
         XOFF=0.
         NX1=1
         NX2=INT(A+0.5)
      ELSEIF(IAXIS.EQ.2)THEN
         A1(2)=1.
         A2(3)=1.
         YOFF=0.
         NY1=1
         NY2=INT(A+0.5)
      ELSEIF(IAXIS.EQ.3)THEN
         A1(3)=1.
         A2(1)=1.
         ZOFF=0.
         NZ1=1
         NZ2=INT(A+0.5)
      ENDIF
      NAT=0
      DO 1200 JZ=NZ1,NZ2
         Z=REAL(JZ)+ZOFF
         DO 1100 JY=NY1,NY2
            Y=REAL(JY)+YOFF
            DO 1000 JX=NX1,NX2
               X=REAL(JX)+XOFF
               CALL TESTCYL(X,Y,Z,A,B,IAXIS,IN)
               IF(IN)THEN
                  NAT=NAT+1
                  IF(NAT.GT.MXNAT)THEN
                      CALL ERRMSG('FATAL','TARCYL',' NAT.GT.MXNAT')
                  ENDIF
                  IX(NAT)=JX
                  IY(NAT)=JY
                  IZ(NAT)=JZ
               ENDIF
 1000       CONTINUE
 1100   CONTINUE
 1200 CONTINUE
C
C NLAY = number of layers in cylinder
C NFAC = number of atoms in slice
      NLAY=INT(A+0.5)
      NFAC=NAT/NLAY
C REFF2=EFFECTIVE RADIUS**2 OF DISK
      REFF2=REAL(NFAC)/PI
C ASPR=ASPECT RATIO (LENGTH/DIAMETER)
      ASPR=0.5*REAL(NLAY)/SQRT(REFF2)
      WRITE(CDESCR,FMT='(A,I7,A,I4,A,I4,A,F7.4)')
     &    ' Cyl.prism, NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY,
     &    ' asp.ratio=',ASPR
C
      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)A,B,IAXIS,NAT,A1,A2
         DO 40 JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IX(JX),IY(JX),IZ(JX)
   40    CONTINUE
         CLOSE(UNIT=IOSHP)
      ENDIF
      RETURN
 6100 FORMAT(' IAXIS=',I2,/,' A=',F7.3,/,' B=',F7.3,/,' NAT=',I4)
 9020 FORMAT(' >TARELL  cylindrical target: Length=',F8.4,
     &' Diameter=',F8.4,' IAXIS=',I1,/,
     &I7,' = NAT',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &'     JA  IX  IY  IZ')
 9030 FORMAT(I7,3I4)
      END
      SUBROUTINE TESTCYL(X,Y,Z,A,B,IAXIS,IN)
C Scalar arguments:
      INTEGER IAXIS
      LOGICAL IN
      REAL A,B,X,Y,Z
C Local scalars:
      REAL R2,S
      IN=.FALSE.
      IF(IAXIS.EQ.1)THEN
         S=X
         R2=Y*Y+Z*Z
      ELSEIF(IAXIS.EQ.2)THEN
         S=Y
         R2=X*X+Z*Z
      ELSEIF(IAXIS.EQ.3)THEN
         S=Z
         R2=X*X+Y*Y
      ENDIF
      IF(R2.GT.0.25*B*B)RETURN
      IF(S.LT.1.)RETURN
      IF(S.GT.INT(A+.5))RETURN
      IN=.TRUE.
      RETURN
      END

      SUBROUTINE TARELL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
C***********************************************************************
C Routine to construct ellipsoid from "atoms"
C Input:
C        AX=(x-length)/d    (d=lattice spacing)
C        AY=(y-length)/d
C        AZ=(z-length)/d
C        IOSHP=device number for "target.out" file
C             =-1 to suppress printing of "target.out"
C        MXNAT=dimensioning information (max number of atoms)
C Output:
C        A1(1-3)=(1,0,0)=unit vector defining target axis 1 in Target Frame
C        A2(1-3)=(0,1,0)=unit vector defining target axis 2 in Target Frame
C        NAT=number of atoms in target
C        IX(1-NAT)=x/d for atoms of target
C        IY(1-NAT)=y/d
C        IZ(1-NAT)=z/d
C        CDESCR=description of target (up to 67 characters)
C B.T.Draine, Princeton Univ. Obs.
C History:
C 91/01/05 (BTD): Changed I4 -> I7 when printing NAT
C 91/04/22 (BTD): Added AX,AY,AZ and A1,A2 to WRITE(IOSHP,FMT=9020)
C 91/05/01 (BTD): TEMPORARY hack to be able to reproduce all arrays
C                 used by Draine 88 (in which XOFF=YOFF=ZOFF=0.5 always).
C 92/10/26 (BTD): Added code to set string CDESCR
C                 Added one call to WRIMSG
C 93/01/01 (BTD): Changed method for determining target axes A1, A2
C                 [previous criteria -- A1 = longest dimension,
C                 A2 = intermediate dimension -- created confusion]
C                 new rule: A1=(1,0,0), A2=(0,1,0)
C 93/03/12 (BTD): Changed CDESCR*(*) -> CDESCR*67
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
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT)
      REAL AX,AY,AZ
      REAL A1(3),A2(3)
C** Local variables:
      INTEGER JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2
      REAL AX2,AY2,AZ2,R,RYZ2,RZ2,X,XOFF,Y,YOFF,Z,ZOFF
C***********************************************************************
C
C Routine to construct pseudo-ellipsoidal target aray.
C With occupied array sites contained within ellipsoidal surface
C defined by (X/AX)**2+(Y/AY)**2+(Z/AZ)**2=0.25
C Ideal volume V=(pi/6)*AX*AY*AZ
C
C Dipoles are located at sites
C (x,y,z)=(I+XOFF,J+YOFF,Z+KOFF), I,J,K=integers
C                                 XOFF,YOFF,ZOFF=constants
C
C For sphere: call with AX=AY=AZ
C For spheroid: call with AX=AY (or AY=AZ or AX=AZ)
C B.T.Draine, Princeton Univ. Obs., 88/8/12
C
C
C Criterion for choosing XOFF:
C If AX is close to an even integer, take XOFF=1/2
C If AX is close to an odd integer, take XOFF=0
      JX=INT(AX+.5)
      IF(2*(JX/2).EQ.JX)THEN
         XOFF=0.5
      ELSE
         XOFF=0.
      ENDIF
C Similar criterion for YOFF:
      JY=INT(AY+.5)
      IF(2*(JY/2).EQ.JY)THEN
         YOFF=0.5
      ELSE
         YOFF=0.
      ENDIF
C Similar criterion for ZOFF:
      JZ=INT(AZ+.5)
      IF(2*(JZ/2).EQ.JZ)THEN
         ZOFF=0.5
      ELSE
         ZOFF=0.
      ENDIF
C
C*** temporary hack for experimental purposes:
C Note: this was commented out 93/03/11 by BTD -- cannot remember why
C       or when this experimental change was made
C      XOFF=0.5
C      YOFF=0.5
C      ZOFF=0.5
C 
      LMX1=-INT(.5*AX+.5)
      LMX2=INT(.5*AX-.25)
      LMY1=-INT(.5*AY+.5)
      LMY2=INT(.5*AY-.25)
      LMZ1=-INT(.5*AZ+.5)
      LMZ2=INT(.5*AZ-.25)
      AX2=AX*AX
      AY2=AY*AY
      AZ2=AZ*AZ
      NAT=0
C
C Specify target axes A1 and A2
C Previous convention: A1 along longest dimension
C                      A2 along intermediate dimension
C New convention: A1=(1,0,0) in target frame
C                 A2=(0,1,0) in target frame
      DO 0010 JX=1,3
         A1(JX)=0.
         A2(JX)=0.
 0010 CONTINUE
C Old code:
C      IF(AX.GE.AY.AND.AX.GE.AZ)THEN
C          A1(1)=1.
C          IF(AY.GE.AZ)THEN
C              A2(2)=1.
C          ELSE
C              A2(3)=1.
C          ENDIF
C      ELSEIF(AY.GT.AX.AND.AY.GE.AZ)THEN
C          A1(2)=1.
C          IF(AX.GE.AZ)THEN
C              A2(1)=1.
C          ELSE
C              A2(3)=1.
C          ENDIF
C      ELSEIF(AZ.GT.AX.AND.AZ.GT.AY)THEN
C          A1(3)=1.
C          IF(AX.GE.AY)THEN
C              A2(1)=1.
C          ELSE
C              A2(2)=1.
C          ENDIF
C      ENDIF
C New code:
      A1(1)=1.
      A2(2)=1.
C             
C Determine list of occupied sites.
      DO 0300 JZ=LMZ1,LMZ2
         Z=REAL(JZ)+ZOFF
         RZ2=Z*Z/AZ2
         IF(RZ2.LT.0.25)THEN
            DO 0200 JY=LMY1,LMY2
               Y=REAL(JY)+YOFF
               RYZ2=RZ2+Y*Y/AY2
               IF(RYZ2.LT.0.25)THEN
                  DO 0100 JX=LMX1,LMX2
                     X=REAL(JX)+XOFF
                     R=RYZ2+X*X/AX2
                     IF(R.LT.0.25)THEN
C Site is occupied:
                        NAT=NAT+1
                        IX(NAT)=JX
                        IY(NAT)=JY
                        IZ(NAT)=JZ
                     ENDIF
 0100              CONTINUE
               ENDIF
 0200       CONTINUE
         ENDIF
 0300 CONTINUE
C
      IF(NAT.GT.MXNAT)THEN
         CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ')
      ENDIF
C***********************************************************************
C Write target description into string CDESCR
      IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN
         WRITE(CDESCR,FMT='(A,I7,A)')' Spherical target containing',
     &      NAT,' dipoles'
      ELSE
         WRITE(CDESCR,FMT='(A,I7,A)')' Ellipsoidal target containing',
     &      NAT,' dipoles'
      ENDIF
C***********************************************************************
      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)AX,AY,AZ,NAT,XOFF,YOFF,ZOFF,A1,A2
         DO 40 JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IX(JX),IY(JX),IZ(JX)
   40    CONTINUE
         CLOSE(UNIT=IOSHP)
      ENDIF

      RETURN
 9020 FORMAT(' >TARELL  ellipsoidal grain; AX,AY,AZ=',3F8.4,/,
     &I7,3F8.4,' = NAT,XOFF,YOFF,ZOFF',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &'     JA  IX  IY  IZ')
 9030 FORMAT(I7,3I4)
      END
      SUBROUTINE TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR1,SHPAR2,SHPAR3,NAT,IX,IY,IZ,ICOMP,IDVOUT,
     &                  MXN3,NAT3)
C
C Arguments:
      CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6
      INTEGER IDVOUT,IDVSHP,IOSHP,MXNAT,NAT,MXN3,NAT3
      REAL A1(3),A2(3),SHPAR1,SHPAR2,SHPAR3
      INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT),ICOMP(MXNAT,3)
      EXTERNAL ERRMSG,REASHP,TARCYL,TARELL,TARHEX,TARREC,TARTET
C
C Local variables:
      INTEGER J,J1,L,L1,L2,L3,LC
C
C***********************************************************************
C Given:
C    CSHAPE = string describing target geometry; will determine which
C             "shape" routine is invoked.  Current options are:
C             RCTNGL, ELLIPS, CYLNDR, HEXGON, TETRAH, UNICYL, ANIELL,
C             TWOELL, TWOAEL, FRMFIL, THRELL, THRAEL
C    SHPAR1,SHPAR2,SHPAR3 = three parameters defining target geometry
C    MXNAT = limit on largest allowed value of NAT
C    MXN3  = 3*MXNAT
C    IDVOUT = 1 to print out information on target shape
C            -1 to suppress printing information on target shape
C    CFLSHP = name of target shape file if CSHAPE='FRMFIL'
C    IDVSHP = device number to use to read target shape if
C             CSHAPE='FRMFIL'
C Returns:
C    A1,A2 = two 3-vectors defining initial target orientation
C    CDESCR = string describing target
C    NAT = number of dipoles in target
C    NAT3 = 3*NAT
C    IX(1-NAT),IY(1-NAT),IZ(1-NAT)=integers fixing x,y,z coordinates of
C         occupied sites in target
C    ICOMP(1-NAT,1-3)=composition for E in x,y,z direction at each site
C
C P.J.Flatau, Colorado State University
C B.T.Draine, Princeton Univ. Observatory
C History:
C 90.11.08 (BTD): Changed DO...ENDDO to DO #...# CONTINUE
C 90.12.02 (BTD): Changed ordering of data in ICOMP
C 90.12.14 (BTD): Changed option MKRECT -> RCTNGL
C 91.05.02 (BTD): Added IDVOUT to argument list
C                 Added IDVOUT to argument list for subr. REASHP
C 91.05.23 (BTD): Added A1,A2 to arg. list for subr. REASHP
C 91.11.12 (BTD): Added IDVSHP to arg. list for TARGET and REASHP
C 92.09.21 (BTD): Added option UNIELL (uniaxial ellipsoid)
C 93.01.07 (BTD): Added option TWOELL (two touching ellipsoids,
C                 consisting of materials 1 and 2)
C 92.01.07 (BTD): Added option TWOAEL (two touching anisotropic
C                 ellipsoids, 1st with diel.fnc.1,2,3 along x,y,z
C                 2nd with diel.fnc.4,5,6 along x,y,z
C 92.01.19 (BTD): Added option THRELL (three touching anisotropic
C                 ellipsoids, 1st with diel.fnc.1,2,3 along x,y,z
C                 2nd with diel.fnc.4,5,6 along x,y,z
C                 3rd with diel.fnc.7,8,9 along x,y,z
C 93.03.12 (BTD): Changed CDESCR*60 -> CDESCR*67
C                 Moved WRITE(IDVOUT,9010) from DDSCAT to TARGET
C 93.12.14 (BTD): Corrected handling of THRELL (had assigned
C                 ICOMP=1 to 50% of sites, ICOMP=2 to other 50%)
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C********************************************************************
      IF(CSHAPE.EQ.'RCTNGL') THEN
         CALL TARREC(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C Homogeneous target:
         DO 100 J=1,NAT
            DO 95 l=1,3
               ICOMP(J,L)=1
   95       CONTINUE
  100    CONTINUE
      ELSEIF(CSHAPE.EQ.'ELLIPS')THEN
         CALL TARELL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C Homogeneous target:
         DO 110 J=1,NAT
            DO 105 L=1,3
               ICOMP(J,L)=1
  105       CONTINUE
  110    CONTINUE
      ELSEIF(CSHAPE.EQ.'CYLNDR')THEN
         CALL TARCYL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C Homogeneous target:
         DO 120 J=1,NAT
            DO 115 L=1,3
               ICOMP(J,L)=1
  115       CONTINUE
  120    CONTINUE
      ELSEIF(CSHAPE.EQ.'HEXGON')THEN
         CALL TARHEX(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C Homogeneous target:
         DO 130 J=1,NAT
            DO 125 L=1,3
               ICOMP(J,L)=1
  125       CONTINUE
  130    CONTINUE
      ELSEIF(CSHAPE.EQ.'TETRAH')THEN
         CALL TARTET(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C Homogeneous target:
         DO 140 J=1,NAT
            DO 135 L=1,3
               ICOMP(J,L)=1
  135       CONTINUE
  140    CONTINUE
      ELSEIF(CSHAPE.EQ.'UNICYL')THEN
         CALL TARCYL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C SHPAR3=1,2,3 for cylinder axis in x,y,z directions
C Assume c axis to be parallel to cylinder axis
C Assume component 1 to be dielectric const. for E par. c axis
C Assume component 2 to be dielectric const. for E perp. c axis
         LC=INT(SHPAR3)
         DO 150 J=1,NAT
            DO 145 L=1,3
                ICOMP(J,L)=2
 145        CONTINUE
            ICOMP(J,LC)=1
 150     CONTINUE
C
      ELSEIF(CSHAPE.EQ.'ANIELL')THEN
         CALL TARELL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
C Dielectric constants 1,2,3 for E along directions a1,a2,a3=a1xa2
C L1=1,2,3 for axis a1 in x,y,z directions (in Target Frame)
C L2=1,2,3 for axis a2 in x,y,z directions (in Target Frame)
C L3=1,2,3 for axis a3=a1xa2 in x,y,z directions (in Target Frame)
         L1=A1(1)+2*A1(2)+3*A1(3)+1.E-6
         L2=A2(1)+2*A2(2)+3*A2(3)+1.E-6
         L3=(A1(2)*A2(3)-A1(3)*A2(2))+2*(A1(3)*A2(1)-A1(1)*A2(3))
     &      +3*(A1(1)*A2(2)-A1(2)*A2(1))+1.E-6
         DO 160 J=1,NAT
            ICOMP(J,L1)=1
            ICOMP(J,L2)=2
            ICOMP(J,L3)=3
 160     CONTINUE
C
C TWOELL -> two touching identical ellipsoids, of materials 1 and 2
C           second ellipsoid is displaced from first in x-direction
      ELSEIF(CSHAPE.EQ.'TWOELL')THEN
         CALL TAR2EL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
         J1=NAT/2
         DO 165 J=1,J1
            ICOMP(J,1)=1
            ICOMP(J,2)=1
            ICOMP(J,3)=1
 165     CONTINUE
         DO 166 J=J1+1,NAT
            ICOMP(J,1)=2
            ICOMP(J,2)=2
            ICOMP(J,3)=2
 166     CONTINUE
C
C TWOAEL -> two touching identical ellipsoids with anisotropic
C           functions.
C           second ellipsoid is displaced from first in x-direction
C           1st ellipsoid has dielectric fcns 1,2,3 in x,y,z dirs.
C           2nd ellipsoid has dielectric fcns 4,5,6 in x,y,z dirs.
      ELSEIF(CSHAPE.EQ.'TWOAEL')THEN
         CALL TAR2EL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
         J1=NAT/2
         DO 175 J=1,J1
            ICOMP(J,1)=1
            ICOMP(J,2)=2
            ICOMP(J,3)=3
 175     CONTINUE
         DO 176 J=J1+1,NAT
            ICOMP(J,1)=4
            ICOMP(J,2)=5
            ICOMP(J,3)=6
 176     CONTINUE
C
C THRELL -> three touching identical ellipsoids, of materials 1,2,3
C           second ellipsoid is displaced from first in +x-direction
C           third ellipsoid is displaced from second in +x-direction
      ELSEIF(CSHAPE.EQ.'THRELL')THEN
         CALL TAR3EL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
         J1=NAT/3
         DO 185 J=1,J1
            ICOMP(J,1)=1
            ICOMP(J,2)=1
            ICOMP(J,3)=1
 185     CONTINUE
         DO 186 J=J1+1,2*J1
            ICOMP(J,1)=2
            ICOMP(J,2)=2
            ICOMP(J,3)=2
 186     CONTINUE
         DO 187 J=2*J1+1,NAT
            ICOMP(J,1)=3
            ICOMP(J,2)=3
            ICOMP(J,3)=3
 187     CONTINUE
C
C THRAEL -> three touching identical ellipsoids with anisotropic
C           functions.
C           second ellipsoid is displaced from first in x-direction
C           1st ellipsoid has dielectric fcns 1,2,3 in x,y,z dirs.
C           2nd ellipsoid has dielectric fcns 4,5,6 in x,y,z dirs.
C           3rd ellipsoid has dielectric fcns 7,8.9 in x,y,z dirs.
      ELSEIF(CSHAPE.EQ.'THRAEL')THEN
         CALL TAR3EL(A1,A2,SHPAR1,SHPAR2,SHPAR3,CDESCR,IOSHP,MXNAT,NAT,
     &               IX,IY,IZ)
         NAT3=3*NAT
         J1=NAT/3
         DO 195 J=1,J1
            ICOMP(J,1)=1
            ICOMP(J,2)=2
            ICOMP(J,3)=3
 195     CONTINUE
         DO 196 J=J1+1,2*J1
            ICOMP(J,1)=4
            ICOMP(J,2)=5
            ICOMP(J,3)=6
 196     CONTINUE
         DO 197 J=2*J1+1,NAT
            ICOMP(J,1)=7
            ICOMP(J,2)=8
            ICOMP(J,3)=9
 197     CONTINUE
C
C add any other specific shape routines here
C
      ELSEIF(CSHAPE.EQ.'FRMFIL') THEN
         CALL REASHP(CFLSHP,IDVOUT,IDVSHP,A1,A2,CDESCR,MXNAT,NAT,
     $                IX,IY,IZ,ICOMP,MXN3,NAT3)
      ENDIF
C Check to see that arrays are large enough
      IF(NAT.GT.MXNAT)CALL ERRMSG('FATAL','TARGET',
     &                ' MXNAT < NAT; must increase MXNAT in DDSCAT')
C Write out target description
      WRITE(IDVOUT,FMT=9010)CDESCR,NAT
C
      RETURN
 9010 FORMAT(' >TARGET:',A,/,
     &       7X,I7,' = NAT0 = number of dipoles in target')
      END
      SUBROUTINE TARHEX(A1,A2,A,B,AORI,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
C Scalar arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT)
      REAL A,AORI,B
C Array arguments:
      REAL A1(3),A2(3)
C Local scalars:
      CHARACTER CMSGNM*70
      INTEGER IAXIS,IFACE,IORI,JA,JX,JY,JZ,
     $    NA,NB,NC,NFAC,NLAY,NMX,NMY,NMZ
      LOGICAL IN
      REAL ACM,ASPR,BCM,BEFF,CCM,PI,X,Y,Z
C***********************************************************************
C Routine to construct hexagonal prism from "atoms"
C Input:
C        A=prism length/d
C        B=width of one hexagonal edge/d
C        AORI=1.,2.,3.,4.,5.,6. to specify orientation of prism
C        IOSHP=device number for "target.out" file
C             =-1 to suppress printing of "target.out"
C        MXNAT=dimensioning information (max number of atoms)
C Output:
C        A1(1-3)=unit vector defining target axis 1 (parallel to prism)
C        A2(1-3)=unit vector defining target axis 2 (perp. to A1)
C        NAT=number of atoms in target
C        IX(1-NAT)=x/d for atoms of target
C        IY(1-NAT)=y/d
C        IZ(1-NAT)=z/d
C        CDESCR=string describing target (up to 67 chars.)
C
C B.T.Draine, Princeton Univ. Obs., 87/4/3
C
C History:
C 90/12/02 (BTD): IOSHP is now output device for target.out
C                 (and is now closed before returning)
C 90/12/14 (BTD): Rewritten so that atoms are always on integral sites
C                 but CM coords may be either integral or half-integral
C 91/01/05 (BTD): Change I4 -> I7 when printing NAT.
C 91/04/22 (BTD): Added XV,YV,ZV and A1,A2 to WRITE(IOSHP,FMT=9020)
C 92/09/09 (BTD+PJF): in output, NLAY was too large by factor 2 and NFAC
C                 too small by factor of 2: now fixed.
C 93/03/12 (BTD): Changed CDESCR*(*) -> CDESCR*67
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      PI=4.*ATAN(1.)
      IORI=AORI+1.e-6
C IORI =1 for symmetry axis along X axis, face normal to Y axis.
C       2                         X                      Z
C       3                         Y                      Z
C       4                         Y                      X
C       5                         Z                      X
C       6                         Z                      Y
C Determine target axis vectors A1 and A2
C Convention: A1 will be along prism symmetry axis
C             A2 will be perp. to A1 along vertex direction
      DO 0100 JA=1,3
          A1(JA)=0.
          A2(JA)=0.
 0100 CONTINUE
      IF(IORI.EQ.1)THEN
          A1(1)=1.
          A2(3)=1.
          IAXIS=1
          IFACE=2
      ELSEIF(IORI.EQ.2)THEN
          A1(1)=1.
          A2(2)=1.
          IAXIS=1
          IFACE=3
      ELSEIF(IORI.EQ.3)THEN
          A1(2)=1.
          A2(1)=1.
          IAXIS=2
          IFACE=3
      ELSEIF(IORI.EQ.4)THEN
          A1(2)=1.
          A2(3)=1.
          IAXIS=2
          IFACE=1
      ELSEIF(IORI.EQ.5)THEN
          A1(3)=1.
          A2(2)=1.
          IAXIS=3
          IFACE=1
      ELSEIF(IORI.EQ.6)THEN
          A1(3)=1.
          A2(1)=1.
          IAXIS=3
          IFACE=2
      ENDIF
C
C IAXIS=1,2,3 for symmetry axis along x,y,z axis
C IFACE=1,2,3 for "face" normal to x,y,z axis
C A = prism length
C B = prism "diameter"
C       B/2 = length of one side of hexagon
C       B= distance between opposite vertices=max.diameter of hexagon
C       for hexagon, area=(3/8)*SQRT(3)*B**2
C       for hex prism, volume=A*area=(3/8)*SQRT(3)*A*B**2
C Now determine limits for testing x,y,z values
C Along axis, run from 1 to NA=INT(A+0.5)
C Perp. to axis, run from 1 to NB=INT(B+0.5)
      NA=INT(A+0.5)
      NB=INT(B+0.5)
      NC=INT(.5*SQRT(3.)*B+0.5)
C ACM="center" of figure in axial direction
C BCM="center" of figure in B direction
C CCM="center" of figure in C direction
      ACM=REAL(NA)/2.+0.5
      BCM=REAL(NB)/2.+0.5
      CCM=REAL(NC)/2.+0.5
      IF(IAXIS.EQ.1)THEN
          NMX=NA
          NMY=NB
          NMZ=NB
          IF(IFACE.EQ.2)NMY=NC
          IF(IFACE.EQ.3)NMZ=NC
      ELSEIF(IAXIS.EQ.2)THEN
          NMY=NA
          NMX=NB
          NMZ=NB
          IF(IFACE.EQ.1)NMX=NC
          IF(IFACE.EQ.3)NMZ=NC
      ELSEIF(IAXIS.EQ.3)THEN
          NMZ=NA
          NMX=NB
          NMY=NB
          IF(IFACE.EQ.1)NMX=NC
          IF(IFACE.EQ.2)NMY=NC
      ENDIF
      NAT=0
      write(0,*)'nmx,nmy,nmz=',nmx,nmy,nmz
      write(0,*)'acm,bcm,ccm=',acm,bcm,ccm
      write(0,*)'iaxis,iface=',iaxis,iface
      DO 1200 JZ=1,NMZ
          Z=REAL(JZ)
          DO 1100 JY=1,NMY
              Y=REAL(JY)
              DO 1000 JX=1,NMX
                  X=REAL(JX)
                  CALL TESTHEX(IAXIS,IFACE,ACM,BCM,CCM,A,B,X,Y,Z,IN)
                  IF(IN)THEN
                      NAT=NAT+1
                      IX(NAT)=JX
                      IY(NAT)=JY
                      IZ(NAT)=JZ
                  ENDIF
 1000         CONTINUE
 1100      CONTINUE
 1200 CONTINUE
C Now compute some geometric properties of pseudo-hexagonal prism
C   NLAY = number of layers = effective length of prism/d
C   NFAC = number of atoms in one hexagonal face
      NLAY=INT(A)
      NFAC=NAT/NLAY
C BEFF = effective length of hexagonal size/d
C        computed for hexagon of area NFAC*d**2
C                                     =(3/2)*SQRT(3)*BEFF**2
C Note: for hexagon, vertex-vertex diameter = 2*BEFF
      BEFF=.6204032*SQRT(FLOAT(NFAC))
C ASPR = effective aspect ratio of target = length/(2*BEFF)
      ASPR=.5*FLOAT(NLAY)/BEFF
C Note: string CDESCR will be printed by subr. TARGET
      WRITE(CDESCR,FMT='(A,I7,A)')' Hexagonal prism of NAT=',NAT,
     &                             ' dipoles'
C Here print any additional target info which is desired by using
C subroutine WRIMSG
      WRITE(CMSGNM,FMT='(A,3F7.4,A)')
     $    ' A1 = (',A1,') = hexagon symmetry axis'
      CALL WRIMSG('TARHEX',CMSGNM)
      WRITE(CMSGNM,FMT='(A,3F7.4,A)')
     $    ' A2 = (',A2,') = center-vertex axis'
      CALL WRIMSG('TARHEX',CMSGNM)
      WRITE(CMSGNM,FMT='(A,I7,A,I4,A,I4,A,F7.4)')
     $    ' NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY,' aspect ratio=',ASPR
      CALL WRIMSG('TARHEX',CMSGNM)
      WRITE(CDESCR,FMT='(A,I7,A,I4,A,I4,A,F7.4)')
     $    ' Hex prism,NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY,
     $    ' asp.ratio=',ASPR
C NOW STORE ATOM COORDINATES IN FILE
      IF(IOSHP.GT.0)THEN
          OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
C*** 3 lines of general description allowed:
          WRITE(IOSHP,9020)A,B,NAT,IAXIS,IFACE
C***
          DO 2000 JX=1,NAT
              WRITE(IOSHP,FMT=9030)JX,IX(JX),IY(JX),IZ(JX)
 2000     CONTINUE
          CLOSE(UNIT=IOSHP)
      ENDIF
      RETURN
 9020 FORMAT(' Hexagonal Prism:  A=',F7.4,' B=',F7.4,/,
     &I7,'= NAT. Orientation: IAXIS=',I1,' IFACE=',I2,/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &'     JA  IX  IY  IZ')
 9030 FORMAT(I7,3I4)
      END
      SUBROUTINE TESTHEX(IAXIS,IFACE,ACM,BCM,CCM,A,B,X,Y,Z,IN)
C Arguments:
      INTEGER IAXIS,IFACE
      REAL ACM,BCM,CCM,A,B,X,Y,Z
      LOGICAL IN
C Local variables:
      REAL AX,Q,U,V
C
C*** ACM,BCM,CCM=location of center in A,B,C directions
C    A direction is along axis
C    B direction is vertex to vertex of hexagon
C    C direction is face to face of hexagon
C
      IF(IAXIS.EQ.1)THEN
          AX=X
          IF(IFACE.EQ.2)THEN
              U=Y
              V=Z
          ELSE
              U=Z
              V=Y
          ENDIF
      ELSEIF(IAXIS.EQ.2)THEN
          AX=Y
          IF(IFACE.EQ.1)THEN
              U=X
              V=Z
          ELSE
              U=Z
              V=X
          ENDIF
      ELSEIF(IAXIS.EQ.3)THEN
          AX=Z
          IF(IFACE.EQ.1)THEN
              U=X
              V=Y
          ELSE
              U=Y
              V=X
          ENDIF
      ENDIF
      IN=.FALSE.
C*** Test along axis:
      IF(2.*ABS(AX-ACM).GT.A)RETURN
C*** Test along C direction
C .4330127=SQRT(3)/4
      IF(ABS(U-CCM).GT.0.4330127*B)RETURN
C*** Now test whether closer to CM than line bounding edge
C .5773503=1/SQRT(3)
      Q=0.5*ABS(U-CCM)+0.8660254*ABS(V-BCM)
      IF(Q.GT.0.4330127*B)RETURN
      IN=.TRUE.
      RETURN
      END
