      SUBROUTINE TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT,
     &                  SHPAR,NAT,IXYZ,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),SHPAR(6)
      INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3)
      EXTERNAL ERRMSG,REASHP,TARCYL,TARELL,TARHEX,TARREC,TARTET
C
C Local variables:
      INTEGER BLKSIZ,IPRINAX,J,J1,L,NBLOCKS
      INTEGER XYZB(3,100)
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, CONELL, BLOCKS,
C             DW1996, TWOSPH
C    IOSHP = device number for "target.out" file
C          = -1 to suppress printing of "target.out"
C    SHPAR(1-6) = up to 6 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    IXYZ(1-NAT,1-3)=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 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) to
C                 allow up to 6 shape parameters
C                 Added call to TARCEL to generate two concentric
C                 ellipsoids
C 94.03.21 (BTD): Changed "CONCEL" to "CONELL" for consistency with
C                 REAPAR
C 95.12.11 (BTD): Modified to handle "BLOCKS" as target option, using
C                 subroutine TARBLOCKS
C 96.01.04 (BTD): Modified to handle "DW1996" as target option
C 96.01.25 (BTD): Modified to handle "TWOSPH" as target option
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 96.01.29 (BTD): Added variable IPRINAX to control choice of axes A1,A2
C                 in option BLOCKS
C 97.06.04 (BTD): Corrected bug in handling of option "UNICYL" -- was
C                 not setting composition correctly after changing REAPAR
C                 so that no longer had choice of cylinder axis.
C                 Cylinder axis is always in x direction in target frame.
C                 Thanks for Mike Wolff for calling attention to the
C                 bug.
C end history
C
C Copyright (C) 1993,1994,1995,1996 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,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
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,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
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,SHPAR(1),SHPAR(2),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
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,SHPAR(1),SHPAR(2),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
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,SHPAR(1),CDESCR,IOSHP,MXNAT,NAT,IXYZ)
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,SHPAR(1),SHPAR(2),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
C Cylinder axis is in x direction in Target Frame.
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
         DO 150 J=1,NAT
            ICOMP(J,1)=1
            ICOMP(J,2)=2
            ICOMP(J,3)=2
 150     CONTINUE
C
      ELSEIF(CSHAPE.EQ.'ANIELL')THEN
         CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
C Dielectric constants 1,2,3 for E along directions x,y,z in Target Frame
         DO 160 J=1,NAT
            ICOMP(J,1)=1
            ICOMP(J,2)=2
            ICOMP(J,3)=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,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
         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,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
         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,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
         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,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP,
     &               MXNAT,NAT,IXYZ)
         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 CONELL -> two concentric ellipsoids with isotropic dielectric functions
C           1 (inner) and 2 (outer)
      ELSEIF(CSHAPE.EQ.'CONELL')THEN
         CALL TARCEL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5),
     &               SHPAR(6),CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP)
C
C TAR2SP -> two touching spheroids
C
      ELSEIF(CSHAPE.EQ.'TWOSPH')THEN
         CALL TAR2SP(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5),
     &               SHPAR(6),CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP)
C
C BLOCKS -> target is constructed from identical cubic blocks
C           locations of blocks, and size of blocks (lattice units)
C           is specified in file 'blocks.par'
C           file structure:
C              target description
C              IPRINAX = 0 (a1=100,a2=010) or 1 (a1,a2=principal axes)
C              N = number of blocks
C              SIZE = width of one block, in lattice units
C              X Y Z = location of block 1 (in units of block width)
C              ...
C              X Y Z =    "     "    "   N  "   "    "    "     "
C
      ELSEIF(CSHAPE.EQ.'BLOCKS')THEN
         OPEN(UNIT=IDVSHP,FILE='blocks.par',STATUS='OLD',
     &        FORM='FORMATTED')
         READ(IDVSHP,FMT='(A67)')CDESCR
         READ(IDVSHP,*)IPRINAX
         READ(IDVSHP,*)NBLOCKS
         READ(IDVSHP,*)BLKSIZ
         DO J=1,NBLOCKS
            READ(IDVSHP,*)XYZB(1,J),XYZB(2,J),XYZB(3,J)
         ENDDO
         CLOSE(IDVSHP)
         CALL TARBLOCKS(A1,A2,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP,
     &                  CDESCR,MXNAT,NAT,IXYZ,ICOMP)
C
C Homogeneous target:
C
         DO J=1,NAT
            DO L=1,3
               ICOMP(J,L)=1
            ENDDO
         ENDDO
C
C DW1996 -> the 13-cube target studied by Draine & Weingartner 1996
C           SHPAR(1)=width per block in lattice units
C
      ELSEIF(CSHAPE.EQ.'DW1996')THEN
         IPRINAX=1
         NBLOCKS=13
         BLKSIZ=SHPAR(1)
         XYZB(1,1)=0.
         XYZB(2,1)=1.
         XYZB(3,1)=0.
         XYZB(1,2)=1.
         XYZB(2,2)=1.
         XYZB(3,2)=0.
         XYZB(1,3)=0.
         XYZB(2,3)=2.
         XYZB(3,3)=0.
         XYZB(1,4)=1.
         XYZB(2,4)=2.
         XYZB(3,4)=0.
         XYZB(1,5)=0.
         XYZB(2,5)=1.
         XYZB(3,5)=1.
         XYZB(1,6)=1.
         XYZB(2,6)=1.
         XYZB(3,6)=1.
         XYZB(1,7)=0.
         XYZB(2,7)=2.
         XYZB(3,7)=1.
         XYZB(1,8)=1.
         XYZB(2,8)=2.
         XYZB(3,8)=1.
         XYZB(1,9)=2.
         XYZB(2,9)=1.
         XYZB(3,9)=0.
         XYZB(1,10)=2.
         XYZB(2,10)=2.
         XYZB(3,10)=0.
         XYZB(1,11)=0.
         XYZB(2,11)=0.
         XYZB(3,11)=1.
         XYZB(1,12)=0.
         XYZB(2,12)=0.
         XYZB(3,12)=2.
         XYZB(1,13)=0.
         XYZB(2,13)=1.
         XYZB(3,13)=2.
         CALL TARBLOCKS(A1,A2,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP,
     &                  CDESCR,MXNAT,NAT,IXYZ,ICOMP)
         CDESCR='13-cube target used by Draine & Weingartner 1996'
C Homogeneous target:
         DO J=1,NAT
            DO L=1,3
               ICOMP(J,L)=1
            ENDDO
         ENDDO
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,
     &                IXYZ,ICOMP,MXN3,NAT3)
      ELSE
         CALL ERRMSG('FATAL','TARGET',
     &               ' Uknown shape requested: check ddscat.par')
      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
C Before returning, compute:
C
      NAT3=3*NAT
      RETURN
 9010 FORMAT(' >TARGET:',A,/,
     &       7X,I7,' = NAT0 = number of dipoles in target')
      END
      SUBROUTINE TARHEX(A1,A2,A,B,CDESCR,IOSHP,MXNAT,NAT,IXYZ)
C Scalar arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 IXYZ(MXNAT,3)
      REAL A,B
C Array arguments:
      REAL A1(3),A2(3)
C Local scalars:
      CHARACTER CMSGNM*70
      INTEGER JA,JX,JY,JZ,NA,NB,NC,NFAC,NLAY,NMX,NMY,NMZ
      LOGICAL IN
      REAL ACM,ASPR,BCM,BEFF,CCM,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        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 (1,0,0) defining target axis 1 (prism axis)
C        A2(1-3)=unit vector (0,1,0) defining target axis 2 (normal to
C                one of the rectangular faces.
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=string describing target (up to 67 chars.)
C
C B.T.Draine, Princeton Univ. Obs., 87.04.03
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 95.07.13 (BTD): Simplify: restrict to case where prism axis is in
C                 x direction and y axis is normal to a rectangular
C                 face.
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C 96.02.23 (BTD): Delete PI (unused)
C 97.04.28 (BTD): Initialize NAT=0 (failure to do so reported by
C                 Michael Wolff).
C end history
C
C Copyright (C) 1993,1995,1996,1997 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      DO 0100 JA=1,3
         A1(JA)=0.
         A2(JA)=0.
 0100 CONTINUE
      A1(1)=1.
      A2(2)=1.
C
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
      NMX=NA
      NMY=NC
      NMZ=NB
      NAT=0
      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(ACM,BCM,CCM,A,B,X,Y,Z,IN)
               IF(IN)THEN
                  NAT=NAT+1
                  IXYZ(NAT,1)=JX
                  IXYZ(NAT,2)=JY
                  IXYZ(NAT,3)=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,A1,A2
C***
          DO 2000 JX=1,NAT
              WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3)
 2000     CONTINUE
          CLOSE(UNIT=IOSHP)
      ENDIF
      RETURN
 9020 FORMAT(' Hexagonal Prism:  A=',F7.4,' B=',F7.4,/,
     &I7,'= NAT. Prism axis along x, one rect. face. normal to y',/,
     &3F9.4,' = A_1 vector',/,
     &3F9.4,' = A_2 vector',/,
     &'     JA  IX  IY  IZ')
 9030 FORMAT(I7,3I4)
      END
      SUBROUTINE TESTHEX(ACM,BCM,CCM,A,B,X,Y,Z,IN)
C Arguments:
      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
      AX=X
      U=Y-CCM
      V=Z-BCM
      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).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)+0.8660254*ABS(V)
      IF(Q.GT.0.4330127*B)RETURN
      IN=.TRUE.
      RETURN
      END
      SUBROUTINE TARREC(A1,A2,XV,YV,ZV,CDESCR,IOSHP,MXNAT,NAT,IXYZ)
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 (1,0,0) defining target axis 1 in Target Frame
C        A2(1-3)=unit vector (0,1,0) 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
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 95.07.12 (BTD): Changed definition of axes A1,A2.  They used to be
C                 longest,next longest dimention.  They are now
C                 fixed to always be (1,0,0) and (0,1,0).
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C
C end history
C
C Copyright (C) 1993,1995,1996 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 IXYZ(MXNAT,3)
      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:
      INTEGER IFIX
      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 New convention: A1 will be (1,0,0) in target frame
C                 A2 will be (0,1,0) in target frame
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
      A1(1)=1.
      A2(2)=1.
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
               IXYZ(NAT,1)=JX
               IXYZ(NAT,2)=JY
               IXYZ(NAT,3)=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,IXYZ(JA,1),IXYZ(JA,2),IXYZ(JA,3)
   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,CDESCR,IOSHP,MXNAT,NAT,IXYZ)
C
C** Arguments:
      CHARACTER CDESCR*67
      INTEGER IOSHP,MXNAT,NAT
      INTEGER*2 IXYZ(MXNAT,3)
      REAL AX
      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        IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms in target
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.08 (BTD): Changed DO...ENDDO to DO #...# CONTINUE
C 90.11.09 (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 95.07.13 (BTD): Eliminated unused arguments AY,AZ
C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
C end history
C
C Copyright (C) 1993,1995,1996 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
                     IXYZ(NAT,1)=I
                     IXYZ(NAT,2)=J
                     IXYZ(NAT,3)=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,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3)
   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,DTIME)
C
C     timeit_vms
C 94.06.20 (PJF) add dtime to the formal parameters
C 96.06.19 (PJF) modified to suppress printing if CMSGTM='noprint'
C This version of TIMEIT is for VMS systems.
C
C Copyright (C) 1993,1995 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)
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            CALL LIB$SHOW_TIMER(IAD,2)
            WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE VERSION(CSTAMP)
C Arguments:
      CHARACTER CSTAMP*22
      CSTAMP='DDSCAT.5a.8 [97.7.24]'
C*******************************************************************************
C History of major changes to DDSCAT package:
C 90.12.21 (BTD): Added code to DDSCAT to create output files
C                 'qtable' (summary of Q values)
C                 'mtable' (summary of diel. func. for material 1)
C 91.08.14 (BTD): Add QBKSCA to arg. list for GETFML
C                 Modify code to print out Q_bk
C 91.08.15 (BTD): Provide different headings for qtable depending on
C                 whether IORTH=1 or 2 (add statement 9043 to DDSCAT).
C 91.09.17 (BTD): Remove calls to ALPHA and EVALA in DDSCAT (needed to 
C                 move these to subroutine GETFML since alpha from Lattice
C                 Dispersion Relation depends on propagation direction
C                 and polarization state)
C                 Add CALPHA,CXEPS,ICOMP,MXCOMP to argument list for
C                 GETFML (information required by ALPHA)
C 91.11.12 (BTD): Added variable IDVSHP to argument list for TARGET
C                 (device number for REASHP to use in reading shape file)
C 93.01.15 (BTD): Modified DDSCAT to divided output file qtable into 2 files:
C                 qtable (containing Q_ext,Q_abs,Q_sca,g,Q_bk) and
C                 qtable2 (containing Q_pha, Q_pol, Q_cpol)
C 93.01.16 (BTD): Modify DDSCAT so that qtable, qtable2, and mtable are closed
C                 after writing, and reopened for each new write
C 93.01.20 (BTD): Add MXNX, MXNY, MXNZ to argument list for subroutine
C                 EXTEND to permit checking of target size against
C                 maximum allowed dimensions
C 93.03.11 (BTD): Deleted all code associated with unused variable
C                 CMACHN (originally included to identify machine/OS)
C 93.03.12 (BTD): Changed CDESCR*60 -> CDESCR*67
C                 Moved WRITE(IDVOUT,9010) to subr. TARGET
C 93.06.02 (BTD): Corrected error in FORMAT statement 9045 in DDSCAT
C 93.06.03 (BTD): Corrected error in computation of angle-averaged
C                 backscattering cross section QBKSUM(JO) in DDSCAT (had
C                 neglected to reset sum to zero for each new target).
C 93.09.28 (BTD): Add "fix" in DDSCAT.f to work around Sun compiler/OS 
C                 bug (see description below.
C 93.12.15 (BTD): Corrected calculation of PPOL in DDSCAT.f  Added 
C                 comments on correspondence between our notation and
C                 elements of 2x2 amplitude scattering matrix and 4x4 
C                 Mueller scattering matrix.
C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) to
C                 allow up to 6 shape parameters in DDSCAT, REAPAR, 
C                 TARGET.
C 94.06.20 (PJF): Version 4c.1:
C                 Modifications to various routines to support use of
C                 GPFA FFT code of Temperton by specifying option
C                 NEWTMP in ddscat.par.  Changes made in
C                 cxfft3, eself, extend, reapar
C                 Add FFT timing code TSTFFT to distribution.
C                 Change made to timeit_xxx routines for compatibility
C                 with TSTFFT
C 94.12.20 (BTD): Add subroutine VERSION to consolidate log of
C                 significant changes to DDSCAT package.
C 95.04.07 (BTD): REASHP: Changed CDESCR*60->CDESCR*67 for compatibility
C                         with calling routine TARGET
C                 TARCYL: corrected error in code
C 95.05.15 (BTD): REASHP: Corrected READ statement
C                         added module for inhomogeneous/anisotropic
C                         targets.
C 95.05.26 (BTD): REAPAR: corrected construction of orthogonal
C                         polarization when E01 is complex (e.g.,
C                         circular polarization).
C 95.06.14 (BTD): REAPAR: added code to check that user is not requesting
C                         computation of f_ml when using other than
C                         linearly polarized E01
C 95.06.14 (BTD): Version 5a.1: first development
C                 Differs from version 4d in computing:
C                 full radiation pressure force (transverse components,
C                 in additional to longitudinal component computed
C                 previously);
C                 Torque due to scattering and absorption of radiation
C                 as measured by "vector torque cross section".
C                 Changes made in:
C                 DDSCAT
C                 GETFML (including argument list)
C                 SCAT (including argument list)
C                 VERSION (CSTAMP is now CHARACTER*22, so that a date
C                    can be attached to version number)
C 95.06.15 (BTD): Added CXE to argument list of SCAT
C                 Added new control parameters
C                 CMDSOL (to choose solution method)
C                 CMDTRQ (to choose whether or not to compute torques)
C 95.06.20 (PJF+: Major modification of DDSCAT and GETFML to support
C           BTD): use of modular iterative solvers (e.g., CCGPACK and
C                 PIM routines).
C               : add CMDSOL to argument list to specify method
C                 of solution, with "supported" options
C                 PETRKP = Petravic&Kuo-Petravic CCG routine (from 
C                          CCGPACK, coded by P.J. Flatau based on 
C                          earlier code by B.T. Draine)
C                 PBCGST = Preconditioned Biconjugate Gradient with 
C                          Stabilization (from PIM package, coded by
C                          R. Dias da Cunha and T. Hopkins)
C               : add CMDTRQ to argument list of REAPAR and GETFML
C 95.07.21 (BTD): Major revision to subroutine SCAT
C                 Added additional scratch arrays to DDSCAT,GETFML,SCAT
C 95.07.24 (BTD): Corrected errors in SCAT
C 95.08.14 (BTD): Modified GETFML and PROGRESS to include
C                 COMMON/NORMRES/ERRSCAL to allow renormalization of
C                 error |Ax-b| relative to |b|.
C 96.01.05 (BTD): Version 5a.2:
C                 Differs from version 5a.1 in including
C                 shape options BLOCKS and DW1996, including
C                 automatic calculation of principal axes A1 and A2
C                 with largest and second-largest moment of inertia
C                 eigenvalue and eigenvector problem is solved using
C                 routines BALANC, BALBAK, EIGV, ELMHS0, ELTRN0, HQR2,
C                 IPMPAR, and SPMPAR taken from the Naval Surface 
C                 Warfare Center Library of Mathematics Subroutines.
C 96.01.25 (BTD): Added shape option TWOSPH (routine TAR2SP, and
C                 modifications to TARGET,REAPAR).
C 96.02.23 (BTD): Subroutine PRINAXIS now carries out calculation of
C                 principal axes and eigenvalues of moment of inertia
C                 for arbitary grain shape (if invoked by appropriate
C                 "target generation" subroutine), and prints out
C                 3 principal axes and eigenvalues
C 96.10.17 (BTD): Version 5a.3: 
C                 Added NAT0 to argument list for routine ALPHA.
C                 This required modification of ALPHA and GETFML.
C 96.10.18 (BTD): Changed NEWTMP to GPFAFT (for Generalized Prime
C                 Factor Algorithm for Fourier Transform).
C 96.11.06 (BTD): Version 5a.4:
C                 Removed code from DDSCAT and moved to subroutines
C                 GETMUELLER and WRITESCA
C                 Changed definition of scattering angle phi.
C                 Previously, phi was measured from plane containing
C                 incident k vector and Re(CXE01)
C                 Henceforth, phi is measured from Lab x,y plane
C                 (incident radiation is along x axis).
C                 This involved changes to SCAVEC
C                 Modified GETMUELLER so that for IORTH=2 it
C                 automatically computes the complete Mueller scattering
C                 matrix for each target orientation, and for the
C                 orientational average.
C                 Modified WRITESCA so that it writes out selected
C                 elements of the Mueller matrix.
C 96.11.14 (PJF): Add timers(mxtimers), ntimers to DDSCAT and getfml
C                 timers contains information about CPU times of 
C                 various parts of the code and informations needed
C                 for timming (e.g. number of iterations)
C                 Further changes of "WRITESCA" to include binary
C                 file option. This should be of use, for example, 
C                 for users of IDL.
C                 Added "cbinflag" to DDSCA, reapar. Added
C                 "cbinfile", "iobin" and "open" unformatted file in DDSCAT
C                 "close(iobin) in DDSCAT.
C                 Removed "getset" routine. Hardwired "ioerr" in
C                 "errmsg", and "ioout" in "wrimsg". 
C                 Combined all "rotation" code in rot.2
C                 Added dd.pro package for DDSCAT output postprocessing.
C                 The package includes IDL reading code of binary
C                 DDSCAT file and Bohren and Huffman Mie routine coded in IDL.
C                 Added 10,000 lines of LAPACK code to solve 3x3
C                 eigenvalue problem. Removed NSWC l(eispack)
C                 driver because (a) it contains ald r1mach
C                 tedious to support code (b) because LAPACK is
C                 more modern. 
C                 Add several *.ps files to doc and BibTeX
C                 bibliography of recent DDA references.
C 96.11.15 (PJF): Added NetCDF binary option to "writesca.f"
C                 One needs "include" in writesca.f which is non-standard
C                 in FORTRAN77 but will probably work on all architectures.
C                 Also -lnetcdf flag is needed during the loading step
C                 if NetCDF is being used. Otherwise  "empty" NetCDF routines
C                 should be substittuted to satisfy loader (c.f. dummynetcdf.f)
C                 "ncclose" is executed from DDSCAT.f instead of writesca.f
C                 NetCDF is available from http://www.unidata.ucar.edu/
C                 for most popular computers. 
C                 Added "writebin.f" and "writecdf.f" for 5a7
C 96.11.21 (BTD): moved NCCLOS to WRITENET so that all dependence on
C                 netCDF is through subroutine WRITENET
C                 This required changes to DDSCAT,WRITESCA, and
C                 WRITENET, and adding additional call to WRITESCA
C                 from DDSCAT, and additional call to WRITENET from
C                 WRITESCA
C                 added variable IDVOUT2 to COMMON/NORMERR to permit
C                 device number for standard output to be passed to
C                 subroutine PROGRESS
C 96.12.02 (BTD): Corrected error in DDSCAT.f : SMORI was not
C                 initialized, and hence was incorrect except (possibly)
C                 for first wavelength and radius.
C                 Also eliminated some superfluous variables from
C                 DDSCAT.f
C 97.04.28 (BTD): Corrected bug in TARHEX: NAT was not initialized to 0
C 97.06.04 (BTD): Corrected bug in TARGET, affecting only target option
C                 "UNICYL" (target composition was not being set
C                 correctly).
C 97.07.24 (BTD): Corrected bug in GETMUELLER: was using PHI rather
C                 than PHIN in computation of amplitude scattering
C                 matrix and Mueller matrix from f_ml computed by
C                 GETFML.  Required change in GETMUELLER argument,
C                 list, with corresponding change required in DDSCAT.f
C                 Also corrected bugs in GETMUELLER pertaining to
C                 computation of Mueller matrix elements S_14 and
C                 S_31.  Note that because S_31 was incorrect, PPOL
C                 (evaluated in WRITESCA) was numerically incorrect.
C end history
C
C Copyright (C) 1993,1994,1995,1996,1997 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      RETURN
      END
      SUBROUTINE WRITEBIN(CBINFLAG,CBINFILE,IOBIN,CSTAMP,CDESCR,CMDFFT,
     +                    CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET,
     +                    MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH,
     +                    NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV,
     +                    NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ,
     +                    MXCOMP,BETMID,BETMXD,THTMID,THTMXD,
     +                    PHIMID,PHIMXD,ICTHM,IPHIM,TOL,WAVEA,AEFFA,
     +                    THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD,
     +                    IORI,AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR,
     +                    CXE01R,CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G,
     +                    QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,SM,PHIN,THETAN,
     +                    TIMERS)
C***********************************************************************
C Subroutine  "writebin" writes FORTRAN unformatted binary file
C to file cbinfile ('dd.bin')
C
C Original version created by P.J.Flatau, University of California, SD
C History:
C 96.11.15 (PJF): First version.
C
C end history
C Copyright (C) 1996 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C     .. Scalar Arguments ..
      REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD,
     +     THTMID,THTMXD,TOL,WAVE,XX
      INTEGER ICTHM,IOBIN,IORI,IORTH,IPHIM,IRAD,IWAV,MXBETA,MXCOMP,MXN3,
     +        MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA,MXTHET,MXTIMERS,
     +        MXWAV,NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI,NRAD,NSCAT,
     +        NTHETA,NTIMERS,NWAV,NX,NY,NZ
      CHARACTER CALPHA*6,CBINFLAG*6,CMDFFT*6,CMDTRQ*6,CSHAPE*6,
     +          CBINFILE*14,CSTAMP*22,CDESCR*67
C     .. Array Arguments ..
      COMPLEX CXE01R(3),CXE02R(3),CXEPS(MXCOMP),CXRFR(MXCOMP)
      REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),G(2),PHI(MXPHI),
     +     PHIN(MXSCA),QABS(2),QAV(17),QBKSCA(2),QEXT(2),QPHA(2),
     +     QSCAT(2),QTRQAB(3,2),QTRQSC(3,2),SM(MXSCA,4,4),THETA(MXTHET),
     +     THETAN(MXSCA),TIMERS(MXTIMERS),WAVEA(MXWAV)
      INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3)
C     .. Local Scalars ..
      INTEGER I,IENTRY,IS,J
C     .. External Subroutines ..
      EXTERNAL WRIMSG
C     .. Data statements ..
      DATA IENTRY/0/

      IF (IENTRY.EQ.0) THEN
          CALL WRIMSG('WRITEBIN','First call to binary write')
          IENTRY = 1
C General run information:
          WRITE (IOBIN) CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,CMDTRQ
C mostly dimensions:
          WRITE (IOBIN) IORTH,NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,
     +      NWAV,NRAD,NTHETA,NPHI,NBETA,NTIMERS
C maximum dimension (not really needed)
          WRITE (IOBIN) MXNX,MXNY,MXNZ,MXCOMP,MXSCA
          WRITE (IOBIN) BETMID,BETMXD,THTMID,THTMXD,PHIMID,PHIMXD
C dump scattering angles
          WRITE (IOBIN) (PHIN(I),I=1,NSCAT), (THETAN(I),I=1,NSCAT)
C misc run information
          WRITE (IOBIN) ICTHM,IPHIM,TOL
          WRITE (IOBIN) (WAVEA(I),I=1,NWAV), (AEFFA(I),I=1,NRAD),
     +      (THETA(I),I=1,NTHETA), (BETA(I),I=1,NBETA),
     +      (PHI(I),I=1,NPHI), (A1(I),I=1,3), (A2(I),I=1,3),
     +      (ICOMP(I),I=1,NAT3), (IXYZ0(I,1),I=1,NAT),
     +      (IXYZ0(I,2),I=1,NAT), (IXYZ0(I,3),I=1,NAT)

c endif  "if(cbinflag.eq.'USEBIN') then"
      END IF

C
C binary write of  for particular "iwav, irad, itheta, ibeta iphi"
C Let us try to keep only  one version of the binary file. Thus, in
C case of IORTH=1 or 'DOTORQ' some values are "missing"
C flatau -->draine  do we need smori for "ppol" do we need cxf_ij ?
C
      IF (IORI.NE.0 .AND. CBINFLAG.EQ.'ALLBIN') THEN
          WRITE (IOBIN) IWAV,IRAD,IORI, (TIMERS(I),I=1,NTIMERS)
          WRITE (IOBIN) AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR,CXE01R,
     +      CXE02R, (CXRFR(I),I=1,NCOMP), (CXEPS(I),I=1,NCOMP)
          WRITE (IOBIN) QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB,QTRQSC
C write sm as 16 independent matrices (IDL handles 7 subscripts only)
          DO 10 I = 1,4
              DO 10 J = 1,4
                  WRITE (IOBIN) (SM(IS,I,J),IS=1,NSCAT)
   10     CONTINUE
C endif "if(iori. ne. 0 .and. cbinflag.eq.'ALLCDF') then"
      END IF

      IF (IORI.EQ.0) THEN
C iori=0 case
Cflatau ---> note. I need to add code for orientational averages
C but may be I can recalculate it from sm for iwav, irad, itheta,
C ---> dump weigths
      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 History:
C 96.11.14 (PJF) Remove "getset" and hardwire" ioout"
C 96.11.20 (BTD) change IOOUT to IDVOUT
C***********************************************************************
C     .. Scalar Arguments ..
      CHARACTER CMSGNM* (*), CSUBRT* (*)
C     ..
C     .. Local Scalars ..
      INTEGER IDVOUT
      SAVE IDVOUT
C     ..
C Note: on some systems (e.g., Solaris) IDVOUT=0 generates unbuffered
C output to "standard output".  On other systems IDVOUT=0 may not be
C valid; then set IDVOUT=6 to get "standard output", probably buffered.
C
      DATA IDVOUT/0/
      WRITE (IDVOUT,FMT=9000) CSUBRT, CMSGNM
 9000 FORMAT (' >',A,' ',A)
      RETURN
      END
      SUBROUTINE WRITESCA(MXNX,MXNY,MXNZ,NX,NY,NZ,ICOMP,IXYZ0,MXNAT,
     &                    MXN3,NAT,NAT3,WAVEA,MXWAV,NWAV,AEFFA,MXRAD,
     &                    NRAD,ITHETA,IBETA,IPHI,THETA,MXTHET,NTHETA,
     &                    BETA,MXBETA,NBETA,PHI,MXPHI,NPHI,TIMERS,
     &                    MXTIMERS,NTIMERS,CBINFILE,IOBIN,CBINFLAG,
     &                    CNETFILE,IDNC,CNETFLAG,ICTHM,ILIN10,ILIN12,
     &                    IORI,IORTH,IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0,
     &                    NCOMP,NORI,NSCAT,CALPHA,CDESCR,CFLAVG,CFLSCA,
     &                    CMDFFT,CMDTRQ,CSHAPE,CSTAMP,AEFF,AK1,AKR,
     &                    BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD,
     &                    THTMID,THTMXD,TOL,WAVE,XX,A1,A2,PHIN,QABS,
     &                    QABSUM,QBKSCA,QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM,
     &                    QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB,QTRQABSUM,
     &                    QTRQSC,QTRQSCSUM,S1111,S2121,SM,SMORI,THETAN,
     &                    CX1121,CXE01,CXE01R,CXE02,CXE02R,CXEPS,CXRFR,
     &                    CXF11,CXF21,CXF12,CXF22)

C flatau note to Draine ---> there is a bug in PPOL (calculated from smori)
C Bruce --- phin looks funny. Check.
Cflatau ---> Bruce add writting qtable, qtable2, mtable in this subroutine?

C History:
C 96.11.14 (PJF) added binary write options (see also info in
C                "versions.f" about IDL postprocessing)
C                possible "chackpointing" could be done here "close(iobin)"
C                open(iobin,file=cbinfile,form='unformatted',access='append')
C                (which is non-standard Fortran 77 (but standard in F90)
C                open(iobin,file=cbinfile,form='unformatted',status='old')
C 96.11.21 (BTD) added call to WRITENET if called with CNETFLAG='NCCLOS'
C                in order to get call to NCCLOS out of DDSCAT.f
C end history
C
C     .. Scalar Arguments ..
      REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD,
     &     THTMID,THTMXD,TOL,WAVE,XX
      INTEGER IBETA,ICTHM,IDNC,ILIN10,ILIN12,IOBIN,IORI,IORTH,IPHI,
     &        IPHIM,IRAD,ITHETA,IWAV,MXBETA,MXCOMP,MXN3,MXNAT,MXNX,MXNY,
     &        MXNZ,MXPHI,MXRAD,MXSCA,MXTHET,MXTIMERS,MXWAV,NAT,NAT0,
     &        NAT3,NBETA,NCOMP,NORI,NPHI,NRAD,NSCAT,NTHETA,NTIMERS,NWAV,
     &        NX,NY,NZ
      CHARACTER CALPHA*6,CBINFLAG*6,CMDFFT*6,CMDTRQ*6,CNETFLAG*6,
     &          CSHAPE*6,CFLAVG*13,CBINFILE*14,CFLSCA*14,CNETFILE*14,
     &          CSTAMP*22,CDESCR*67
C     .. Array Arguments ..
      COMPLEX CX1121(MXSCA),CXE01(3),CXE01R(3),CXE02(3),CXE02R(3),
     &        CXEPS(MXCOMP),CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),
     &        CXF22(MXSCA),CXRFR(MXCOMP)
      REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),PHI(MXPHI),
     &     PHIN(MXSCA),QABS(2),QABSUM(2),QBKSCA(2),QBKSUM(2),QEXSUM(2),
     &     QEXT(2),QPHA(2),QPHSUM(2),QSCAG(3,2),QSCAT(2),QSCGSUM(3,2),
     &     QSCSUM(2),QTRQAB(3,2),QTRQABSUM(3,2),QTRQSC(3,2),
     &     QTRQSCSUM(3,2),S1111(MXSCA),S2121(MXSCA),SM(MXSCA,4,4),
     &     SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA),
     &     TIMERS(MXTIMERS),WAVEA(MXWAV)
      INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3)
C     .. Local Scalars ..
      COMPLEX CXVAR1
      REAL CX01R,DEGRAD,PHIND,PPOL,RVAR1,RVAR2,RVAR3,THETND
      INTEGER I,J,ND
C     .. Local Arrays ..
      REAL G(2),QAV(17)
C     .. External Subroutines ..
      EXTERNAL NAMER,WRITEBIN,WRITENET
C     .. Intrinsic Functions ..
      INTRINSIC CONJG,SQRT
C     .. Save statement ..
      SAVE DEGRAD
C     .. Data statements ..
      DATA DEGRAD/57.29578/
C
C If CNETFLAG = 'NCCLOS' then call WRITESCA to
C                             call WRITENET to
C                             call NCCLOS to close netCDF files
      IF(CNETFLAG.EQ.'NCCLOS')THEN
         CALL WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT,
     &                 CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET,MXBETA,
     &                 MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH,NX,NY,NZ,
     &                 NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV,NRAD,NTHETA,
     &                 NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ,MXCOMP,
     &                 BETMID,BETMXD,THTMID,THTMXD,PHIMID,PHIMXD,ICTHM,
     &                 IPHIM,TOL,WAVEA,AEFFA,THETA,BETA,PHI,A1,A2,
     &                 ICOMP,IXYZ0,IWAV,IRAD,ITHETA,IBETA,IPHI,IORI,
     &                 AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR,CXE01R,
     &                 CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G,QBKSCA,
     &                 QPHA,QAV,QTRQAB,QTRQSC,SM,SMORI,PHIN,THETAN,
     &                 TIMERS)
         RETURN
      ENDIF
C If IORI > 0, then write out scattering properties for this specific
C target orientation.
      IF(IORI.GT.0)THEN
C*** Compute g=<cos(theta)>
          DO 12 J=1,IORTH
              G(J)=QSCAG(1,J)/QSCAT(J)
   12     CONTINUE
          CALL NAMER(IWAV,IRAD,IORI,CFLSCA,CFLAVG)
          OPEN(UNIT=8,FILE=CFLSCA,STATUS='UNKNOWN')
          WRITE(8,FMT=9030)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0
          WRITE(8,FMT=9032)AEFF,WAVE,XX
          WRITE(8,FMT=9031)(CXRFR(ND),CXEPS(ND),ND,ND=1,NCOMP)
          WRITE(8,FMT=9033)TOL,ICTHM,IPHIM,A1,A2
          WRITE(8,FMT=9035)AKR,CXE01R,CXE02R
          WRITE(8,FMT=9040)BETAD,THETAD,PHID
          WRITE(8,FMT=9050)QEXT(1),QABS(1),QSCAT(1),G(1),QBKSCA(1),
     &      QPHA(1)
          IF(IORTH.EQ.1)THEN
Cflatau  assign "missing value" value if there is only one polarization
              DO 13 I=1,17
                  QAV(I)=-999.
   13         CONTINUE
              QEXT(2)=-999.
              QABS(2)=-999.
              QSCAT(2)=-999.
              QBKSCA(2)=-999.
              QPHA(2)=-999.
              G(2)=-999.
              DO 14 J=1,3
                  QSCAG(J,2)=-999.
   14         CONTINUE
C endif "IF(IORTH.EQ.1)THEN"
          END IF


          IF(IORTH.EQ.2)THEN
              QAV(1)=.5*(QEXT(1)+QEXT(2))
              QAV(2)=.5*(QABS(1)+QABS(2))
              QAV(3)=.5*(QSCAT(1)+QSCAT(2))
              QAV(4)=.5*(QSCAG(1,1)+QSCAG(1,2))/QAV(3)
              QAV(5)=.5*(QBKSCA(1)+QBKSCA(2))
              QAV(6)=.5*(QPHA(1)+QPHA(2))
              QAV(7)=QEXT(1)- QEXT(2)
              QAV(8)=QPHA(1)- QPHA(2)
              DO 16 J=1,3
                  QAV(8+J)=.5*(QSCAG(J,1)+QSCAG(J,2))
   16         CONTINUE
              IF(CMDTRQ.EQ.'DOTORQ')THEN
                  DO 17 J=1,3
                      QAV(11+J)=.5*(QTRQAB(J,1)+QTRQAB(J,2))
                      QAV(14+J)=.5*(QTRQSC(J,1)+QTRQSC(J,2))
   17             CONTINUE
              ENDIF

              WRITE(8,FMT=9055)QEXT(2),QABS(2),QSCAT(2),G(2),
     &          QBKSCA(2),QPHA(2),(QAV(J),J=1,8)
          ENDIF

          WRITE(8,FMT=9150)(QSCAG(J,1),J=1,3)
          IF(IORTH.EQ.2)THEN
              WRITE(8,FMT=9155)(QSCAG(J,2),J=1,3),(QAV(J),J=9,11)
          ENDIF

          IF(CMDTRQ.EQ.'DOTORQ')THEN
              WRITE(8,FMT=9250)(QTRQAB(J,1),J=1,3),
     &          (QTRQSC(J,1),J=1,3)
              IF(IORTH.EQ.2)WRITE(8,FMT=9255)(QTRQAB(J,2),J=1,3),
     &            (QTRQSC(J,2),J=1,3),(QAV(J),J=12,17)

          ELSE
CFLATAU --->DRAINE is j=1 set for IORTH=1 ???? in QTRQAB, QTRQSC ??
C add missing value code (set to -999.)
              DO 19 I=1,3
                  DO 21 J=1,2
                      QTRQAB(I,J)=-999.
                      QTRQSC(I,J)=-999.
   21             CONTINUE
   19         CONTINUE
C endif "IF(CMDTRQ.EQ.'DOTORQ')THEN"
          ENDIF

          IF(IORTH.EQ.1)THEN
              WRITE(8,FMT=9060)
              DO 22 ND=1,NSCAT
C convert scattering angles to degrees
                  PHIND=DEGRAD*PHIN(ND)
                  THETND=DEGRAD*THETAN(ND)
                  RVAR1=CONJG(CXF11(ND))*CXF11(ND)
                  RVAR2=CONJG(CXF21(ND))*CXF21(ND)
                  CXVAR1=CONJG(CXF11(ND))*CXF21(ND)
                  WRITE(8,FMT=9070)ND,THETND,PHIND,RVAR1,RVAR2,CXVAR1
   22         CONTINUE

          ELSE IF(IORTH.EQ.2)THEN
              WRITE(8,FMT=9080)
              DO 27 ND=1,NSCAT
C convert scattering angles to degrees
                  PHIND=DEGRAD*PHIN(ND)
                  THETND=DEGRAD*THETAN(ND)
C
C Compute PPOL = degree of linear polarization of scattered light
C for incident unpolarized light
C
                  PPOL=SQRT(SM(ND,2,1)**2+SM(ND,3,1)**2)/
     &                   SM(ND,1,1)
                  WRITE(8,FMT=9090)THETND,PHIND,SM(ND,1,1),SM(ND,2,1),
     &              SM(ND,3,1),SM(ND,4,1),SM(ND,1,2),SM(ND,1,3),PPOL
   27         CONTINUE
C endif "IF(IORTH.EQ.1)THEN"
          ENDIF

          CLOSE(8)
C
      ELSE
C This module writes out results from orientational averaging.
C (arrive here when called with IORI=0)
C
          OPEN(UNIT=8,FILE=CFLAVG,STATUS='UNKNOWN')
          WRITE(8,FMT=9030)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0
          WRITE(8,FMT=9032)AEFF,WAVE,XX
          WRITE(8,FMT=9031)(CXRFR(ND),CXEPS(ND),ND,ND=1,NCOMP)
          WRITE(8,FMT=9033)TOL,ICTHM,IPHIM,A1,A2
          RVAR1=AK1
          RVAR2=0.
          RVAR3=0.
          WRITE(8,FMT=9037)RVAR1,RVAR2,RVAR3,CXE01,CXE02
          WRITE(8,FMT=9042)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA,
     &      PHIMID,PHIMXD,NPHI,NORI,IORTH
C
          RVAR1=QSCGSUM(1,1)/QSCSUM(1)
          WRITE(8,FMT=9050)QEXSUM(1),QABSUM(1),QSCSUM(1),RVAR1,
     &      QBKSUM(1),QPHSUM(1)
          IF(IORTH.EQ.1)THEN
              WRITE(8,FMT=9150)(QSCGSUM(J,1),J=1,3)
              IF(CMDTRQ.EQ.'DOTORQ')WRITE(8,
     &            FMT=9250)(QTRQABSUM(J,1),J=1,3),
     &            (QTRQSCSUM(J,1),J=1,3)
              WRITE(8,FMT=9060)
              DO 28 ND=1,NSCAT
C Convert scattering angles to degrees
                  PHIND=DEGRAD*PHIN(ND)
                  THETND=DEGRAD*THETAN(ND)
                  WRITE(8,FMT=9070)ND,THETND,PHIND,S1111(ND),
     &              S2121(ND),CX1121(ND)
   28         CONTINUE
C
C*** Write orientationally-averaged Q values (except Qpha) to 'qtable':
              OPEN(UNIT=10,FILE='qtable',STATUS='OLD')
              DO 31 ND=1,ILIN10
                  READ(10,FMT=*)
   31         CONTINUE
              RVAR1=QSCGSUM(1,1)/QSCSUM(1)
              WRITE(10,FMT=9056)AEFF,WAVE,QEXSUM(1),QABSUM(1),
     &          QSCSUM(1),RVAR1,QBKSUM(1)
              ILIN10=ILIN10+1
              CLOSE(10)
              OPEN(UNIT=12,FILE='qtable2',STATUS='OLD')
              DO 32 ND=1,ILIN12
                  READ(12,FMT=*)
   32         CONTINUE
              WRITE(12,FMT=9057)AEFF,WAVE,QPHSUM(1)
              ILIN12=ILIN12+1
              CLOSE(12)

          ELSE IF(IORTH.EQ.2)THEN
              RVAR1=QSCGSUM(1,2)/QSCSUM(2)
              QAV(1)=.5*(QEXSUM(1)+QEXSUM(2))
              QAV(2)=.5*(QABSUM(1)+QABSUM(2))
              QAV(3)=.5*(QSCSUM(1)+QSCSUM(2))
              QAV(4)=.5*(QSCGSUM(1,1)+QSCGSUM(1,2))/QAV(3)
              QAV(5)=.5*(QBKSUM(1)+QBKSUM(2))
              QAV(6)=.5*(QPHSUM(1)+QPHSUM(2))
              QAV(7)=QEXSUM(1)- QEXSUM(2)
              QAV(8)=QPHSUM(1)- QPHSUM(2)
              DO 38 J=1,3
                  QAV(8+J)=.5*(QSCGSUM(J,1)+QSCGSUM(J,2))
   38         CONTINUE
              IF(CMDTRQ.EQ.'DOTORQ')THEN
                  DO 41 J=1,3
                      QAV(11+J)=.5*(QTRQABSUM(J,1)+QTRQABSUM(J,2))
                      QAV(14+J)=.5*(QTRQSCSUM(J,1)+QTRQSCSUM(J,2))
   41             CONTINUE
              ENDIF

              WRITE(8,FMT=9055)QEXSUM(2),QABSUM(2),QSCSUM(2),RVAR1,
     &          QBKSUM(2),QPHSUM(2),(QAV(J),J=1,8)
              WRITE(8,FMT=9150)(QSCGSUM(J,1),J=1,3)
              WRITE(8,FMT=9155)(QSCGSUM(J,2),J=1,3),(QAV(J),J=9,11)
              IF(CMDTRQ.EQ.'DOTORQ')THEN
                  WRITE(8,FMT=9250)(QTRQABSUM(J,1),J=1,3),
     &              (QTRQSCSUM(J,1),J=1,3)
                  WRITE(8,FMT=9255)(QTRQABSUM(J,2),J=1,3),
     &              (QTRQSCSUM(J,2),J=1,3),(QAV(J),J=12,17)
              ENDIF
C
C*** Write orientationally-averaged Q values (except Qpha) to 'qtable':
              OPEN(UNIT=10,FILE='qtable',STATUS='OLD')
              DO 46 ND=1,ILIN10
                  READ(10,FMT=*)
   46         CONTINUE
              WRITE(10,FMT=9056)AEFF,WAVE,QAV(1),QAV(2),QAV(3),QAV(4),
     &          QAV(5)
              ILIN10=ILIN10+1
              CLOSE(10)
              OPEN(UNIT=12,FILE='qtable2',STATUS='OLD')
              DO 47 ND=1,ILIN12
                  READ(12,FMT=*)
   47         CONTINUE
              WRITE(12,FMT=9057)AEFF,WAVE,QAV(6),QAV(7),QAV(8)
              ILIN12=ILIN12+1
              CLOSE(12)
C
              WRITE(8,FMT=9080)
              DO 48 ND=1,NSCAT
C Convert scattering angles to degrees
                  PHIND=DEGRAD*PHIN(ND)
                  THETND=DEGRAD*THETAN(ND)
C Compute PPOL = degree of linear polarization of scattered light
C for incident unpolarized light (Bohren & Huffman p. 382)
                  PPOL=SQRT(SMORI(ND,2,1)**2+SMORI(ND,3,1)**2)/
     &                   SMORI(ND,1,1)
                  WRITE(8,FMT=9090)THETND,PHIND,SMORI(ND,1,1),
     &              SMORI(ND,2,1),SMORI(ND,3,1),SMORI(ND,4,1),
     &              SMORI(ND,1,2),SMORI(ND,1,3),PPOL

   48         CONTINUE

          ENDIF

          CLOSE(8)
      ENDIF

C write FORTRAN unformatted file
      IF(CBINFLAG.NE.'NOTBIN')THEN
          CALL WRITEBIN(CBINFLAG,CBINFILE,IOBIN,CSTAMP,CDESCR,CMDFFT,
     &                  CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET,MXBETA,
     &                  MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH,NX,NY,NZ,
     &                  NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV,NRAD,NTHETA,
     &                  NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ,MXCOMP,
     &                  BETMID,BETMXD,THTMID,THTMXD,PHIMID,PHIMXD,ICTHM,
     &                  IPHIM,TOL,WAVEA,AEFFA,THETA,BETA,PHI,A1,A2,
     &                  ICOMP,IXYZ0,IWAV,IRAD,IORI,AEFF,WAVE,XX,AK1,
     &                  BETAD,THETAD,PHID,AKR,CX01R,CXE02R,CXRFR,CXEPS,
     &                  QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,
     &                  SM,PHIN,THETAN,TIMERS)
C endif "if(cbinflag.ne.'NOTCDF') then"
      ENDIF

      IF(CNETFLAG.NE.'NOTCDF')THEN
          CALL WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT,
     &                  CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET,MXBETA,
     &                  MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH,NX,NY,NZ,
     &                  NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV,NRAD,NTHETA,
     &                  NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ,MXCOMP,
     &                  BETMID,BETMXD,THTMID,THTMXD,PHIMID,PHIMXD,ICTHM,
     &                  IPHIM,TOL,WAVEA,AEFFA,THETA,BETA,PHI,A1,A2,
     &                  ICOMP,IXYZ0,IWAV,IRAD,ITHETA,IBETA,IPHI,IORI,
     &                  AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR,CXE01R,
     &                  CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G,QBKSCA,
     &                  QPHA,QAV,QTRQAB,QTRQSC,SM,SMORI,PHIN,THETAN,
     &                  TIMERS)
C endif "if(cnetflag.ne.'NOTCDF') then "
      ENDIF

      RETURN

 9030 FORMAT (' DDSCAT --- ',A,/,' TARGET ---',A,/,' ',A,
     &       ' --- method of solution ',/,' ',A,
     &       ' --- prescription for polarizabilies',/,' ',A,
     &       ' --- shape ',/,' NAT0 = ',I7,' = number of dipoles')
 9031 FORMAT (' n= (',F7.4,' , ',F7.4,'),  epsilon= (',F8.4,' , ',F7.4,
     &       ')  for material',I2)
 9032 FORMAT ('  AEFF=',F10.5,' =',' effective radius (physical units)',
     &       /,'  WAVE=',F10.5,' = wavelength (physical units)',/,
     &       '   K*A=',F10.5,' = 2*pi*aeff/lambda')
 9033 FORMAT ('   TOL=',1P,E10.3,' = error tolerance for CCG method',/,
     &       ' ICTHM=',I3,' = theta values used in comp. of Qsca,g',/,
     &       ' IPHIM=',I3,' = phi values used in comp. of Qsca,g',/,1X,
     &       0P,'(',F8.5,2F9.5,') = target axis A1 in Target Frame',/,
     &       1X,'(',F8.5,2F9.5,') = target axis A2 in Target Frame')
 9035 FORMAT (1X,'(',F8.5,2F9.5,') = k vector (latt. units) in TF',/,1X,
     &       3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 1 in TF',/,1X,
     &       3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 2 in TF')
 9037 FORMAT (1X,'(',F8.5,2F9.5,
     &       ') = k vector (latt. units) in Lab Frame',/,1X,
     &       3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 1 in LF',/,1X,
     &       3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 2 in LF')
 9040 FORMAT (' BETA =',F7.3,' = rotation of target around A1',/,
     &       ' THETA=',F7.3,' = angle between A1 and k',/,'  PHI =',
     &       F7.3,' = rotation of A1 around k')
 9042 FORMAT (2F8.3,' = beta_min, beta_max ;  NBETA =',I2,/,2F8.3,
     &       ' = theta_min, theta_max; NTHETA=',I2,/,2F8.3,
     &       ' = phi_min, phi_max   ;   NPHI =',I2,/,
     &       ' Results averaged over ',I3,' target orientations',/,
     &       '                   and ',I3,' incident polarizations')
 9050 FORMAT (10X,'Qext',6X,'Qabs',7X,'Qsca',5X,'g=<cos>',4X,'Qbk',7X,
     &       'Qpha',/,1X,'JO=1: ',1P,3E10.3,E11.3,2E10.3)
 9055 FORMAT (1X,'JO=2: ',1P,3E10.3,E11.3,2E10.3,/,1X,'mean: ',3E10.3,
     &       E11.3,2E10.3,/,1X,'Qpol= ',E10.3,35X,'dQpha=',E10.3)
 9056 FORMAT (1P,E10.4,E11.4,3E10.3,E11.3,E10.3)
 9057 FORMAT (1P,E10.4,E11.4,3E11.3)
 9060 FORMAT ('**** Selected scattering directions ',
     &       ' [note: incident pol state 1 is FIXED!]',/,
     &       ' ND THETA   PHI <|f11|^2> <|f21|^2> Re<f11*f21>',
     &       ' Im<f11*f21>')
 9070 FORMAT (I3,2F6.1,1P,2E10.3,2E11.3)
C 9080 FORMAT(' ***************** Selected scattering directions *****',
C     &       '************',/,
C     &       ' ND THETA   PHI <|f11|^2> <|f21|^2> <|f12|^2> ',
C     &       '<|f22|^2>  DEPOLR    PPOL')
 9080 FORMAT (12X,'** Mueller matrix elements for selected scattering ',
     &       'directions **',/,
     &       'theta   phi   S_11       S_21       S_31       S_41',
     &       '       S_12       S_13      Pol.')
 9090 FORMAT (F5.1,F6.1,1P,E10.3,1P,5E11.3,0P,F9.5)
 9150 FORMAT (8X,'Qsca*g(1)  Qsca*g(2)  Qsca*g(3)',/,1X,'JO=1:',1P,
     &       3E11.3)
 9155 FORMAT (1X,'JO=2:',1P,3E11.3,/,1X,'mean:',1P,3E11.3)
 9250 FORMAT (8X,'Qtrqab(1)  Qtrqab(2)  Qtrqab(3)  ',
     &       'Qtrqsc(1)  Qtrqsc(2)  Qtrqsc(3)',/,1X,'JO=1:',1P,6E11.3)
 9255 FORMAT (1X,'JO=2:',1P,6E11.3,/,1X,'mean:',1P,6E11.3)

      END

