SUBROUTINE TAR2EL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
C***********************************************************************
C Routine to construct two touching ellipsoids from "atoms"
C Input:
C        AX=(x-length of one ellipsoid)/d    (d=lattice spacing)
C        AY=(y-length of one ellipsoid)/d
C        AZ=(z-length of one ellipsoid)/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
C Note: atoms 1 - NAT/2 are in first ellipsoid
C             NAT/2+1 - NAT are in second ellipsoid
C
C B.T.Draine, Princeton Univ. Obs.
C History:
C 92/01/07 (BTD) adapted from tarell.f
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:
      CHARACTER CMSGNM*70
      INTEGER JA,JX,JXMAX,JXMIN,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NAT2
      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
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
C
      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
      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
      A1(1)=1.
      A2(2)=1.
C             
C Determine list of occupied sites in first ellipsoid
C
      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
                        IF(NAT.GT.1)THEN
                           IF(JX.GT.JXMAX)JXMAX=JX
                           IF(JX.LT.JXMIN)JXMIN=JX
                        ELSE
                           JXMAX=JX
                           JXMIN=JX
                        ENDIF
                     ENDIF
 0100             CONTINUE
               ENDIF
 0200       CONTINUE
         ENDIF
 0300 CONTINUE
C
C Now create duplicate second ellipsoid
      JXMAX=JXMAX-JXMIN+1
      DO 0400 JA=1,NAT
         IX(NAT+JA)=IX(JA)+JXMAX
         IY(NAT+JA)=IY(JA)
         IZ(NAT+JA)=IZ(JA)
 0400 CONTINUE
      NAT=2*NAT
      IF(NAT.GT.MXNAT)THEN
          CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ')
      ENDIF
C***********************************************************************
C Write target description into string CDESCR
      NAT2=NAT/2
      IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN
         WRITE(CDESCR,FMT='(A,I7,A)')' Two spheres, each containing',
     &      NAT2,' dipoles'
      ELSE
         WRITE(CDESCR,FMT='(A,I7,A)')' Two ellipsoids, each containing',
     &      NAT2,' dipoles'
      ENDIF
C***********************************************************************
      CMSGNM=CDESCR
      CALL WRIMSG('TARELL',CMSGNM)
      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 TAR3EL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
C***********************************************************************
C Routine to construct three collinear touching ellipsoids from "atoms"
C Input:
C        AX=(x-length of one ellipsoid)/d    (d=lattice spacing)
C        AY=(y-length of one ellipsoid)/d
C        AZ=(z-length of one ellipsoid)/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 Note: atoms 1       - NAT/3   are in first ellipsoid
C             NAT/3+1 - 2*NAT/3 are in second ellipsoid
C             2*NAT/3+1 - NAT   are in third ellipsoid
C Ellipsoids are displaced from one another in x-direction
C First ellipsoid is at smallest x values, third at largest
C B.T.Draine, Princeton Univ. Obs.
C
C History:
C 93/01/20 (BTD) adapted from tar2el.f
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:
      CHARACTER CMSGNM*70
      INTEGER JA,JX,JXMAX,JXMIN,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NAT2
      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
C
      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
      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
      A1(1)=1.
      A2(2)=1.
C             
C Determine list of occupied sites in first ellipsoid
C
      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
                        IF(NAT.GT.1)THEN
                           IF(JX.GT.JXMAX)JXMAX=JX
                           IF(JX.LT.JXMIN)JXMIN=JX
                        ELSE
                           JXMAX=JX
                           JXMIN=JX
                        ENDIF
                     ENDIF
 0100             CONTINUE
               ENDIF
 0200       CONTINUE
         ENDIF
 0300 CONTINUE
C
C Now create duplicate second and third ellipsoids
      JXMAX=JXMAX-JXMIN+1
      DO 0400 JA=1,NAT
         IX(NAT+JA)=IX(JA)+JXMAX
         IY(NAT+JA)=IY(JA)
         IZ(NAT+JA)=IZ(JA)
         IX(2*NAT+JA)=IX(JA)+2*JXMAX
         IY(2*NAT+JA)=IY(JA)
         IZ(2*NAT+JA)=IZ(JA)
 0400 CONTINUE
      NAT=3*NAT
      IF(NAT.GT.MXNAT)THEN
         CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ')
      ENDIF
C***********************************************************************
C Write target description into string CDESCR
      NAT2=NAT/3
      IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN
         WRITE(CDESCR,FMT='(A,I7,A)')' Three spheres, each containing',
     &      NAT2,' dipoles'
      ELSE
         WRITE(CDESCR,FMT='(A,I7,A)')' Three ellipsoids, each ',
     &      'containing',NAT2,' dipoles'
      ENDIF
C***********************************************************************
      CMSGNM=CDESCR
      CALL WRIMSG('TARELL',CMSGNM)
      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 TARREC(A1,A2,XV,YV,ZV,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
C***********************************************************************
C Routine to construct rectangular prism from "atoms"
C Input:
C        XV=(x-length)/d    (d=lattice spacing)
C        YV=(y-length)/d
C        ZV=(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)=unit vector defining target axis 1 (longest dimension)
C        A2(1-3)=unit vector defining target axis 2 (intermediate "   )
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
C B.T.Draine, Princeton Univ. Obs.
C History:
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 93/03/12 (BTD): Specified CDESCR*67, and now use it.
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:
      REAL XV,YV,ZV
      INTEGER IOSHP,MXNAT,NAT
      CHARACTER CDESCR*67
      INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT)
      REAL A1(3),A2(3)
C
C Local variables:
      INTEGER JA,JX,JY,JZ,NX,NY,NZ
      CHARACTER CMSGNM*70
C
C External subroutines:
      EXTERNAL ERRMSG
C
C Intrinsic functions:
      INTRINSIC IFIX
C***********************************************************************
      NX=IFIX(XV+0.5)
      NY=IFIX(YV+0.5)
      NZ=IFIX(ZV+0.5)
      WRITE(CMSGNM,FMT='(A,3I4)')
     &      ' Rectangular prism; NX,NY,NZ=',NX,NY,NZ
      WRITE(CDESCR,FMT='(A,3I4)')
     &      ' Rectangular prism; NX,NY,NZ=',NX,NY,NZ
C
C Specify target axes A1 and A2
C Convention: A1 will be along longest dimension
C             A2 will be along intermediate dimension
      DO 0100 JA=1,3
         A1(JA)=0.
         A2(JA)=0.
 0100 CONTINUE
      IF(NX.GE.NY.AND.NX.GE.NZ)THEN
         A1(1)=1.
         IF(NY.GE.NZ)THEN
            A2(2)=1.
         ELSE
            A2(3)=1.
         ENDIF
      ELSEIF(NY.GT.NX.AND.NY.GE.NZ)THEN
         A1(2)=1.
         IF(NX.GE.NZ)THEN
            A2(1)=1.
         ELSE
            A2(3)=1.
         ENDIF
      ELSEIF(NZ.GT.NX.AND.NZ.GT.NY)THEN
         A1(3)=1.
         IF(NX.GE.NY)THEN
            A2(1)=1.
         ELSE
            A2(2)=1.
         ENDIF
      ENDIF
C               
C Now populate lattice:
      NAT=0
      DO 30 JZ=1,NZ
         DO 20 JY=1,NY
            DO 10 JX=1,NX
               NAT=NAT+1
               IF(NAT.GT.MXNAT)THEN
                  CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT')
               ENDIF
               IX(NAT)=JX
               IY(NAT)=JY
               IZ(NAT)=JZ
   10       CONTINUE
   20    CONTINUE
   30 CONTINUE
      IF(NAT.GT.MXNAT)THEN
         CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT ')
      ENDIF
      WRITE(CMSGNM,FMT='(A, I7, A)')
     &      '  Rectangular prism of NAT=',NAT,' dipoles'
      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)XV,YV,ZV,NAT,A1,A2
         DO 40 JA=1,NAT
            WRITE(IOSHP,FMT=9030)JA,IX(JA),IY(JA),IZ(JA)
   40    CONTINUE
         CLOSE (UNIT=IOSHP)
      ENDIF
      RETURN
 9020 FORMAT (' >TARREC   rectangular prism; AX,AY,AZ=',3F8.4,/,
     &I7,' = NAT',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &'     JA  IX  IY  IZ')
 9030 FORMAT(I7,3I4)
      END
      SUBROUTINE TARTET(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ)
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 I,IMAX,IMIN,J,JMAX,JMIN,JX,K,KMAX,KMIN
      REAL FY,S,X,XOFF,Y,YMAX,YMIN,YOFF,Z,ZMAX,ZMAX0,ZOFF
C***********************************************************************
C
C Routine to construct regular tetrahedral target array
C one of the faces is parallel to the y-z plane
C one of the edges of this face is parallel to x-y plane
C
C Input:
C        AX = length of one edge of tetrahedron
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 (1,0,0) (along one axis of tetrahedron)
C        A2(1-3) = unit vector (0,1,0)
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 Occupied array sites are those within tetrahedral surface
C Size:
C    S=length of each side of tetrahedron
C    Volume = S**3/(6*sqrt(2))
C Orientation:
C    Center of mass is at origin.
C    One face is parallel to yz plane.
C    Projection onto yz plane has vertex on y axis.
C    S=length of one side (in lattice units)
C    Vertices are at
C    A=( sqrt(3/8),         0,    0)*S
C    B=(-sqrt(1/24), sqrt(1/3),    0)*S
C    C=(-sqrt(1/24),-sqrt(1/12),-1/2)*S
C    D=(-sqrt(1/24),-sqrt(1/12), 1/2)*S
C Length in x direction = S*sqrt(2/3)
C           y           = S*sqrt(3)/2
C           z           = S
C Angle(AOB)=arccos(-1/3)=109.4712 deg.
C OA=OB=OC=OD=sqrt(3/8)*S
C Occupied sites are assumed to be located at
C (X,Y,Z)=(I+XOFF,J+YOFF,K+ZOFF)
C where I,J,K are integers, and XOFF,YOFF,ZOFF are constants.
C Program sets XOFF,YOFF,ZOFF depending on choice of parameter S.
C
C Criterion for choosing XOFF:
C    Base of tetrahedron is located at x = -sqrt(1/24)*S
C    Let IMIN be value of I for this plane
C    Choose XOFF so that IMIN+XOFF = -sqrt(1/24)*S + 0.5
C    with -0.5 < XOFF < 0.5
C Criterion for choosing YOFF:
C    One edge of tetrahedron is located at y= -S/(2*sqrt(3))
C    Let JMIN be value of J for this line
C    Choose YOFF so that JMIN+YOFF = -S/sqrt(12) + 0.5
C Criterion for choosing ZOFF:
C    One edge of tetrahedron is parallel to z axis
C    Choose ZOFF in order to have number of dipoles along
C    this edge as close as possible to S
C    e.g., if S=odd integer, then take ZOFF=0 to place dipoles
C          at integral locations
C          if S=even integer, then take ZOFF=1/2 to place dipoles
C          at half-integral locations
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/ 9 (BTD): Corrected error in location of center of mass.
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
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      S=AX
C Set XOFF (and IMIN,IMAX):
      IMIN=-INT(S*SQRT(1./24.))
      XOFF=0.5-S*SQRT(1./24.)-IMIN
      IMAX=IMIN+INT(S*SQRT(2./3.)+0.5)-1
C Set YOFF (and JMIN,JMAX):
      JMIN=-INT(S/SQRT(12.))
      YOFF=0.5-S/SQRT(12.)-JMIN
      JMAX=JMIN+INT(S*SQRT(.75)+0.5)-1
C Set ZOFF (and KMIN,KMAX):
C Determine whether S is closest to even or odd integer.
C  (Temporarily let KMIN be integer which S is closest to)
      KMIN=INT(S+0.5)
C If KMIN is even, then ZOFF=0.5
C If KMIN is odd, then ZOFF=0.
      ZOFF=0.
      IF(KMIN-2*(KMIN/2).EQ.0)ZOFF=0.5
      KMIN=-INT(0.5*S+ZOFF)
      KMAX=KMIN+INT(S+0.5)-1
C
C Determine list of occupied sites.
      NAT=0
      DO 30 I=IMIN,IMAX
         X=REAL(I)+XOFF
C YMAX=largest value of Y which can occur for this X value
C YMIN=smallest value of Y which can occur for this X value
C ZMAX0=largest value of Z which can occur for this X value
         YMAX=S*SQRT(3./16.)-X/SQRT(2.)
         YMIN=-S*SQRT(3./64.)+X/SQRT(8.)
         ZMAX0=3.*S/8.-X*SQRT(3./8.)
         DO 20 J=JMIN,JMAX
            Y=REAL(J)+YOFF
            IF(Y.GE.YMIN.AND.Y.LE.YMAX)THEN
C ZMAX=largest value of Z which can occur for this (X,Y)
               FY=(Y-YMIN)/(YMAX-YMIN)
               ZMAX=(1.-FY)*ZMAX0
               DO 10 K=KMIN,KMAX
                  Z=REAL(K)+ZOFF
                  IF(ABS(Z).LE.ZMAX)THEN
C Site is occupied:
                     NAT=NAT+1
                     IF(NAT.GT.MXNAT)THEN
                        CALL ERRMSG('FATAL','TARTET',
     &                              ' NAT.GT.MXNAT ')
                     ENDIF
                     IX(NAT)=I
                     IY(NAT)=J
                     IZ(NAT)=K
                  ENDIF
  10           CONTINUE
            ENDIF
  20     CONTINUE
  30  CONTINUE
C
C*** Specify vectors A1 and A2 which will define target orientation
C    A1 is initially along x-axis
C    A2 is initially along y-axis
      A1(1)=1.
      A1(2)=0.
      A1(3)=0.
      A2(1)=0.
      A2(2)=1.
      A2(3)=0.
C
      WRITE(CDESCR,FMT='(A,I7,A)')' Tetrahedron of NAT=',NAT,' dipoles'
      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)S,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(' >TARTET  tetrahedral grain: S=',F8.3,/,
     &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 TIMEIT(CMSGTM)
C
C     timeit_vms
C
C This version of TIMEIT is for VMS systems.
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 CMSGTM* (*)
C     ..
C     .. Local Scalars ..
      REAL DTIME, OLDTIM, TIME
      INTEGER IAD
      CHARACTER CSTA*3, CMSGNM*70
C     ..
C     .. External Subroutines ..
      EXTERNAL LIB$INIT_TIMER, LIB$SHOW_TIMER, WRIMSG
C     ..
C     .. Data statements ..
      DATA CSTA/'ONE'/
C     ..
      IF (CSTA .EQ. 'ONE') THEN
          CSTA = 'TWO'
          IAD = 0
          CALL LIB$INIT_TIMER(IAD)
      ELSE IF (CSTA .EQ. 'TWO') THEN
          WRITE (CMSGNM,FMT='(A, A)') 'TIMING RESULTS OF: ', CMSGTM
          CALL WRIMSG('TIMEIT',CMSGNM)
          CSTA = 'ONE'
          CALL LIB$SHOW_TIMER(IAD,2)
      END IF
      RETURN
      END
      SUBROUTINE WRIMSG(CSUBRT,CMSGNM)
C Standard procedure for writing messages
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 CMSGNM* (*), CSUBRT* (*)
C     ..
C     .. Local Scalars ..
      INTEGER IOOUT
C     ..
C     .. External Subroutines ..
      EXTERNAL GETSET
C     ..
      CALL GETSET('GET','IOOUT',IOOUT)
      WRITE (IOOUT,FMT=9000) CSUBRT, CMSGNM
 9000 FORMAT (' >',A,' ',A)
      RETURN
      END
      SUBROUTINE XNORM3(A,XN)
C auxilary for "rotation" routines
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 ..
      REAL XN
C     ..
C     .. Array Arguments ..
      REAL A(3)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC SQRT
C     ..
      XN = SQRT(A(1)**2+A(2)**2+A(3)**2)
      RETURN
      END