      SUBROUTINE TAR2EL(A1,A2,AX,AY,AZ,DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,
     &                  ICOMP)
      IMPLICIT NONE
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3)
      REAL AX,AY,AZ
      REAL A1(3),A2(3),DX(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 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        DX(1-3)=x,y,z lattice spacing/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 dipoles in target
C        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for dipole locations
C        ICOMP(1-NAT,1-3)=x,y,z composition at dipole locations
C                        = 1 for locations in first ellipsoid
C                          2 for locations in second ellipsoid
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 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 97.12.26 (BTD): Added DX(3) to argument list to support nonuniform
C                 lattice spacing.
C                 *** Note: additional code modifications required
C                     to support noncubic lattice!
C 98.03.06 (BTD): If called with noncubic lattice (not yet supported)
C                 call ERRMSG and halt with fatal error.
C 98.03.07 (BTD): Modified to write DX to "target.out" for compatibiity
C                 with REASHP
C 00.11.02 (BTD): Modified to return ICOMP, and to print ICOMP to
C                 target.out
C end history
C
C Copyright (C) 1993,1996,1997,1998,2000 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C Current version of TAR2EL is restricted to cubic lattices
C
      IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN
         CALL ERRMSG('FATAL','TAR2EL',
     &               ' tar2el does not support noncubic lattice')
      ENDIF
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 JX=1,3
         A1(JX)=0.
         A2(JX)=0.
      ENDDO
      A1(1)=1.
      A2(2)=1.
C             
C Determine list of occupied sites in first ellipsoid
C
      JXMAX=LMX1
      JXMIN=LMX2
      DO JZ=LMZ1,LMZ2
         Z=REAL(JZ)+ZOFF
         RZ2=Z*Z/AZ2
         IF(RZ2.LT.0.25)THEN
            DO JY=LMY1,LMY2
               Y=REAL(JY)+YOFF
               RYZ2=RZ2+Y*Y/AY2
               IF(RYZ2.LT.0.25)THEN
                  DO 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
                        IXYZ(NAT,1)=JX
                        IXYZ(NAT,2)=JY
                        IXYZ(NAT,3)=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
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDDO

      IF(2*NAT.GT.MXNAT)THEN
          CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ')
      ENDIF

C Set composition of first ellipsoid

      DO JX=1,3
         DO JA=1,NAT
            ICOMP(JA,JX)=1
         ENDDO
      ENDDO
         

C Now create duplicate second ellipsoid

      JXMAX=JXMAX-JXMIN+1
      DO JA=1,NAT
         IXYZ(NAT+JA,1)=IXYZ(JA,1)+JXMAX
         IXYZ(NAT+JA,2)=IXYZ(JA,2)
         IXYZ(NAT+JA,3)=IXYZ(JA,3)
      ENDDO
      NAT=2*NAT

C Set composition of second ellipsoid

      DO JX=1,3
         DO JA=NAT/2+1,NAT
            ICOMP(JA,JX)=2
         ENDDO
      ENDDO

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,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &                           ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         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',/,
     &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/,
     &'     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END
      SUBROUTINE TAR2SP(A1,A2,A_1,B_1,A_2,B_2,PHI,PRINAX,DX,
     &                   CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP)
      IMPLICIT NONE
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3)
      REAL A_1,A_2,B_1,B_2,PHI,PRINAX
      REAL A1(3),A2(3),DX(3)
C** Local variables:
      CHARACTER CMSGNM*70
      INTEGER JA,JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,
     &   NAT1,NAT2
      REAL AX2,AY2,AZ2,COSPHI,SINPHI,U,V,
     &   R,RX2,RYZ2,RZ2,X,XC1,XC2,Y,YC1,YC2,Z,ZC1,ZC2
C***********************************************************************
C Routine to construct two touching spheroids from "atoms"
C    First spheroid has symmetry axis in y direction.
C    Second spheroid is displaced in x direction, with symmetry
C    axis in yz plane, in direction ey*cos(phi)+ez*sin(phi).
C    Separation between centroids=(B_1+B_2)/2
C   
C Input:
C        A_1=length/d of 1st spheroid along symm.axis  (d=lattice spacing)
C        A_2=length/d of 2nd spheroid along symm.axis
C        B_1=diameter/d of 1st spheroid
C        B_2=diameter/d of 2nd spheroid
C        PHI=angle (deg) specifying relative orientation of spheroids
C        PRINAX=0. to set A1=(1,0,0),A2=(0,1,0) in Target Framee
C              =1. to set A1,A2=principal axes of largest and second 
C                               largest moment of inertia 
C        DX(1-3)=(dx/d,dy/d,dz/d) where dx,dy,dz=lattice spacing in x,y,z
C            directions, and d=(dx*dy*dz)**(1/3)=effective lattice spacing
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        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target
C        ICOMP(1-NAT,1-3)=1 for sites within first spheroid
C                        =2 for sites within second spheroid
C        CDESCR=description of target (up to 67 characters)
C
C B.T.Draine, Princeton Univ. Obs.
C History:
C 96.01.25 (BTD): First written.
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 96.01.29 (BTD): Modified to add PRINAX to argument list, and
C                 call to PRINAXIS to determine principal axes if
C                 PRINAX=1.
C 97.12.26 (BTD): Added DX(3) to argument list to support nonuniform
C                 lattice spacing.
C                 *** Note: additional code modifications required
C                     to support noncubic lattice!
C 98.03.06 (BTD): If called with noncubic lattice (not yet supported)
C                 call ERRMSG and halt with fatal error.
C 98.03.07 (BTD): Modified to write DX to file "target.out" for
C                 compatibility with REASHP.
C 99.06.30 (BTD): Revised procedure for constructing target.
C                 With new prescription,
C                 should obtain a symmetric structure when called
C                 with A_1,B_1 = A_2,B_2 and PHI=0 or 90
C                 Note that we have not yet made modifications 
C                 necessary to support noncubic lattice.
C 00.11.02 (BTD): Modified to write ICOMP to target.out
C end history
C
C Copyright (C) 1996,1997,1998,1999,2000 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C Current version of TAR2SP is restricted to cubic lattices
C
      IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN
         CALL ERRMSG('FATAL','TAR2SP',
     &               ' tar2sp does not yet support noncubic lattice')
      ENDIF
C Routine to construct target array representing two spheroids in
C contact.
C line connecting spheroid centers is parallel to x
C symmetry axis of first spheroid is parallel to y
C symmetry axis of second spheroid is in y-z plane, at angle PHI to y
C
C First spheroid:
C Array sites contained within spheroidal surface defined by
C    ((X-XC1)/B_1)**2+((Y-YC1)/A_1)**2+((Z-ZC1)/B_1)**2=0.25
C Ideal volume V=(pi/6)*A_1*B_1*B_1
C
C Second spheroid:
C    ((X-XC2)/B_2)**2+((U-UC1)/A_1)**2+((V-VC1)/B_1)**2=0.25
C where
C    U  =  Y*cos(PHI)  +  Z*sin(PHI)
C    V  = -Y*sin(PHI)  +  Z*cos(PHI)
C   UC1 = YC1*cos(PHI) + ZC1*sin(PHI)
C   VC1 =-YC1*sin(PHI) + ZC1*cos(PHI)
C    
C Dipoles are located at sites
C (x,y,z)=(I,J,K)*d, I,J,K=integers
C
C Criterion for choosing XC1:
C Point where spheroids contact each other should be midway between
C lattice points.  Thus take XC1+0.5*B_1 = 0.5
C or XC1 = 0.5(1.-B_1)
C
C Criterion for choosing YC1 and ZC1:
C If XC1 is close to an integer, take YC1=ZC1=0
C If XC1 is close to a half-integer, take YC1=ZC1=0.5
C
C Note that other criteria are obviously possible.  For example,
C could have chosen YC1 and ZC1 to maximize number of lattice
C sites falling within spheroidal surfaces.
C
      XC1=0.5*(1.-B_1)
      IF(XC1-INT(XC1).LT.0.25.OR.XC1-INT(XC1).GT.0.75)THEN
         YC1=0.
         ZC1=0.
      ELSE
         YC1=0.5
         ZC1=0.5
      ENDIF
C
      NAT=0
C             
C Determine list of occupied sites in first spheroid
C
      LMX1=INT(XC1-.5*B_1)
      LMX2=INT(XC1+.5*B_1)
      LMY1=INT(YC1-.5*A_1)
      LMY2=INT(YC1+.5*A_1)
      LMZ1=INT(ZC1-.5*B_1)
      LMZ2=INT(ZC1+.5*B_1)
      AX2=B_1*B_1
      AY2=A_1*A_1
      AZ2=B_1*B_1
      DO 0300 JZ=LMZ1,LMZ2
         Z=REAL(JZ)-ZC1
         RZ2=Z*Z/AZ2
         IF(RZ2.LT.0.25)THEN
            DO 0200 JY=LMY1,LMY2
               Y=REAL(JY)-YC1
               RYZ2=RZ2+Y*Y/AY2
               IF(RYZ2.LT.0.25)THEN
                  DO 0100 JX=LMX1,LMX2
                     X=REAL(JX)-XC1
                     R=RYZ2+X*X/AX2
                     IF(R.LT.0.25)THEN
                        NAT=NAT+1
                        IXYZ(NAT,1)=JX
                        IXYZ(NAT,2)=JY
                        IXYZ(NAT,3)=JZ
                     ENDIF
 0100             CONTINUE
               ENDIF
 0200       CONTINUE
         ENDIF
 0300 CONTINUE
      NAT1=NAT
C
C Determine occupied sites in second spheroid
C
      XC2=XC1+.5*(B_1+B_2)
      YC2=YC1
      ZC2=ZC1
      LMX1=INT(XC2-.5*B_2)
      LMX2=INT(XC2+.5*B_2)
C
C This is not time-consuming, so no need to optimize choices of 
C LMY1,LMY2,LMZ1,LMZ2
C
      LMY1=INT(YC1-.5*MAX(A_2,B_2))
      LMY2=INT(YC1+.5*MAX(A_2,B_2))
      LMZ1=INT(ZC1-.5*MAX(A_2,B_2))
      LMZ2=INT(ZC1+.5*MAX(A_2,B_2))
      SINPHI=SIN(3.1415927*PHI/180.)
      COSPHI=COS(3.1415927*PHI/180.)
C
C transform (y,z) -> (u,v) with u=y*cos(phi)-z*sin(phi)
C                               v=y*sin(phi)+z*cos(phi)
C
      DO 0600 JX=LMX1,LMX2
         X=REAL(JX)-XC2
         RX2=(X/B_2)**2
         IF(RX2.LT.0.25)THEN
            DO 0500 JY=LMY1,LMY2
               Y=REAL(JY)-YC2
               DO 0400 JZ=LMZ1,LMZ2
                  Z=REAL(JZ)-ZC2
                  U=(COSPHI*Y-SINPHI*Z)/A_2
                  V=(SINPHI*Y+COSPHI*Z)/B_2
                  R=RX2+U**2+V**2
                  IF(R.LT.0.25)THEN
C Site is occupied:
                     NAT=NAT+1
                     IXYZ(NAT,1)=JX
                     IXYZ(NAT,2)=JY
                     IXYZ(NAT,3)=JZ
                  ENDIF
 0400          CONTINUE
 0500       CONTINUE
         ENDIF
 0600 CONTINUE
      NAT2=NAT-NAT1
      IF(NAT.GT.MXNAT)THEN
          CALL ERRMSG('FATAL','TAR2SP',' NAT.GT.MXNAT ')
      ENDIF
C
      DO JX=1,3
         DO JA=1,NAT1
            ICOMP(JA,JX)=1
         ENDDO
         DO JA=NAT1+1,NAT
            ICOMP(JA,JX)=2
         ENDDO
      ENDDO

C Specify target axes A1 and A2
C If PRINAX=0. then
C    A1=(1,0,0) in target frame
C    A2=(0,1,0) in target frame
C If PRINAX=1. then
C    A1,A2 are principal axes of largest,second largest moment of
C    inertia

      IF(PRINAX.LE.0.)THEN
         DO JX=1,3
            A1(JX)=0.
            A2(JX)=0.
         ENDDO
         A1(1)=1.
         A2(2)=1.
      ELSE
         CALL PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,DX,A1,A2)
      ENDIF

C***********************************************************************
C Write target description into string CDESCR

      WRITE(CDESCR,FMT='(A,I7,A,I7,A,I7,A)')
     &      ' Spheroids with',NAT1,' +',NAT2,'=',NAT,' dipoles'

C***********************************************************************

      CMSGNM=CDESCR
      CALL WRIMSG('TAR2SP',CMSGNM)
      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)A_1,B_1,A_2,B_2,PHI,NAT,NAT1,NAT2,A1,A2,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &                           ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         CLOSE(UNIT=IOSHP)
      ENDIF

      RETURN
 9020 FORMAT(' >TAR2SP: touching spheroids, A_1,B_2,A_2,B_2,PHI=',5F8.4,
     &/,3I7,' = NAT,NAT1,NAT2',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/,
     &'     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END
      SUBROUTINE TAR3EL(A1,A2,AX,AY,AZ,DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,
     &                  ICOMP)
      IMPLICIT NONE
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3)
      REAL AX,AY,AZ
      REAL A1(3),A2(3),DX(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 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        DX(1-3)=lattice spacing(dx,dy,dz)/d [d=(dx*dy*dz)**(1/3)]
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 dipoles in target
C        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for target dipole locations
C        ICOMP(1-NAT,1-3)=x,y,z composition for dipole locations
C                        = 1 for locations in first ellipsoid
C                          2                  second
C                          3                  third
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 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 97.12.26 (BTD): Added DX(3) to argument list to support nonuniform
C                 lattice spacing.
C                 *** Note: additional code modifications required
C                     to support noncubic lattice!
C 98.03.06 (BTD): If called with noncubic lattice (not yet supported)
C                 call ERRMSG and halt with fatal error.
C 98.03.07 (BTD): Modify to write DX to file "target.out" for
C                 compatibility with REASHP
C 00.11.02 (BTD): Add ICOMP to argument list
C                 Set ICOMP=1,2,3 for first,second,third ellipsoid
C                 Write ICOMP to target.out
C end history
C
C Copyright (C) 1993,1996,1997,1998,2000 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************

C Current version of TAR3EL is restricted to cubic lattices

      IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN
         CALL ERRMSG('FATAL','TAR3EL',
     &               ' tar3el does not support noncubic lattice')
      ENDIF

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 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 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 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

      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 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 JX=1,3
         A1(JX)=0.
         A2(JX)=0.
      ENDDO
      A1(1)=1.
      A2(2)=1.

C Determine list of occupied sites in first ellipsoid

      JXMAX=LMX1
      JXMIN=LMX2
      DO JZ=LMZ1,LMZ2
         Z=REAL(JZ)+ZOFF
         RZ2=Z*Z/AZ2
         IF(RZ2.LT.0.25)THEN
            DO JY=LMY1,LMY2
               Y=REAL(JY)+YOFF
               RYZ2=RZ2+Y*Y/AY2
               IF(RYZ2.LT.0.25)THEN
                  DO 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
                        IXYZ(NAT,1)=JX
                        IXYZ(NAT,2)=JY
                        IXYZ(NAT,3)=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
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDDO
      DO JX=1,3
         DO JA=1,NAT
            ICOMP(JA,JX)=1
         ENDDO
      ENDDO

      IF(3*NAT.GT.MXNAT)THEN
         CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ')
      ENDIF

C Now create duplicate second and third ellipsoids

      JXMAX=JXMAX-JXMIN+1
      DO JA=1,NAT
         IXYZ(NAT+JA,1)=IXYZ(JA,1)+JXMAX
         IXYZ(NAT+JA,2)=IXYZ(JA,2)
         IXYZ(NAT+JA,3)=IXYZ(JA,3)
         IXYZ(2*NAT+JA,1)=IXYZ(JA,1)+2*JXMAX
         IXYZ(2*NAT+JA,2)=IXYZ(JA,2)
         IXYZ(2*NAT+JA,3)=IXYZ(JA,3)
      ENDDO
      NAT=3*NAT
      DO JX=1,3
         DO JA=NAT/3+1,2*NAT/3
            ICOMP(JA,JX)=2
         ENDDO
         DO JA=2*NAT/3,NAT
            ICOMP(JA,JX)=3
         ENDDO
      ENDDO

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,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &                           ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         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',/,
     &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/,
     &'     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END

      SUBROUTINE TARBLOCKS(A1,A2,DX,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP,
     &                     CDESCR,MXNAT,NAT,IXYZ,ICOMP)
      IMPLICIT NONE
C Arguments:
      INTEGER BLKSIZ,IPRINAX,NBLOCKS,IOSHP,MXNAT,NAT
      INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3)
      INTEGER XYZB(3,100)
      REAL A1(3),A2(3),DX(3)
      CHARACTER CDESCR*67
C Local variables:
      INTEGER JB,JD,JJX,JJY,JJZ,JOFF,JOFFX,JOFFY,JOFFZ,JX,JY,JZ
C***********************************************************************
C
C Subroutine to construct a target out of "cubic" blocks.
C Blocks are only cubic for case of a cubic lattice DX=(1.,1.,1.)
C In case of a noncubic lattice, "blocks" are rectangular, with an equal
C number (=BLKSIZ) of dipole spacings in each direction.
C
C Given:
C DX(1-3)     =(dx/d,dy/d,dz/d) where dx,dy,dz = x,y,z lattice spacing,
C              and d = (dx*dy*dz)**(1/3) = effective lattice spacing
C NBLOCKS     =number of cubic blocks in target
C BLKSIZ      =(integer) ratio of block width/dipole spacing
C              (there will be BLKSIZ**3 dipoles per block)
C XYZB(1-3,J) =x,y,z coordinates of block a, in units of block width
C IPRINAX     =0. to return A1=(1,0,0), A2=(0,1,0)
C             =1. to return A1,A2=principal axes with largest, 2nd largest
C              moment of inertia
C IOSHP       =device number for output file "target.out"
C             =-1 to suppress generation of file "target.out"
C MXNAT       =limit on largest allowed value of NAT=number of dipoles
C
C Returns:
C NAT         =number of dipoles in target
C IXYZ(J,1-3) =x,y,z location of dipole J on lattice
C ICOMP(J,1-3)=composition for dipole J; x,y,z directions
C A1(1-3)     =principal axis A1 in target frame (normalized)
C A2(1-3)     =principal axis A2 in target frame (normalized)
C CDESCR      =string describing target
C
C where A1 = principal axis with largest moment of inertia
C       A2 = principal axis with second-largest moment of inertia
C
C B.T. Draine, Princeton Univ. Observatory, 95.12.11
C History:
C 96.01.02 (BTD) modified to compute eigenvalues and eigenvectors of
C                moment of inertia tensor, and to define vectors A1 and 
C                A2 to be along eigenvectors with largest and second
C                largest eigenvalues, respectively.
C 96.01.26 (BTD) Replaced IX,IY,IZ by IXYZ
C                Removed principal axis calculation to separate
C                subroutine PRINAXIS
C 96.01.29 (BTD) Added variable IPRINAX to allow control over definition
C                of axes A1,A2
C 96.02.23 (BTD) Remove computation of principal axes and eigenvalues
C                as this is now carried out by routine PRINAXIS.
C 97.12.26 (BTD) Add DX(3) to argument list to support nonuniform
C                lattice spacing.
C                Add DX to argument list of PRINAXIS.
C 98.03.07 (BTD) Add DX to file target.out for compatibility with
C                new version of REASHP
C End history.
C Copyright (C) 1996,1997,1998 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************

C For the moment, let us assume that blocks do not overlap.

      JOFF=1-((BLKSIZ+1)/2)
      JD=0
      DO JB=1,NBLOCKS
         JOFFX=JOFF+BLKSIZ*XYZB(1,JB)
         JOFFY=JOFF+BLKSIZ*XYZB(2,JB)
         JOFFZ=JOFF+BLKSIZ*XYZB(3,JB)
         DO JX=1,BLKSIZ
            JJX=JX+JOFFX
            DO JY=1,BLKSIZ
               JJY=JY+JOFFY
               DO JZ=1,BLKSIZ
                  JJZ=JZ+JOFFZ
                  JD=JD+1
                  IXYZ(JD,1)=JJX
                  IXYZ(JD,2)=JJY
                  IXYZ(JD,3)=JJZ

C Homogeneous target:

                  ICOMP(JD,1)=1
                  ICOMP(JD,2)=1
                  ICOMP(JD,3)=1
               ENDDO
            ENDDO
         ENDDO
      ENDDO
      NAT=JD

      IF(IPRINAX.NE.0)THEN

C Call PRINAXIS to (1) compute principal axes
C                  (2) compute eigenvalues of moment of inertia tensor

         CALL PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,DX,A1,A2)
      ELSE

C Do not compute principal axes: simply set 
C a1=(1,0,0) and a2=(0,1,0) in target frame.

         A1(1)=1.
         A1(2)=0.
         A1(3)=0.
         A2(1)=0.
         A2(2)=1.
         A2(3)=0.

      ENDIF

C***********************************************************************
C Write target description into string CDESCR
C Note: this will replace whatever CDESCR was read in from the
C 'blocks.par' file. If you want to retain that description,
C then comment out the following statement:
C      WRITE(CDESCR,FMT='(A,I7,A,I3,A)')'Target containing',
C     &      NAT,' dipoles arranged in ',NBLOCKS,' cubic blocks'
C***********************************************************************

      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)NAT,A1,A2,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &                           ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         CLOSE(UNIT=IOSHP)
      ENDIF   
      RETURN
 9020 FORMAT(' >TARBLOCKS: cubic building blocks',/,
     &I7,' = NAT',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/,
     &'     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END
      SUBROUTINE TARCEL(A1,A2,AX,AY,AZ,BX,BY,BZ,DX,
     &                  CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP)
      IMPLICIT NONE
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3)
      REAL AX,AY,AZ,BX,BY,BZ
      REAL A1(3),A2(3),DX(3)
C** Local variables:
      INTEGER JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NIN
      REAL AX2,AY2,AZ2,BX2,BY2,BZ2,R,RR,RYZ2,RZ2,X,XOFF,Y,YOFF,Z,ZOFF
C***********************************************************************
C
C Routine to construct target consisting of two materials, with outer
C surface an allipsoid of dimensions AX,AY,AZ,
C and core/mantle interface a concentric ellipsoid of dimensions BX,BY,BZ
C Input:
C        AX=(x-length)/d of outer ellipsoid  (d=lattice spacing)
C        AY=(y-length)/d "   "     "
C        AZ=(z-length)/d "   "     "
C        BX=(x-length)/d of inner ellipsoid
C        BY=(y-length)/d "   "     "
C        BZ=(z-length)/d "   "     "
C        DX(1-3)=lattice spacing (dx/d,dy/d,dz/d), d=(dx*dy*dz)**(1/3)
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        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target
C        ICOMP(1-NAT,1-3)=1 for sites within inner ellipsoid,
C                        =2 for sites between inner and outer ellipsoids
C        CDESCR=description of target (up to 67 characters)
C
C B.T.Draine, Princeton Univ. Obs.
C History:
C 94.01.26 (AP) : Adapted from TARELL by Antonio Peimbert, Princeton U.
C 94.01.27 (BTD): Modified for compatibility with TARGET, using
C                 AX,AY,AZ,BX,BY,BZ to transfer shape parameters read by
C                 subroutine REAPAR
C 94.03.21 (BTD): Corrected format statement
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 97.12.26 (BTD): Added DX(3) to argument list to support nonuniform
C                 lattice spacing.
C                 *** Note: additional code modifications required
C                     to support noncubic lattice!
C 98.03.06 (BTD): If called with noncubic lattice (not yet supported)
C                 call ERRMSG and halt with fatal error.
C 98.03.07 (BTD): Modify to write DX to file "target.out" for
C                 compatibility with REASHP
C 00.11.02 (BTD): Remove XOFF,YOFF,ZOFF from target.out
C                 cosmetic changes
C end history
C
C Copyright (C) 1994,1996,1997,1998 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C Current version of TARCEL is restricted to cubic lattices
C
      IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN
         CALL ERRMSG('FATAL','TARCEL',
     &               ' tarcel does not support noncubic lattice')
      ENDIF
C
C Inner ellipsoid is defined by
C      (X/AX)**2+(Y/AY)**2+(Z/AZ)**2 = 0.25
C Outer ellipsoid is defined by
C      (X/BX)**2+(Y/BY)**2+(Z/BZ)**2 = 0.25
C Material within inner ellipsoid is of composition 1
C Material between inner and outer ellipsoids is of composition 2
C
C Ideal volume V=(pi/6)*BX*BY*BZ
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 concentric spheres: call with AX=AY=AZ, BX=BY=BZ
C For spheroids: call with AX=AY (or AY=AZ or AX=AZ) and BX=BY (or BY=BZ
C   or BX=BZ)
C
C***********************************************************************
C
C Check that input parameters have "inner" ellipsoid smaller than
C "outer" ellipsoid:
      IF(AX.LT.BX)CALL ERRMSG('FATAL','TARCEL',' AX < BX ')
      IF(AY.LT.BY)CALL ERRMSG('FATAL','TARCEL',' AY < BY ')
      IF(AZ.LT.BZ)CALL ERRMSG('FATAL','TARCEL',' AZ < BZ ')
C Criterion for choosing XOFF: try to optimize outer ellipsoid surface
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

      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
      BX2=BX*BX
      BY2=BY*BY
      BZ2=BZ*BZ

      NAT=0
      NIN=0

C Specify target axes A1 and A2
C   A1=(1,0,0) in target frame
C   A2=(0,1,0) in target frame

      DO JX=1,3
         A1(JX)=0.
         A2(JX)=0.
      ENDDO
      A1(1)=1.
      A2(2)=1.

C Determine list of occupied sites.

      DO JZ=LMZ1,LMZ2
         Z=REAL(JZ)+ZOFF
         RZ2=Z*Z/AZ2
         IF(RZ2.LT.0.25)THEN
            DO JY=LMY1,LMY2
               Y=REAL(JY)+YOFF
               RYZ2=RZ2+Y*Y/AY2
               IF(RYZ2.LT.0.25)THEN
                  DO 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
                        RR=Z*Z/BZ2+Y*Y/BY2+X*X/BX2
                        IF(RR.LT.0.25)THEN
C within inner ellipsoid:
                           ICOMP(NAT,1)=1
                           ICOMP(NAT,2)=1
                           ICOMP(NAT,3)=1
                           NIN=NIN+1
                        ELSE
C between inner and outer ellipsoids:
                           ICOMP(NAT,1)=2
                           ICOMP(NAT,2)=2
                           ICOMP(NAT,3)=2
                        ENDIF
                        IXYZ(NAT,1)=JX
                        IXYZ(NAT,2)=JY
                        IXYZ(NAT,3)=JZ
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDDO

      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,I6,A,I7,A)')' Concentric ellipsoids with',
     &      NIN,' inner ',NAT,' total dipoles'
      ENDIF

C***********************************************************************

      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)AX,AY,AZ,BX,BY,BZ,NAT,NIN,
     &      A1,A2,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &         ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         CLOSE(UNIT=IOSHP)
      ENDIF
      RETURN
 9020 FORMAT(' >TARCEL: concentric ellipsoids; AX,AY,AZ=',3F8.4,
     &       ' BX,BY,BZ=',3F8.4,/,
     &       I7,I7,' = NAT,NIN',/,
     &       3F9.4,' = A_1 vector',/,
     &       3F9.4,' = A_2 vector',/,
     &       3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/,
     &       '     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END
      SUBROUTINE TARCYL(A1,A2,A,B,DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP)
      IMPLICIT NONE
C Scalar arguments:
      CHARACTER CDESCR*67
      INTEGER MXNAT,NAT
      REAL A,B
C Array arguments:
      INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3)
      REAL A1(3),A2(3),DX(3)
C Local scalars:
      INTEGER IOSHP,JX,JY,JZ,
     &    NB,NFAC,NLAY,
     &    NX2,NY1,NY2,NZ1,NZ2
      REAL ASPR,PI,REFF2,Y2M,YCM,Z2,ZCM
C***********************************************************************
C Purpose: to construct cylindrical target from "atoms", with cylinder
C axis along x axis.
C Input:
C        A = cylinder length (in units of lattice spacing d)
C        B = cylinder diameter (in units of lattice spacing d)
C        DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z
C                directions, and d=(dx*dy*dz)**(1/3)=effective lattice
C                spacing
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        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms in target
C        ICOMP(1-NAT,1-3)=1 (composition identifier)
C
C B.T.Draine, Princeton Univ. Obs., 87.04.03
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 95.04.07 (BTD): Corrected error in code which affected target
C                 construction for 4.5<B<5.5, 6.5<B<7.5, ...
C                 In these cases code was construction target with
C                 diameter 4, 6, ... dipoles
C                 This error was discovered by Tao Du
C 96.01.25 (BTD): Modified for conformity with DDSCAT.5a
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 97.12.26 (BTD): Added DX(1-3) to argument list to support noncubic
C                 lattice.
C 98.02.11 (BTD): Make changes to generate target for noncubic lattice.
C                 Take opportunity to simplify and streamline.
C 98.04.27 (BTD): Delete unused variables IN,Y2,ROFF,X,Z,IAXIS
C                 Declare variabels Y2M
C                 Delete IAXIS from WRITE(IOSHP,FMT=9020) statement
C 00.11.02 (BTD): Add ICOMP to argument list, return with ICOMP=1
C                 for occupied sites.
C                 write ICOMP to target.out
C end history
C Copyright (C) 1993,1995,1996,1997,1998,2000 
C               B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************

      PI=4.*ATAN(1.)

C Cylinder axis is along x axis in target frame.

      DO JX=1,3
         A1(JX)=0.
         A2(JX)=0.
      ENDDO
      A1(1)=1.
      A2(2)=1.

C A = cylinder length/d
C B = cylinder diameter/d

C Determine centroid XCM,YCM,ZCM

C If A/d_x is near an odd integer, XCM=0
C If A/d_x is near an even integer, XCM=0.5

C If B/d_y is near an odd or even integer, YCM=0 or 0.5
C If B/d_z is near an odd or even integer, ZCM=0 or 0.5

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)

      NX2=INT(A/DX(1)+0.5)
      NB=INT(B/DX(2)+0.5)
      IF(2*(NB/2).LT.NB)THEN
         YCM=0.
         NY1=-NB/2
         NY2=NB/2
      ELSE
         YCM=0.5
         NY1=-NB/2+1
         NY2=NB/2
      ENDIF
      NB=INT(B/DX(3)+0.5)
      IF(2*(NB/2).LT.NB)THEN
         ZCM=0.
         NZ1=-NB/2
         NZ2=NB/2
      ELSE
         ZCM=0.5
         NZ1=-NB/2+1
         NZ2=NB/2
      ENDIF
      
      DO JZ=NZ1,NZ2
         Z2=((JZ-ZCM)*DX(3))**2
         Y2M=((0.5*B)**2-Z2)/DX(2)**2
         DO JY=NY1,NY2
            IF((JY-YCM)**2.LT.Y2M)THEN
               DO JX=1,NX2
                  NAT=NAT+1
                  IF(NAT.GT.MXNAT)THEN
                     CALL ERRMSG('FATAL','TARCYL',' NAT.GT.MXNAT')
                  ENDIF
                  IXYZ(NAT,1)=JX
                  IXYZ(NAT,2)=JY
                  IXYZ(NAT,3)=JZ
               ENDDO
            ENDIF
         ENDDO
      ENDDO

C Set composition

      DO JZ=1,3
         DO JX=1,NAT
            ICOMP(JX,JZ)=1
         ENDDO
      ENDDO

C NLAY = number of layers in cylinder
C NFAC = number of atoms in slice

      NLAY=NX2
      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

      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)A,B,NAT,A1,A2,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &                           ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         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,/,
     &I7,' = NAT',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &3F9.6,' = (d_x,d_y,d_z)/d',/,
     &'     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END

      SUBROUTINE TARELL(A1,A2,AX,AY,AZ,DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,
     &                  ICOMP)
      IMPLICIT NONE
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3)
      REAL AX,AY,AZ
      REAL A1(3),A2(3),DX(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 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        DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z
C                directions, and d=(dx*dy*dz)**(1/3)=effective lattice
C                spacing
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        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target
C        CDESCR=description of target (up to 67 characters)
C        ICOMP(1-NAT,1-3)=1 (composition identifier)

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 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 97.11.02 (BTD): Added DX to argument list, and modified to allow use
C                 of noncubic lattice
C 98.02.11 (BTD): Modified to write DX to file "target.out"
C 00.10.19 (BTD): Corrected to reduce length of string written to
C                 CDESCR for ellipsoidal grain
C 00.11.02 (BTD): Add ICOMP to argument list, return ICOMP=1
C                 Write ICOMP to target.out
C end history
C
C Copyright (C) 1993,1996,1997,1998,2000 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************

C Routine to construct pseudo-ellipsoidal target aray.
C With occupied array sites contained within ellipsoidal surface
C defined by (X/AX*d)**2+(Y/AY*d)**2+(Z/AZ*d)**2=0.25
C Ideal volume V=(pi/6)*AX*AY*AZ*d**3
C where d = effective lattice spacing

C Dipoles are located on lattice at sites
C (x,y,z)=(I+XOFF,J+YOFF,Z+KOFF), I,J,K=integers
C                                 XOFF,YOFF,ZOFF=constants

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.08.12

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/DX(1)+.5)
      IF(2*(JX/2).EQ.JX)THEN
         XOFF=0.5
      ELSE
         XOFF=0.
      ENDIF
C Similar criterion for YOFF:
      JY=INT(AY/DX(2)+.5)
      IF(2*(JY/2).EQ.JY)THEN
         YOFF=0.5
      ELSE
         YOFF=0.
      ENDIF
C Similar criterion for ZOFF:
      JZ=INT(AZ/DX(3)+.5)
      IF(2*(JZ/2).EQ.JZ)THEN
         ZOFF=0.5
      ELSE
         ZOFF=0.
      ENDIF

      LMX1=-INT(.5*AX/DX(1)+.5)
      LMX2=INT(.5*AX/DX(1)-.25)
      LMY1=-INT(.5*AY/DX(2)+.5)
      LMY2=INT(.5*AY/DX(2)-.25)
      LMZ1=-INT(.5*AZ/DX(3)+.5)
      LMZ2=INT(.5*AZ/DX(3)-.25)
      AX2=AX*AX
      AY2=AY*AY
      AZ2=AZ*AZ
      NAT=0

C Specify target axes A1 and A2
C Convention: A1=(1,0,0) in target frame
C             A2=(0,1,0) in target frame

      DO JX=1,3
         A1(JX)=0.
         A2(JX)=0.
      ENDDO
      A1(1)=1.
      A2(2)=1.

C Determine list of occupied sites.

      DO JZ=LMZ1,LMZ2
         Z=(REAL(JZ)+ZOFF)*DX(3)
         RZ2=Z*Z/AZ2
         IF(RZ2.LT.0.25)THEN
            DO JY=LMY1,LMY2
               Y=(REAL(JY)+YOFF)*DX(2)
               RYZ2=RZ2+Y*Y/AY2
               IF(RYZ2.LT.0.25)THEN
                  DO JX=LMX1,LMX2
                     X=(REAL(JX)+XOFF)*DX(1)
                     R=RYZ2+X*X/AX2
                     IF(R.LT.0.25)THEN
C Site is occupied:
                        NAT=NAT+1
                        IXYZ(NAT,1)=JX
                        IXYZ(NAT,2)=JY
                        IXYZ(NAT,3)=JZ
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDDO

      IF(NAT.GT.MXNAT)THEN
         CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ')
      ENDIF

C Set composition

      DO JX=1,3
         DO JZ=1,NAT
            ICOMP(JZ,JX)=1
         ENDDO
      ENDDO

C***********************************************************************
C Write target description into string CDESCR
      IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN
         WRITE(CDESCR,FMT='(A,I7,A,F6.3,F6.3,F6.3,A)')
     &         ' Sphere,',NAT,' dipoles,',
     &         DX(1),DX(2),DX(3),'=x,y,z lattice spacing'
      ELSE
         WRITE(CDESCR,FMT='(A,I7,A,F6.3,F6.3,F6.3,A)')
     &         ' Ellipsoid,',NAT,' dipoles,',
     &         DX(1),DX(2),DX(3),'=x,y,z lattice spacing'
      ENDIF
C***********************************************************************
      IF(IOSHP.GE.0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         WRITE(IOSHP,FMT=9020)AX,AY,AZ,NAT,A1,A2,DX
         DO JX=1,NAT
            WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3),
     &                           ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3)
         ENDDO
         CLOSE(UNIT=IOSHP)
      ENDIF
      RETURN
 9020 FORMAT(' >TARELL  ellipsoidal grain; AX,AY,AZ=',3F8.4,/,
     &I7,' = NAT',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/,
     &'     JA  IX  IY  IZ ICOMP(x,y,z)')
 9030 FORMAT(I7,3I4,3I2)
      END
