      SUBROUTINE PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,DX,A1,A2)
      IMPLICIT NONE
C
C Arguments:
C
      INTEGER MXNAT,NAT
      INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3)
      REAL A1(3),A2(3),DX(3)
C
C Local variables:
C
      INTEGER LWORK
      PARAMETER(LWORK=12)
      INTEGER 
     &   INFO,
     &   J,J1,J2,J3,JA,JD,
     &   K,LDA,LDVL,LDVR,N
      REAL DENOM,PI,ZNORM,ZNORM2
      REAL 
     &   A(3,3),A3(3),VL(3,3),WI(3),WORK(LWORK),WR(3),
     &   ZR(3,3)
      DOUBLE PRECISION DSUM
      DOUBLE PRECISION CM(3)
      CHARACTER CMSGNM*70
C***********************************************************************
C Subroutine to determine principal axes of moment of inertia tensor
C of target.
C Target is assumed to be homogeneous.
C Given:
C
C MXNAT       =limit on largest allowed value of NAT=number of dipoles
C NAT         =number of dipoles in target
C IXYZ(J,1-3) =x,y,z location of dipole Jon lattice
C DX(1-3)     =(dx/d,dy/d,dz/d) where (dx,dy,dz)=lattice spacing in
C              x,y,z directions, 
C              and d = (dx*dy*dz)**(1/3) = effective lattice spacing
C ICOMP(J,1-3)=composition for dipole J; x,y,z directions
C
C Returns:
C
C A1(1-3)     =principal axis A1 in target frame (normalized)
C A2(1-3)     =principal axis A2 in target frame (normalized)
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, 96.01.26
C History:
C 96.01.26 (BTD) extracted from routine TARBLOCKS
C                IX,IY,IZ replaced by IXYZ
C 96.02.23 (BTD) corrected typo in 2 format statements
C                add code computing axis a3 and printing out
C                eigenvalue.
C 96.11.01 (PJF) modified to use LAPACK code rather than EIGV
C                and associated routines from NSWC package.
C 96.11.21 (BTD) modified to remove WRITE(0 and replace with
C                calls to WRIMSG
C 97.12.26 (BTD) added DX(3) to argument list, and modified to
C                compute principal axes with unequal lattice
C                spacing.
C 98.04.27 (BTD) deleted unused variables I,ZI
C 98.12.15 (BTD) Appears to be a problem with LAPACK routine SGEEV
C                when compiled using g77 under Linux.
C                Add trap to catch this problem if it occurs,
C                print out warning, and terminate execution.
C 00.06.13 (BTD) corrected error in computation of diagonal elements
C                of moment of inertia tensor (NAT/6 -> NAT/18)
C End history.
C
C Copyright (C) 1993,1994,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
C Compute moment of inertia tensor.
C
      DO JD=1,3
         CM(JD)=0.
      ENDDO
      DO JA=1,NAT
         DO JD=1,3
            CM(JD)=CM(JD)+IXYZ(JA,JD)*DX(JD)
         ENDDO
      ENDDO
      DO JD=1,3
         CM(JD)=CM(JD)/NAT
      ENDDO
C
C Compute moment of inertia tensor
C
      DO J=1,3
         DO K=1,3
            A(J,K)=0.
         ENDDO
      ENDDO
      DSUM=0.
      DO K=1,3
         DO JA=1,NAT
            DSUM=DSUM+((IXYZ(JA,K)-CM(K))*DX(K))**2
         ENDDO
      ENDDO
      DO J=1,3
         A(J,J)=REAL(DSUM+(NAT/18.)*(DX(1)**2+DX(2)**2+DX(3)**2))
      ENDDO
      DO J=1,3
         DO K=1,3
            DSUM=0.
            DO JA=1,NAT
               DSUM=DSUM+(IXYZ(JA,J)-CM(J))*(IXYZ(JA,K)-CM(K))*
     &              DX(J)*DX(K)
            ENDDO
            A(J,K)=A(J,K)-REAL(DSUM)
         ENDDO
      ENDDO
C
C Normalize moment of inertia tensor to units of 0.4*mass*a_eff**2
C
      DENOM=0.4*NAT*(.75*NAT/PI)**(2./3.)
      DO 3500 K=1,3
         DO 3400 J=1,3
            A(J,K)=A(J,K)/DENOM
 3400    CONTINUE
 3500 CONTINUE
C
C Compute eigenvalues and eigenvectors of A(J,K)
C
      N=3
      LDA=3
      LDVR=3
      LDVL=3

C Call LAPACK routine SGEEV to obtain eigenvectors of moment of
C inertia tensor.  Only right eigenvectors are needed.

      CALL SGEEV('V','V',N,A,LDA,WR,WI,VL,LDVL,ZR,
     &           LDVR,WORK,LWORK,INFO)

c  ZR      (output) REAL array, dimension (LDVR,N)
c          If JOBVR = 'V', the right eigenvectors v(j) are stored one
c          after another in the columns of ZR, in the same order
c          as their eigenvalues.
c          If JOBVR = 'N', ZR is not referenced.
c          If the j-th eigenvalue is real, then v(j) = ZR(:,j),
c          the j-th column of ZR.
c          If the j-th and (j+1)-st eigenvalues form a complex
c          conjugate pair, then v(j) = ZR(:,j) + i*ZR(:,j+1) and
c          v(j+1) = ZR(:,j) = i*ZR(:,j+1).

C 98.12.06 BTD It appears that SGEEV or subsidiary LAPACK routines
C              are not executing properly under Linux g77.  This 
C              problem has not been noticed previously, and certainly
C              did not appear with the Solaris f77 compiler.  Perhaps
C              one of the LAPACK routines is failing to SAVE some
C              variable between calls.  In any event, we have 
C              introduced this trap to terminate execution in case
C              of failure of SGEEV.
C
      IF(INFO.NE.0)THEN
         WRITE(CMSGNM,FMT='(A,I4)')
     &      'Error in PRINAXIS: returned from SGEEV with INFO=',INFO
         WRITE(CMSGNM,FMT='(A)')
     &      'SGEEV is part of lapack.f'
         WRITE(CMSGNM,FMT='(A)')
     &      'LAPACK routines may not work under Linux g77'
         WRITE(CMSGNM,FMT='(A)')
     &      'workaround by specifying target axis by hand'
         CALL ERRMSG('FATAL','PRINAXIS',' failure in SGEEV')
      ENDIF

C  
C Normalize the eigenvectors:
C
      DO 4000 J=1,3
         ZNORM2=0.
         DO 3800 K=1,3
            ZNORM2=ZNORM2+ZR(K,J)**2
 3800    CONTINUE
C
C Make first component of each eigenvector positive
C
         ZNORM=SQRT(ZNORM2)
         IF(ZR(1,J).LT.0.)THEN
            ZNORM=-ZNORM
         ELSEIF(ZR(1,J).EQ.0.)THEN
            IF(ZR(2,J).LT.0.)ZNORM=-ZNORM
         ENDIF
         DO 3900 K=1,3
            ZR(K,J)=ZR(K,J)/ZNORM
 3900    CONTINUE
 4000 CONTINUE
C
C Order the eigenvectors by eigenvalue, largest first:
C
      J1=1
      IF(WR(2).GT.WR(1))J1=2
      IF(WR(3).GT.WR(J1))J1=3
      IF(J1.EQ.1)THEN
         J2=2
         IF(WR(3).GT.WR(2))J2=3
      ELSEIF(J1.EQ.2)THEN
         J2=1
         IF(WR(3).GT.WR(1))J2=3
      ELSEIF(J1.EQ.3)THEN
         J2=1
         IF(WR(2).GT.WR(1))J2=2
      ENDIF
      DO 4100 K=1,3
         A1(K)=ZR(K,J1)
         A2(K)=ZR(K,J2)
 4100 CONTINUE
C
C Set J3=3 if J1+J2=1+2=3
C        2         =1+3=4
C        1         =2+3=5
C
      J3=6-J1-J2
C
C Having fixed eigenvectors a_1 and a_2, compute a_3=a_1xa_2
C
      A3(1)=A1(2)*A2(3)-A1(3)*A2(2)
      A3(2)=A1(3)*A2(1)-A1(1)*A2(3)
      A3(3)=A1(1)*A2(2)-A1(2)*A2(1)      
C
C*** enable following lines to print out eigenvalues
      WRITE(CMSGNM,FMT='(A,0PF8.5)')'I_1/.4Ma_eff^2 =',WR(J1)
      CALL WRIMSG('PRINAXIS',CMSGNM)
      WRITE(CMSGNM,FMT='(A,0PF8.5)')'I_2/.4Ma_eff^2 =',WR(J2)
      CALL WRIMSG('PRINAXIS',CMSGNM)
      WRITE(CMSGNM,FMT='(A,0PF8.5)')'I_3/.4Ma_eff^2 =',WR(J3)
      CALL WRIMSG('PRINAXIS',CMSGNM)
      WRITE(CMSGNM,FMT='(A,0P3F9.5,A)')'axis a_1 = (',A1,')'
      CALL WRIMSG('PRINAXIS',CMSGNM)
      WRITE(CMSGNM,FMT='(A,0P3F9.5,A)')'     a_2 = (',A2,')'
      CALL WRIMSG('PRINAXIS',CMSGNM)
      WRITE(CMSGNM,FMT='(A,0P3F9.5,A)')'     a_3 = (',A3,')'
      CALL WRIMSG('PRINAXIS',CMSGNM)
C***
      RETURN
      END



      SUBROUTINE ROT2(A,THETA,RM)
      IMPLICIT NONE
C***********************************************************************
C This subroutine finds the (3 by 3) rotation matrix for the
C given axis (vector) --- a and the rotation angle --- theta.
C On input:
C a(3)    --- axis of rotation  (vector, any length)
C theta   --- angle of rotation (in radians) around vector a
C On output:
C rm(3,3) --- rotation matrix such that new vector v2 is obtained
C             from the old vector v1 by the transformation:
C             v2 = rm \dot v1
C (If more rotations are needed: generate matrices rm with the
C help of this subroutine and (matrix) multiply them before
C performing v2 = rm \dot v1.)
C
C Reference:
C Leubner, C., 1977,  Coordinate-free  rotation operator,
C Am. J. Phys. 47, 727---729.
C
C History:
C Written by PJF, May 1989, for use in the rotation-averaged
C DDA calculations.
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 pre-calculate some stuff:
C     .. Scalar Arguments ..
      REAL THETA
C     ..
C     .. Array Arguments ..
      REAL A(3), RM(3,3)
C     ..
C     .. Local Scalars ..
      REAL ANORM, CT, ST
      INTEGER I, J
C     ..
C     .. Local Arrays ..
      REAL AHAT(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL XNORM3
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC COS, SIN
C     ..
      CT = COS(THETA)
      ST = SIN(THETA)
      CALL XNORM3(A,ANORM)
      AHAT(1) = A(1)/ANORM
      AHAT(2) = A(2)/ANORM
      AHAT(3) = A(3)/ANORM
C  cos(\theta) {\bf 1}  term:
      RM(1,1) = CT
      RM(2,2) = CT
      RM(3,3) = CT
C Skew-term  a \times {\bf 1}:
      RM(1,2) = -ST*AHAT(3)
      RM(1,3) = ST*AHAT(2)
      RM(2,1) = ST*AHAT(3)
      RM(2,3) = -ST*AHAT(1)
      RM(3,1) = -ST*AHAT(2)
      RM(3,2) = ST*AHAT(1)
C aa-dyadic (outer-product):
      DO 20 I = 1, 3
          DO 10 J = 1, 3
              RM(I,J) = RM(I,J) + (1.-CT)*AHAT(J)*AHAT(I)
   10     CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE XNORM3(A,XN)
C auxilary for "rotation" routines
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C     .. Scalar Arguments ..
      REAL XN
C     ..
C     .. Array Arguments ..
      REAL A(3)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC SQRT
C     ..
      XN = SQRT(A(1)**2+A(2)**2+A(3)**2)
      RETURN
      END
      SUBROUTINE MULT3(A,B,C)
C***********************************************************************
C Purpose:
C Given:  A(I,J)= 3x3 matrix
C         B(J,K)= 3x3 matrix
C Returns:C(I,K)= matrix product of A and B
C This is auxilary for "rotation" routines
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
***********************************************************************
C Arguments:
      REAL A(3,3),B(3,3),C(3,3)
C Local variables:
      REAL S
      INTEGER I,J,K
C
      DO 30 J = 1, 3
          DO 20 I = 1, 3
              S = 0.
              DO 10 K = 1, 3
                  S = S + A(I,K)*B(K,J)
   10         CONTINUE
              C(I,J) = S
   20     CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE PROD3(RM,VIN,VOUT)
C Given:  RM(I,J)=3x3 matrix
C         VIN(J)=3-vector
C Returns:VOUT(I)=RM*VIN
C 
C Auxilary for "rotation" routines
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments;
      REAL RM(3,3), VIN(3), VOUT(3)
C Local variables:
      INTEGER I,J
C
      DO 20 I = 1, 3
         VOUT(I) = 0.
         DO 10 J = 1, 3
            VOUT(I) = VOUT(I) + RM(I,J)*VIN(J)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE PROD3C(RM,CVIN,CVOUT)
C Arguments:
      REAL RM(3,3)
      COMPLEX CVIN(3),CVOUT(3)
C Local variables:
      INTEGER I,J
C Given:
C       RM=3x3 rotation matrix
C       CVIN=complex 3-vector
C Computes:
C       CVOUT=RM \dot CVIN
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      DO 20 I=1,3
         CVOUT(I)=(0.,0.)
         DO 10 J=1,3
            CVOUT(I)=CVOUT(I)+RM(I,J)*CVIN(J)
 10      CONTINUE
 20   CONTINUE
      RETURN
      END
      SUBROUTINE PROD3V(RM,VIN,VOUT,NV,NVMAX)
C Arguments:
      INTEGER NV,NVMAX
      REAL RM(3,3),VIN(3,NVMAX),VOUT(3,NVMAX)
C Local variables:
      INTEGER I,J,N
C Given:
C       RM = 3x3 rotation matrix
C       VIN = NV 3-vectors
C Computes:
C       VOUT=RM \dot VIN for each of the NV 3-vectors
C
C Copyright (C) 1993,1996 B.T. Draine and P.J. Flatau
C History
C 96.08.09 (BTD) trivial modification to remove obsolescent construct.
C end history
C This code is covered by the GNU General Public License.
C***********************************************************************
      DO 20 I=1,3
         DO 10 N=1,NV
            VOUT(I,N)=0.
 10      CONTINUE
 20   CONTINUE
C
      DO 50 I=1,3
         DO 40 J=1,3
            DO 30 N=1,NV
               VOUT(I,N)=VOUT(I,N)+RM(I,J)*VIN(J,N)
 30         CONTINUE
 40      CONTINUE
 50   CONTINUE
      RETURN
      END

      SUBROUTINE ROTATE(A1,A2,AK1,CXE01,CXE02,BETA,THETA,PHI,
     &                  EN0R,CXE01R,CXE02R,MXSCA,NSCAT,ENSC,EM1,EM2,
     &                  AKSR,EM1R,EM2R)
      IMPLICIT NONE

C Arguments:

      INTEGER MXSCA,NSCAT
      REAL AK1,BETA,THETA,PHI
      REAL A1(3),A2(3),AKSR(3,MXSCA),
     &     EM1(3,MXSCA),EM1R(3,MXSCA),EM2(3,MXSCA),
     &     EM2R(3,MXSCA),EN0R(3),ENSC(3,MXSCA)
      COMPLEX CXE01(3),CXE02(3),CXE01R(3),CXE02R(3)

C Local variables

      INTEGER I
      REAL BETA0,PHI0,PI,SINTHE,THETA0
      REAL RM1(3,3),RM2(3,3),RM3(3,3),VEC(3)
C***********************************************************************
C Given:
C       A1(1-3)=axis 1 of (original, unrotated) target
C       A2(1-3)=axis 2 of (original, unrotated) target
C       AK1=magnitude of k vector
C       CXE01(1-3)=original incident polarization vector 1
C       CXE02(1-3)=original incident polarization vector 2
C       BETA=angle (radians) through which target is to be rotated
C            around target axis A1
C       THETA=angle (radians)which target axis A1 is to make with
C            Lab x-axis
C       PHI=angle (radians) between Lab x,A1 plane and Lab x,y plane
C       ENSC(1-3,1-NSCAT)=scattering directions in Lab Frame
C       EM1(1-3,1-NSCAT)=scattering polarization vector 1 in Lab Frame
C       EM2(1-3,1-NSCAT)=scattering polarization vector 2 in Lab Frame
C
C Returns:
C       EN0R(1-3)=incident propagation vector in Target Frame
C       CXE01R(1-3)=incident polariz. vector 1 in Target Frame
C       CXE02R(1-3)=incident polariz. vector 2 in Target Frame
C       AKSR(1-3,1-NSCAT)=scattering k vectors in Target Frame
C       EM1R(1-3,1-NSCAT)=scattering pol. vector 1 in Target Frame
C       EM2R(1-3,1-NSCAT)=scattering pol. vector 2 in Target Frame
C
C Purpose:
C       In the Lab Frame, we hold the propagation direction (x-axis)
C       and incident polarization vectors (CXE01 and CXE02) fixed,
C       and rotate target to an orientation specified by the three
C       angles BETA,THETA,PHI.
C       However, computations in the main program DDSCAT are carried out
C       in the Target Frame, in which the target is held fixed and the
C       propagation vectors and polarization vectors are rotated.
C       Hence, this routine computes the incident propagation vector,
C       scattering propagation vectors, and associated polarization
C       vectors in the Target Frame.
C
C B. T. Draine, Princeton Univ. Obs., 89.11.20
C History:
C 91.04.22 (BTD) Eliminated vector A1R and removed superfluous line
C                computing A1R
C 96.11.03 (BTD) Corrected comments
C 03.06.06 (BTD) Corrected calculation of BETA0,PHI0 -- it appears that
C                these were not correct when a_1z < 0
C                                     and/or a_3x < 0
C end history
C
C Copyright (C) 1996,2003 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      PI=4.*ATAN(1.)

c Determine initial orientation THETA0,PHI0,BETA0 of target:
c theory:

c     theta = arccos(a_1x)

c if sin(theta) .ne. 0 , then

c     cos(phi)=a_1y/sin(theta)      sin(phi)=a_1z/sin(theta)
c     if(a_1z/sin(theta) > 0) then phi = arccos(a_1y/sin(theta))
c     if(a_1z/sin(theta) < 0) then phi = 2*pi-arccos(a_1y/sin(theta))

c     cos(beta)=-a_2x/sin(theta)    sin(beta)=a_3x/sin(theta)
c     if(a_3x/sin(theta) > 0) then beta = arccos(-a_2x/sin(theta))
c     if(a_3x/sin(theta) < 0) then beta = 2*pi-arccos(-a_2x/sin(theta))

c note: a_3x = a_1y*a_2z - a_1z*a_2y

c if sin(theta) = 0 , then

c     phi = 0

c     cos(beta)=a_2y                sin(beta)=a_2z
c     if(a_2z) > 0)      then      beta = arccos(a_2y)
c     if(a_2z) < 0)      then      beta = 2*pi-arccos(a_2y)

      THETA0=ACOS(A1(1))
      SINTHE=SIN(THETA0)
      IF(SINTHE.NE.0.)THEN
         IF(A1(3).GE.0.)THEN
            PHI0=ACOS(A1(2)/SINTHE)
         ELSE
            PHI0=2.*PI-ACOS(A1(2)/SINTHE)
         ENDIF
         IF(A1(2)*A2(3)-A1(3)*A2(2).GE.0.)THEN         
            BETA0=ACOS(-A2(1)/SINTHE)
         ELSE
            BETA0=2.*PI-ACOS(-A2(1)/SINTHE)
         ENDIF
      ELSE
         PHI0=0.
         IF(A2(3).GE.0.)THEN
            BETA0=ACOS(A2(2))
         ELSE
            BETA0=2.*PI-ACOS(A2(2))
         ENDIF
      ENDIF

C**** First rotate the target through angle PHI-PHI0 around x axis
C    (i.e., rotate the lab through angle PHI0-PHI around x axis)
C** obtain rotation matrix:

      VEC(1)=1.
      VEC(2)=0.
      VEC(3)=0.
      CALL ROT2(VEC,PHI0-PHI,RM1)

C**** Now rotate the target through angle THETA-THETA0 around axis
C     VEC = xaxis cross a1
C     (i.e., rotate the lab through angle THETA0-THETA around
C     xaxis cross a1 )
C** first compute rotation axis

      VEC(1)=0.
      VEC(2)=-A1(3)
      VEC(3)=A1(2)

C*** For special case where a1=xaxis, we want rotation axis=a1 cross a2

      IF((VEC(2)**2+VEC(3)**2).LT.1.E-10)THEN
         VEC(2)=A1(3)*A2(1)-A1(1)*A2(3)
         VEC(3)=A1(1)*A2(2)-A1(2)*A2(1)
      ENDIF
      CALL ROT2(VEC,THETA0-THETA,RM2)
      CALL MULT3(RM2,RM1,RM3)

C**** Now rotate the target through angle BETA-BETA0 around a1
C     (i.e., rotate the lab through angle BETA0-BETA around a1)

      CALL ROT2(A1,BETA0-BETA,RM1)
      CALL MULT3(RM1,RM3,RM2)

C**** RM2 is now the full rotation matrix
C
C**** Now rotate all required vectors:
C Incident propagation vector:

      VEC(1)=1.
      VEC(2)=0.
      VEC(3)=0.
      CALL PROD3(RM2,VEC,EN0R)

C Incident polarization vectors:

      CALL PROD3C(RM2,CXE01,CXE01R)
      CALL PROD3C(RM2,CXE02,CXE02R)

C Scattering vectors:

      CALL PROD3V(RM2,ENSC,AKSR,NSCAT,MXSCA)
      DO 50 I=1,NSCAT
         AKSR(1,I)=AK1*AKSR(1,I)
         AKSR(2,I)=AK1*AKSR(2,I)
         AKSR(3,I)=AK1*AKSR(3,I)
   50 CONTINUE

C Scattering polarization vectors:

      CALL PROD3V(RM2,EM1,EM1R,NSCAT,MXSCA)
      CALL PROD3V(RM2,EM2,EM2R,NSCAT,MXSCA)
      RETURN
      END



      SUBROUTINE SCAT(AK,AKS,DX,EM1,EM2,E02,CMDTRQ,
     &    CBKSCA,CSCA,CSCAG,CTRQIN,CTRQSC,
     &    CXE,CXE00,CXF1L,CXF2L,CXP,CXSCR1,CXSCR2,CXSCR3,CXSCR4,
     &    ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR,
     &    PHIN,SCRRS1,SCRRS2,THETAN,IXYZ)
      IMPLICIT NONE
C
C Arguments:
      CHARACTER CMDTRQ*6
      REAL CBKSCA,CSCA,E02
      INTEGER ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR
      INTEGER*2 IXYZ(MXNAT,3)
      COMPLEX CXE(NAT,3),CXE00(3),CXF1L(MXSCA),CXF2L(MXSCA),CXP(NAT,3),
     &        CXSCR1(MXN3),CXSCR2(MXN3),CXSCR3(MXN3),CXSCR4(MXNAT,3)
      REAL AK(3),AKS(3,MXSCA),CSCAG(3),CTRQIN(3),CTRQSC(3),DX(3),
     &     EM1(3,MXSCA),EM2(3,MXSCA),
     &     PHIN(MXSCA),SCRRS1(NAT,3),SCRRS2(MXNAT),THETAN(MXSCA)
C
C Local variables:
      COMPLEX CXI,CXSCL1,CXSCL2,CXSCL3
      REAL AK2,AK3,AKK,CONST,COSPHI,COSTH,DCOSTH,DPHI,EINC,ES2,
     &     PI,PHI,RRR,SINPHI,SINTH,TERM,W,WTH
      INTEGER ICOSTH,IPHI,IPHM,J,K,ND
      COMPLEX CXBS(3),CXES(3),CXTRM(3)
      REAL AKS0(3),AKS1(3),AKS2(3),AKSN(3),CTRQSCTF(3),XYZCM(3)
C
C Intrinsic functions:
      INTRINSIC COS,EXP,MIN,SIN,SQRT,CONJG,REAL
C
C Data statements:
      DATA CXI/(0.,1.)/
C***********************************************************************
C
C SCAT computes energy radiated by array of oscillating dipoles
C and corresponding scattering properties if oscillations are in
C response to incident E field.
C
C Given:
C     AK(1-3) = (kx,ky,kz)*d = where (kx,ky,kz)=incident k vector
C               d = (dx*dy*dz)**(1/3) = effective lattice spacing
C     AKS(3,MXSCA)=scattered k vectors
C     DX(1-3) = (dx,dy,dz)/d where dx,dy,dz=x,y,z lattice spacing
C     CMDTRQ = 'DOTORQ' to calculate torques
C            = 'NOTORQ' to skip calculation of torques
C     EM1(3,MXSCA)=pol.vector 1 for each scattering direction
C     EM2(3,MXSCA)=pol.vector 2 for each scattering direction
C     E02 = |E_0|^2 , where E_0 = incident complex E field
C     NAT = number of dipoles
C     NAT3 = 3*NAT
C     IXYZ(NAT,1-3) = x/d, y/d, z/d (integers) for dipoles 1-NAT
C     CXE(1-NAT,1-3) = components of E field at each dipole at t=0
C     CXE00(1-3) = (complex) reference polarization vector
C     CXP(1-NAT,1-3) = components of polarization of each dipole
C     ICTHM = number of THETA values for calculation of CSCA and G
C             (ICTHM must be odd; recommend ICTHM > 5*(2*pi*a/lambda)
C              for accuracy)
C     IPHIM = number of PHI values for calculation of CSCA and G
C     NDIR = number of directions at which scattering matrix elements
C            are to be computed
C     THETAN(1-NDIR) : scattering angle (radians);
C                      cos(THETAN) = n0 dot n
C     PHIN(1-NDIR)   : scattering angle (radians);
C                      PHIN is azimuthal angle of n from plane defined
C                      by n0 and Re(CXE00)
C
C Returns:
C     CBKSCA=differential scattering cross section for theta=pi
C            (lattice units)
C     CSCA = scattering cross section (lattice units)
C     CSCAG(1) = CSCA*<cos(theta)> , where theta=scattering angle
C     CSCAG(2) = CSCA*<sin(theta)cos(phi)>
C     CSCAG(3) = CSCA*<sin(theta)sin(phi)>
C     CTRQIN(1-3)=cross section for torque on grain due to incident
C                 E and B fields,
C                 for torque in directions x,y,z assuming incident 
C                 radiation is along x direction, 
C                 and "reference pol state" is in y direction.
C     CTRQSC(1-3)=cross section for torque on grain due to scattered
C                 E and B fields,
C                 for torque in directions x,y,z assuming incident 
C                 radiation is along x direction, 
C                 and "reference pol state" is in y direction
C                 Here "torque cross section" is defined so that
C                 torque(1-3)=(momentum flux)*(1/k)*C_torque(1-3)
C                 Total torque cross section=CTRQIN+CTRQSC
C     F1L(NDIR) = f_{1L} for direction NDIR
C     F2L(NDIR) = f_{2L} for direction NDIR
C                where f_{1L} connects incident polarization L to
C                outgoing theta polarization (E in scattering plane)
C                f_{2L} connects incident polarization L to outgoing phi
C                polarization (E perpendicular to scattering plane)
C                This is same f_{ml} notation as used by Draine (1988;
C                Ap.J., 333, 848).
C
C Notes:
C     AKS0(1-3)=k_s*d=scattered k vector (lattice units)
C     AKSN(1-3)=k_s/|k_s|=unit vector in scattering direction
C Scratch vectors introduced to permit vectorization:
C     SCRRS1 = real vector of length.GE.3*NAT
C     SCRRS2 = real vector of length.GE.NAT
C     CXSCR1 = complex scratch array of length 3*MXNAT
C     CXSCR2 = complex scratch array of length 3*MXNAT
C     CXSCR3 = complex scratch array of length 3*MXNAT
C     CXSCR4 = complex scratch array of length 3*MXNAT
C
C B.T.Draine, Princeton Univ. Obs., 87.01.04
C History:
C 88.06.29 (BTD): modified
C 89.11.28 (BTD): modified
C 90.11.02 (BTD): modified to enable use with vacuum sites.
C 90.11.20 (BTD): corrected error in calculation of f_ml (missing
C                 factor of k).
C 90.12.03 (BTD): changed ordering of XYZ0 and CXP
C                 eliminated cxe0 from argument list
C 90.12.10 (BTD): replaced XYZ0 with IXYZ0 to reduce memory use
C 91.08.14 (BTD): add CBKSCA to argument list for SCAT
C                 add code to eliminate sum over phi when theta=0 or pi
C                 add code to calculate CBKSCA
C 94.04.04 (BTD): removed superfluous CMPLX( ) operation from 2 lines
C                 in DO 200 loop. Also added new variable TERM to move
C                 some computation out of last DO loop.
C 95.06.14 (BTD): modified to compute moments <sin(theta)cos(phi)>
C                 and <sin(theta)sin(phi)> of scattered radiation
C                 added code to compute torque on grain relative to
C                 grain centroid XYZCM(1-3)
C                 NOTE: handling of grain centroid is not as efficient
C                 as it should be.  Can be recoded to reduce number
C                 of computations.
C 95.06.15 (BTD): modified to compute absorption contribution to
C                 torque on grain.
C                 added CXE and CTRQIN to argument list
C 95.06.15 (BTD): corrected errors found by Joseph Weingartner
C 95.06.16 (BTD): more corrections
C 95.06.20 (BTD): added argument CMDTRQ to allow torque calculations
C                 to be skipped if not desired.
C 95.06.22 (BTD): added factor CXI in terms appearing in evaluation
C                 of scattered torque
C 95.06.29 (BTD): Added loop to compute AKS0(1-3) for special cases
C                 theta-0 and 180. Error found by J. Weingartner.
C 95.07.20 (BTD): Add scratch vector SCRRS2
C                 Add scratch array CXSCR4 to argument list
C                 Rewrite to make torque computation more efficient:
C                 move as many calculations as possible out of loops 
C                 over scattering directions:
C                 SCRRS1 now used to store coords relative to centroid
C                 SCRRS2 now used to store (k_s dot r_j)
C                 CXSCR4 is needed to store (p_j cross r_j)
C 95.07.21 (BTD): Rewrite (again!) to make the correspondence
C                 to notation of Draine & Weingartner more apparent.
C                 This involved changing various definitions by
C                 factors of k, and change of sign of CXBS
C 95.07.24 (BTD): Corrected errors in evaluation of both scattering
C                 torque and incident torque
C 95.07.28 (BTD): Reenabled computation of CBKSCA
C 97.11.02 (BTD): Modified to allow use of noncubic lattice:
C                 added DX(1-3) to argument list and used DX in\
C                 evaluation of phase factors
C 98.08.10 (BTD): Added DX(1-3) in calculation of scattering amplitudes
C                 FML for selected scattering directions.
C
C Copyright (C) 1993,1994,1995,1997,1998 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      PI=4.*ATAN(1.)
C
C ICTHM,IPHIM determine angular grid used for calculation of total
C scattering cross section and G=<cos(theta)>
C ICTHM=number of theta values for which scattering is evaluated
C       (0 to pi)
C IPHIM=number of phi values for which scattering is evaluated
C       (0 to 2*pi)
C Note: ICTHM must be odd as Simpson's rule integration is employed.
C Check whether ICTHM is odd, and increase by 1 if not.
C
      IF(ICTHM.EQ.2*(ICTHM/2))ICTHM=ICTHM+1
C
C Incident k vector AK defines one axis (for spherical integration)
C Use (real part of) reference polarization vector CXE00 to define
C second coordinate axis.
C Construct third coordinate axis by taking cross product of first
C two axes.
C AKS1, AKS2 = vectors of length equal to AK lying along second and
C              third coordinate axes.
C
      AK2=AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3)
      AKK=SQRT(AK2)
      AK3=AKK*AK2
      AKS1(1)=REAL(CXE00(1))
      AKS1(2)=REAL(CXE00(2))
      AKS1(3)=REAL(CXE00(3))
C
C EINC = magnitude of (real part of) reference electric field
C
      EINC=SQRT(AKS1(1)**2+AKS1(2)**2+AKS1(3)**2)
      IF(EINC.LE.0.)CALL ERRMSG('FATAL','SCAT',
     &                          'CXE00 IS PURE IMAGINARY!')
      RRR=AKK/EINC
      AKS1(1)=RRR*AKS1(1)
      AKS1(2)=RRR*AKS1(2)
      AKS1(3)=RRR*AKS1(3)
      AKS2(1)=(AK(2)*AKS1(3)-AK(3)*AKS1(2))/AKK
      AKS2(2)=(AK(3)*AKS1(1)-AK(1)*AKS1(3))/AKK
      AKS2(3)=(AK(1)*AKS1(2)-AK(2)*AKS1(1))/AKK
C
      DCOSTH=2./REAL(ICTHM-1)
      DPHI=2.*PI/REAL(IPHIM)
      CSCA=0.
      DO K=1,3
         CSCAG(K)=0.
         CTRQSCTF(K)=0.
      ENDDO
C
C CONST = multiplicative factor entering weights
C factor 2/(ICTHM-1) = d(cos(theta))
C factor 2*pi/IPHIM = d(phi)
C factor 3 in denom from Simpson's rule (later mult. by 1,4,2..)
C Use Simpson's rule integration over d(cos(theta))
      CONST=4.*PI/ (3.*REAL((ICTHM-1)*IPHIM))
C Note that for cos(theta)=-1 or 1 we do not need to sum over phi,
C so treat these two directions separately
      IF(CMDTRQ.EQ.'DOTORQ')THEN
C
C***********************************************************************
C
C Additional quantities required for torque computations
C XYZCM(1-3) = coordinates of grain centroid
C SCRRS1(J,K)= r_j = coordinates of dipole J relative to centroid.
C CXSCR3(J) = r_j dot p_j
C CXSCR4(J,1-3) = p_j cross r_j
C
         DO K=1,3
            XYZCM(K)=0.
            DO J=1,NAT
               XYZCM(K)=XYZCM(K)+IXYZ(J,K)*DX(K)
            ENDDO
            XYZCM(K)=XYZCM(K)/NAT
         ENDDO
         DO K=1,3
            DO J=1,NAT
               SCRRS1(J,K)=IXYZ(J,K)*DX(K)-XYZCM(K)
            ENDDO
         ENDDO
         DO J=1,NAT
            CXSCR3(J)=(0.,0.)
            DO K=1,3
               CXSCR3(J)=CXSCR3(J)+SCRRS1(J,K)*CXP(J,K)
            ENDDO
         ENDDO
         DO J=1,NAT
            CXSCR4(J,1)=CXP(J,2)*SCRRS1(J,3)-CXP(J,3)*SCRRS1(J,2)
            CXSCR4(J,2)=CXP(J,3)*SCRRS1(J,1)-CXP(J,1)*SCRRS1(J,3)
            CXSCR4(J,3)=CXP(J,1)*SCRRS1(J,2)-CXP(J,2)*SCRRS1(J,1)
         ENDDO
C***********************************************************************
      ELSE
         DO K=1,3
            DO J=1,NAT
               SCRRS1(J,K)=IXYZ(J,K)*DX(K)
            ENDDO
         ENDDO
      ENDIF
C
C Proceed to sum over scattering angles:
C
      DO 80 ICOSTH=1,ICTHM
         IPHM=IPHIM
         WTH=4.
         IF(ICOSTH.GT.2*(ICOSTH/2))WTH=2.
C
C In special case ICOSTH=1 or ICTHM, we have theta=0 or pi,
C and the angle phi becomes irrelevant.  In this case we
C need only one value of phi, so we set IPHM=1 and adjust weights
C accordingly.
C
         IF(ICOSTH.EQ.1.OR.ICOSTH.EQ.ICTHM)THEN
            WTH=IPHIM
            IPHM=1
         ENDIF
         COSTH=1.-REAL(ICOSTH-1)*DCOSTH
         SINTH=SQRT(1.-MIN(COSTH*COSTH,1.))
C
C Evaluate total weight factor W
C
         W=CONST*WTH
         DO 70 IPHI=1,IPHM
            PHI=DPHI*(REAL(IPHI)-0.5)
            SINPHI=SIN(PHI)
            COSPHI=COS(PHI)
C
C Compute scattered k vector AKS0(1-3)
C Initialize CXES(1-3) = scattered electric field
C
            DO 10 K=1,3
               AKS0(K)=COSTH*AK(K)+SINTH*
     &                 (COSPHI*AKS1(K)+SINPHI*AKS2(K))
               CXES(K)=(0.,0.)
   10       CONTINUE
C
C Compute normalized scattering vector = nhat
C
            DO K=1,3
               AKSN(K)=AKS0(K)/AKK
            ENDDO
C
C Now compute CXES(1-3) = sum_j[p_j-nhat(nhat dot p_j)]exp[-ik_s dot x_j]
C      =(r/k*k)*exp(-ikr+iwt)*E(k_s,r,t)
C      where E(k_s,r,t) is complex electric field vector propagating
C      with wave vector k_s, at distance r from origin
C Notation:
C    If CMDTRQ='DOTORQ', then r_j = position relative to centroid XYZCM
C    If CMDTRQ='NOTORQ', then r_j = original lattice coords, and XYZCM=0.
C Note that since we are not concerned with overall phase shifts for
C evaluation of forces and torques, we do NOT need to worry about change 
C in phase resulting from using two different definitions of r_j when
C computing phase factors exp[-i k_s dot r_j]
C 
            DO 50 J=1,NAT
C
C CXSCR1(J) = nhat dot p(j)
C
               CXSCR1(J)=(0.,0.)
               DO 20 K=1,3
                  CXSCR1(J)=CXSCR1(J)+AKSN(K)*CXP(J,K)
   20          CONTINUE
C
C SCRRS2(J) = nhat dot r(j)
C
               SCRRS2(J)=0.
               DO 30 K=1,3
                  SCRRS2(J)=SCRRS2(J)+SCRRS1(J,K)*AKSN(K)
   30          CONTINUE
C
C CXSCR2(J) = exp(-i*k_s dot r_j) factor for atom j
C
               CXSCR2(J)=EXP(-CXI*AKK*SCRRS2(J))
C
C CXES(1-3) = sum_j[p_j-nhat(nhat dot p_j)]exp[-ik_s dot x_j]
C
               DO 40 K=1,3
                  CXES(K)=CXES(K)+(CXP(J,K)-AKSN(K)*CXSCR1(J))*
     &                    CXSCR2(J)
   40          CONTINUE
   50       CONTINUE
C
            IF(CMDTRQ.EQ.'DOTORQ')THEN
C***********************************************************************
C Additional terms required to compute torque on grain
C Have already computed
C CXES(1-3)=sum_j[p_j - k_s*(k_s dot p_j)/|k|^2]exp[-ik_s dot x_j]
C             where p_j = polarization of dipole j
C                   x_j = position of dipole j rel. to CM
C                   k_s = scattered k vector
C
C CXBS = -nhat cross CXES  = -(r/k*k)*exp(-ikr+iwt)*B(k_s,r,t)
C
               CXBS(1)=AKSN(3)*CXES(2)-AKSN(2)*CXES(3)
               CXBS(2)=AKSN(1)*CXES(3)-AKSN(3)*CXES(1)
               CXBS(3)=AKSN(2)*CXES(1)-AKSN(1)*CXES(2)
C
C CXSCL1 = sum_j p_j dot (r_j cross nhat)*exp(-i*k_s dot r_j)
C        = nhat dot sum_j (p_j cross r_j)*exp(-i*k_s dot r_j)
C recall that have previously evaluated
C CXSCR2(J) = exp(-i*k_s dot r_j) factor for atom j
C CXSCR4(J,K)=p_j cross r_j
C
               CXSCL1=(0.,0.)
               DO K=1,3
                  CXTRM(K)=(0.,0.)
               ENDDO
               DO K=1,3
                  DO J=1,NAT
                     CXTRM(K)=CXTRM(K)+CXSCR2(J)*CXSCR4(J,K)
                  ENDDO
                  CXSCL1=CXSCL1+AKSN(K)*CXTRM(K)
               ENDDO
C
C CXSCL2 = sum_j[(r_j dot p_j) - 
C                 (nhat dot p_j)*((nhat dot r_j) + 2i/k_s)]*
C                 exp(-i*k_s dot r_j)
C Have previously computed
C CXSCR1(J) = nhat dot p(j)
C CXSCR2(J) = exp(-i*k_s dot r_j) factor for atom j
C CXSCR3(J) = r_j dot p_j
C SCRRS2(J) = nhat dot r(j)
               CXSCL2=(0.,0.)
               DO J=1,NAT
                  CXSCL2=CXSCL2+CXSCR2(J)*(
     &                   CXSCR3(J)-CXSCR1(J)*(SCRRS2(J)+2.*CXI/AKK))
               ENDDO
C
C Note following notational correspondence:
C    present    Draine & Weingartner (1996)
C    -------    ---------------------------
C     CXES   =        V_E
C     CXBS   =        V_B
C    CXSCL1  =        S_B
C    CXSCL2  =        S_E
C
C***********************************************************************
            ENDIF
C
C Have completed computation of vector CXES(1-3) and terms to calculate
C angular momentum flux for scattering direction (THETA,PHI).
C CTRQSCTF(K) is "torque cross section" due to "scattered radiation"
C resulting in torque in lattice direction k (i.e., in "target
C frame".  After CTRQSCTF has been evaluated, we will compute CTRQSC in
C the "scattering frame", where the radiation is incident along the x direction,
C and "reference polarization state" is along the y direction.
C
            ES2=0.0
            DO 60 K=1,3
               ES2=ES2+REAL(CXES(K)*CONJG(CXES(K)))
   60       CONTINUE
            CSCA=CSCA+W*ES2
            CSCAG(1)=CSCAG(1)+W*COSTH*ES2
            CSCAG(2)=CSCAG(2)+W*SINTH*COSPHI*ES2
            CSCAG(3)=CSCAG(3)+W*SINTH*SINPHI*ES2
            IF(CMDTRQ.EQ.'DOTORQ')THEN
               DO K=1,3
                  CTRQSCTF(K)=CTRQSCTF(K)-
     &               W*REAL(CONJG(CXBS(K))*CXSCL2+CONJG(CXES(K))*CXSCL1)
               ENDDO
            ENDIF
   70    CONTINUE
C
C Save backscattering cross section:
C
         IF(ICOSTH.EQ.ICTHM)CBKSCA=AK2*AK2*ES2/E02
   80 CONTINUE
C
C***********************************************************************
C
      CSCA=AK2*AK2*CSCA/E02
C
C Convert CSCAG and CTRQSCTF to cross sections measured in lattice units:
C
      DO K=1,3
         CSCAG(K)=AK2*AK2*CSCAG(K)/E02
      ENDDO
      IF(CMDTRQ.EQ.'DOTORQ')THEN
C***********************************************************************
C Compute the two contributions to the torque on the target:
C (1) torque due to E and B fields of incident radiation (CTRQIN)
C (2) torque due to E and B fields of scattered radiation (CTRQSC)
C These are first computed in the target frame, then transformed to
C the lab frame (frame where incident radiation propagates along
C the x axis, and "reference pol. state 1" is along the y axis).
C
         DO K=1,3
            CTRQSCTF(K)=AK2*AK2*CTRQSCTF(K)/E02
         ENDDO
C
C CTRQSCTF(1-3) is the vector torque cross section/|k| measured in the
C target frame.  We now compute CTRQSC(1-3)=the vector torque cross
C section in the "scattering frame", where the radiation is incident along
C the x axis, and the "reference polarization state" (at t=0) is
C linearly polarized along the y axis.
C Note that the following computation introduces an additional factor
C of |k|.
         CTRQSC(1)=AK(1)*CTRQSCTF(1)+AK(2)*CTRQSCTF(2)+AK(3)*CTRQSCTF(3)
         CTRQSC(2)=AKS1(1)*CTRQSCTF(1)+AKS1(2)*CTRQSCTF(2)+
     &             AKS1(3)*CTRQSCTF(3)
         CTRQSC(3)=AKS2(1)*CTRQSCTF(1)+AKS2(2)*CTRQSCTF(2)+
     &             AKS2(3)*CTRQSCTF(3)
C
C Now compute vector torque cross section due to incident radiation:
C
         DO K=1,3
            CTRQIN(K)=0.
         ENDDO
C
C Redefine scratch vector CXSCR1:
C CXSCR1(J)=conjg(p_j) dot E_0j at t=0
C
         DO J=1,NAT
            CXSCR1(J)=0.
            DO K=1,3
               CXSCR1(J)=CXSCR1(J)+CONJG(CXP(J,K))*CXE(J,K)
            ENDDO
         ENDDO
         CXSCL1=0.
         CXSCL2=0.
         CXSCL3=0.
C
C Compute CXSCL1,CXSCL2,CXSCL3=components of the torque due to
C incident E and B fields in the target frame.
C 95.06.22 (BTD): added factor CXI in r cross k term appearing in
C                 evaluations of CXSCL1,CXSCL2,CXSCL2
C 95.07.21 (BTD): corrected sign mistake (+CXI -> -CXI)
C 95.07.24 (BTD): changed def. of CXSCR1 to agree with paper,
C                 changed calc. of CXSCLn to more closely mirror paper,
C                 and changed AKS0 to AK
         DO J=1,NAT
            CXSCL1=CXSCL1+CONJG(CXP(J,2))*CXE(J,3)-
     &                    CONJG(CXP(J,3))*CXE(J,2)-CXI*
     &             (AK(2)*SCRRS1(J,3)-AK(3)*SCRRS1(J,2))*CXSCR1(J)
            CXSCL2=CXSCL2+CONJG(CXP(J,3))*CXE(J,1)-
     &                    CONJG(CXP(J,1))*CXE(J,3)-CXI*
     &             (AK(3)*SCRRS1(J,1)-AK(1)*SCRRS1(J,3))*CXSCR1(J)
            CXSCL3=CXSCL3+CONJG(CXP(J,1))*CXE(J,2)-
     &                    CONJG(CXP(J,2))*CXE(J,1)-CXI*
     &             (AK(1)*SCRRS1(J,2)-AK(2)*SCRRS1(J,1))*CXSCR1(J)
         ENDDO
C
C Compute CTRQIN(1-3)=2*components of torque due to incident E and B
C fields (in lab frame), expressed as a cross section:
C
         CTRQIN(1)=AK(1)*REAL(CXSCL1)+AK(2)*REAL(CXSCL2)+
     &             AK(3)*REAL(CXSCL3)
         CTRQIN(2)=AKS1(1)*REAL(CXSCL1)+AKS1(2)*REAL(CXSCL2)+
     &             AKS1(3)*REAL(CXSCL3)
         CTRQIN(3)=AKS2(1)*REAL(CXSCL1)+AKS2(2)*REAL(CXSCL2)+
     &             AKS2(3)*REAL(CXSCL3)
         DO K=1,3
            CTRQIN(K)=4.*PI*CTRQIN(K)/E02
         ENDDO
      ENDIF
C***********************************************************************
      IF(NDIR.LE.0)RETURN
C
C*** Calculate the 2 matrix elements FM1 for selected directions
C F1L(NDIR) is to the polarization state e1 (E in scattering plane)
C F2L(NDIR) is to the polarization state e2 (E perp. to scatt. plane)
C
C IMPORTANT NOTE: If the incident radiation is not linearly polarized,
C with zero imaginary component to the polarization vector, then the
C F1L and F2L computed here are NOT the "usual" f_ml, which
C describe scattering from linearly polarized states l to linearly
C polarized states m.
C If one wishes to reconstruct the "usual" f_ml even when computing
C with elliptically polarized light, then one must add additional
C code to the routine SCAT to reconstruct the f_ml from the F1L and
C F2L computed here.
C
      TERM=AK3/SQRT(E02)
      DO 230 ND=1,NDIR
         CXF1L(ND)=(0.0,0.0)
         CXF2L(ND)=(0.0,0.0)
         DO 210 J=1,NAT
            CXSCR2(J)=(0.0,0.0)
            CXSCR3(J)=(0.0,0.0)
            DO 200 K=1,3
               CXSCR2(J)=CXSCR2(J)+CXP(J,K)*EM1(K,ND)
               CXSCR3(J)=CXSCR3(J)+CXP(J,K)*EM2(K,ND)
  200       CONTINUE
  210    CONTINUE
         DO 220 J=1,NAT
            CXSCR1(J)=EXP(-CXI*(AKS(1,ND)*IXYZ(J,1)*DX(1)+
     &                          AKS(2,ND)*IXYZ(J,2)*DX(2)+
     &                          AKS(3,ND)*IXYZ(J,3)*DX(3)))
            CXF1L(ND)=CXF1L(ND)+CXSCR2(J)*CXSCR1(J)
            CXF2L(ND)=CXF2L(ND)+CXSCR3(J)*CXSCR1(J)
  220    CONTINUE
C
C   Units: Dipole polarizations are in units of [E]*d**3, where d=lattice
C          spacing, and [E]=dimensions of E field
C          Therefore, up to this point CXF1L, CXF2L are in units
C          of [E]*d**3
C          Now multiply by AK3/SQRT(E02) (units of 1/([E]d**3)) to get
C          dimensionless result.
C
         CXF1L(ND)=CXF1L(ND)*TERM
         CXF2L(ND)=CXF2L(ND)*TERM
  230 CONTINUE
      RETURN
C
      END
      SUBROUTINE SCAVEC(MXSCA,NSCAT,THETAN,PHIN,CXE01,ENSC,EM1,EM2)
      IMPLICIT NONE
C***********************************************************************
C Given:
C      MXSCA=dimension information for ENSC,PHIN,THETAN
C      NSCAT=number of scattering directions
C      THETAN(1-NSCAT)=scattering angles theta
C      PHIN(1-NSCAT)=scattering angles phi
C      CXE01(1-3)=complex polarization vector 1 (phi=0 direction
C                 is defined by x,y plane (Lab Frame), where incident
C                 radiation propagates along the x axis.
C Returns:
C      ENSC(1-3,1-NSCAT)=scattering vectors in Lab Frame
C      EM1(1-3,1-NSCAT)=scattered pol vectors parallel to scat. plane
C                       in Lab Frame
C      EM2(1-3,1-NSCAT)=scattered pol vectors perp. to scat. plane
C                       in Lab Frame
C It is assumed that incident propagation vector is in x-direction
C in Lab Frame
C
C History:
C 96.11.06 (BTD): Changed definition of scattering angle phi
C                 Previously, phi was measured from plane containing
C                 incident k vector (i.e., Lab x-axis) and Re(CXE01)
C                 Henceforth, phi is measured from Lab x,y plane.
C
C Copyright (C) 1993,1996, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C
C Arguments:
      INTEGER MXSCA,NSCAT
      COMPLEX CXE01(3)
      REAL EM1(3,MXSCA), EM2(3,MXSCA), ENSC(3,MXSCA), PHIN(MXSCA),
     &     THETAN(MXSCA)
C
C Local variables:
      INTEGER I
      REAL COSPHI,COSTHE,SINPHI,SINTHE
C
C***********************************************************************
      DO 1000 I=1,NSCAT
         COSPHI=COS(PHIN(I))
         SINPHI=SIN(PHIN(I))
         COSTHE=COS(THETAN(I))
         SINTHE=SIN(THETAN(I))
         ENSC(1,I)=COSTHE
         ENSC(2,I)=SINTHE*COSPHI
         ENSC(3,I)=SINTHE*SINPHI
         EM1(1,I)=-SINTHE
         EM1(2,I)=COSTHE*COSPHI
         EM1(3,I)=COSTHE*SINPHI
         EM2(1,I)=0.
         EM2(2,I)=-SINPHI
         EM2(3,I)=COSPHI
 1000 CONTINUE
      RETURN
      END
