      SUBROUTINE ALPHA(AKR,CALPHA,CXALDS,CXALPH,CXALOF,CXALOS,
     &     CXE0R,CXEPS,CXSC,CXSCR1,CXZC,CXZW,DX,IBETH,IBETH1,
     &     ICOMP,IOCC,IPHI,IPHI1,IXYZ0,JORTH,MYID,MXCOMP,MXNAT,MXN3,
     &     NAT,NAT0,NAT3,NCOMP,NX,NY,NZ,CXRLOC,CSHAPE,SHPAR)
      IMPLICIT NONE
C*** Arguments:
      CHARACTER*6 CALPHA,CSHAPE
      INTEGER IBETH,IBETH1,IPHI,IPHI1,JORTH,MXCOMP,MXNAT,MXN3,MXMEM,
     &     MYID,NAT,NAT0,NAT3,NCOMP,NX,NY,NZ
      COMPLEX CXALOF(NAT,3),CXALPH(NAT,3),CXE0R(3),CXEPS(MXCOMP),
     &     CXALDS(NAT,3),CXALOS(NAT,3),CXSC(NAT,3,3),CXSCR1(NAT,3),
     &     CXZC(NX+1,NY+1,NZ+1,6),CXZW(2*NX,2*NY,2*NZ,3),
     &     CXRLOC(MXCOMP+1,3,3)
      INTEGER*2 ICOMP(NAT,3),IOCC(NAT),IXYZ0(MXNAT,3)
      REAL AKR(3),DX(3),SHPAR(6)
C
C*** Local Variables:
      LOGICAL INIT
      CHARACTER CMSGNM*70
      INTEGER IA,IC,L,M,N,P,INFO,MINVIWRK(3)
      REAL A1(3),A2(3),A3(3),AK1,AK2,B1,B2,EMKD,PI,SUM,
     &     SHAPEL(3,3),SHAPELT(3,3),WREV(3)
      COMPLEX CXI,CXRR,CXTERM,CX33(3,3),
     &     CXALPH2(3,3),CXTMP(3,3),CXTMP2(3,3),
     &     MINVZWRK(3)
      DOUBLE PRECISION ASHAP(3),LSHAP(3)
C
      SAVE INIT,CXI,PI
      DATA CXI/(0.,1.)/
      DATA INIT/.TRUE./
C*** Common:
      CHARACTER*6 CMETHD
      INTEGER DUM1,DUM2,DUM3,DUM4,DUM5,DUM6,DUM7,DUM8,DUM9,
     &     DUM10,DUM11,DUM12
C-----------------------------------------------------------------------
      COMMON /M6/  DUM1, DUM2, DUM3, DUM4, DUM5, DUM6, DUM7, DUM8, 
     &     DUM9, DUM10, DUM11, DUM12, CMETHD
C-----------------------------------------------------------------------
C***********************************************************************
C Given:
C       AKR(1-3)=(kx,ky,kz)*d, where d=effective lattice spacing
C       CALPHA = 'LATTDR' or 'SCLDR' = polarizability prescription
C       CMETHD = determines 3-d FFT routine used by ESELF
C       CSHAPE = descriptor of target shape
C       CXE0R(1-3) = polarization vector in lattice coordinates
C                    (assumed to be normalized)
C       CXEPS(1-NCOMP)=distinct values of dielectric constant
C       CXSC(1-NAT,1-3,1-3) = complex scratch space for SCLDR calculation
C       CXSCR1(1-NAT,1-3)   = complex scratch space for SCLDR calculation
C       CXZC = complex scratch space needed by ESELF
C       CXZW = complex scratch space needed by ESELF
C       DX(1-3)=(dx/d,dy/d,dz/d), where dx,dy,dz=lattice spacings in
C               x,y,z directions, and d=(dx*dy*dz)**(1/3)
C       IBETH= MYID+IBETH1 if first time through combined BETA/THETA 
C              orientation loop
C       IBETH1=starting value of IBETH (see above)
C       ICOMP(1-NAT3)=composition identifier for each lattice site and
C                     direction (ICOMP=0 if lattice site is unoccupied)
C                storage scheme: 
C                1x,2x,...,NATx,1y,2y,...,NATy,1z,2z,...,NATz
C       IOCC(1-NAT) == 0 if site unoccupied
C                   == 1 if site occupied
C       IPHI = IPHI1 if first time through PHI orientation loop
C       IPHI1= starting value of IPHI (see above)
C       JORTH= 1 if first incident polarization state
C       MXCOMP = dimensioning information
C       MXN3 = dimensioning information
C       MYID = parallel process identifier (=0 if only 1 process)
C       NAT  = number of sites in extended target (incl. vacuum sites)
C       NAT3 = 3*number of sites in extended target (incl. vacuum sites)
C       NCOMP= number of different dielectric tensor elements in target
C       NX   = x-dimension of extended target
C       NY   = y-dimension of extended target
C       NZ   = z-dimension of extended target
C       SHPAR(1-6)=target shape parameters
C
C Returns:
C       CXALPH(J,1-3)=(alpha_11,alpha_22,alpha_33)/d^3 for dipole J=1-NAT
C                     where alpha_ij=complex polarizability tensor.
C       CXALOF(J,1-3)=(alpha_23,alpha_31,alpha_12)/d^3 for dipole J=1-NAT
C
C***
C If CALPHA = LATTDR:
C    Compute dipole polarizability using "Lattice Dispersion Relation"
C    of Draine & Goodman (1993,ApJ,March 10).  It is required that
C    polarizability be such that an infinite lattice of such dipoles
C    reproduce the continuum dispersion relation for radiation
C    propagating with direction and polarization of radiation incident
C    on the DDA target.  This is the recommended option.
C
C***
C    Note: CXALPH = polarizability/d^3
C    In the event that ICOMP=0, then we set CXALPH=1.
C***********************************************************************
C B.T.Draine, Princeton Univ. Obs.
C History:
C 90.09.13 (BTD): Corrected error in DO loop limit.
C 90.11.01 (BTD): Special treatment for "vacuum" sites in order
C                 to allow FFT treatment.
C 90.11.02 (BTD): Set CXALPH=(1.,0.) at vacuum sites.
C 91.04.30 (BTD): Modified to include both O[(kd)^2] correction term
C                 from Goedecke & Obrien (1988) in addition to
C                 radiative reaction correction.
C 91.05.07 (BTD): Modified to allow easy choice among three methods
C                 for computing polarizability.
C 91.05.07 (BTD): Added call to WRIMSG to record which method in use
C 91.05.08 (BTD): Added CALPHA to argument list.
C 91.05.09 (BTD): Added printing of kd and magnitude of correction
C                 terms.
C 91.05.13 (BTD): Experiment with use of approximate directional
C                 average for coefficient of CXEPS in Lattice
C                 Disperision Relation theory
C 91.05.13 (BTD): Corrected error in radiative reaction correction
C                 (error presumably introduced during last week)
C 91.05.14 (BTD): Added separate options LDRXYZ and LDRAVG
C 91.05.15 (BTD): Added option LDR000 to omit CXEPS-dependent
C                 correction term
C 91.05.30 (BTD): Changed LDRXYZ -> LDR100
C                 Added option LDR111
C 91.09.12 (BTD): Corrected error in numerical coefficient for
C                 LDRAVG case
C 91.09.17 (BTD): Eliminate options LDR000,LDRXYZ,LDRAVG
C                 Introduce option LATTDR (LATTice Dispersion Relation)
C                 which includes explicit dependence
C                 on direction and polarization state
C                 remove AK1 from argument list
C                 add AKR, CXE0R to argument list
C 92.05.14 (BTD): Introduce option LDRISO to experiment with
C                 using average directional correction term
C                 rather than using correction for incident direction
C                 and polarization
C 92.05.16 (BTD): Correct error in LDRISO option: average SUM should be
C                 3/15 rather than 4/15
C 97.11.02 (BTD): Add DX to argument list to allow use with noncubic
C                 lattices (additional modifications still required!!)
C 97.12.26 (BTD): Add CXALOF to argument list to allow use with noncubic
C                 lattices (additional modifications still required).
C 97.12.30 (BTD): Add code to evaluate alpha for noncubic case,
C                 using subroutine NONCUBIC to evaluate the lattice
C                 sums R0,R1,R2,R3
C 98.01.14 (BTD): Change sign of radiative reaction correction term.
C 98.01.19 (BTD): Change treatment of R_1,C3, and C4
C                 and output value of C4
C                 temporary modification to allow value of C4 to be
C                 read in from file 'alpha.par'
C 98.01.22 (BTD): Minor change in computation of radiative reaction
C                 correction.
C 98.03.10 (BTD): set C4=-C3, disable reading C4 from 'alpha.par'
C 98.04.27 (BTD): declared IC1,IC2,IC3 as integers.
C 99.02.09 (BTD): experiment with change in expression used to calculate
C                 off-diagonal elements of polarizability tensor
C                 experiment with C4=-C3-R1 to "cancel" the contribution
C                 of R1 to the off-diagonal terms
C 03.01.29 (BTD): remove code to calculate alpha for noncubic lattice
C                 (hold for release in future version)
C end history
C Copyright (C) 1993,1996,1997,1998,1999 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      PI=4.*ATAN(1.)
      AK2=AKR(1)*AKR(1)+AKR(2)*AKR(2)+AKR(3)*AKR(3)
      AK1=SQRT(AK2)
      CXRR=-(AK1*AK2/1.5)*CXI
C*** EMKD = |m|*k_0*d , where |m|=refractive index
C                             k_0 = wave vector in vacuo
C                             d = lattice spacing
      EMKD=SQRT(CABS(CXEPS(1)))*AK1
C***********************************************************************
C*** Lattice dispersion relation
C
      IF(CALPHA.EQ.'LATTDR')THEN
         WRITE(CMSGNM,FMT='(A,F6.4)')
     &         'Lattice dispersion relation for |m|k_0d=',EMKD
         CALL WRIMSG('ALPHA ',CMSGNM)
C*** Compute sum (a_j*e_j)^2 , where a_j=unit propagation vector
C                                    e_j=unit polarization vector
         SUM=0.0
         DO L=1,3
            SUM=SUM+(AKR(L)*CABS(CXE0R(L)))**2
         ENDDO
         SUM=SUM/AK2
         B1=-1.8915316*AK2
         B2=(0.1648469-1.7700004*SUM)*AK2
         DO L=1,3
            DO IA=1,NAT
               IC=ICOMP(IA,L)
               IF(IC.GT.0)THEN

C*** First compute Clausius-Mossotti polarizability:

                  CXTERM=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.)

C*** Determine polarizability by requiring that infinite lattice of
C    dipoles have dipersion relation of continuum.

                  CXTERM=CXTERM/(1.+CXTERM*(B1+CXEPS(IC)*B2))

C*** Radiative-reaction correction:

                  CXALPH(IA,L)=CXTERM/(1.+CXTERM*CXRR)

c set off-diagonal terms to zero

                  CXALOF(IA,L)=0.

               ELSEIF(IC.EQ.0)THEN

C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites.

                  CXALPH(IA,L)=1.
                  CXALOF(IA,L)=0.
               ENDIF
            ENDDO
         ENDDO
         RETURN
      ELSE
         WRITE(0,*)'Error: invalid option for subroutine ALPHA'
         STOP
      ENDIF
      END
      SUBROUTINE PETR(X,B,WRK,LDA,IPAR,SPAR,MATVEC,CMATVEC)
      IMPLICIT NONE
C solves Ax=b
C
C Input:
C B(LOCLEN) - RHS (doesn't change on output)
C WRK(LDA, LOCLEN) - scratch array 
C IPAR - defines integer paramteres (see documentation of PIM)
C SPAR - define real parameters
C LDA - leading dimension of arrays 
C MATVEC -- name of procedure for computing y=Ax
C CMATVEC -- name of procedure for computing y= conj(A') x
C
C Output:
C X(LOCLEN)
C 
C Reference:
C    Petravic,M and Kuo-Petravic,G. 1979, JCP, 32,263
C History: 
C    Late 80'  (Coded by Draine)
C    Changes to include "cprod" (PJF+BTD)
C    95.06.01 (PJF) re-written to conform to PIM
C    95.08.11 (BTD) minor editing (comments, etc.)
C
C Parameters:
      INTEGER IACE,IAXI,IGI,IPI,IQI,IR
      PARAMETER (IACE=1,IGI=2,IPI=3,IQI=4,IAXI=5,IR=6)
C Arguments:
      INTEGER LDA
      INTEGER IPAR(*)
      REAL SPAR(*)
      COMPLEX B(*), WRK(LDA,*), X(*)
      EXTERNAL CMATVEC,MATVEC
C
C Local variables
      COMPLEX CDUMMY
      INTEGER IDUMMY,ITNO,L,LOCLEN,MAXIT,STATUS
      REAL
     &   ALPHAI,BETAI1,BNORM,DUMMY,EPSILON,EXITNORM,
     &   GIGI,GI1GI1,QIQI,RHSSTOP
      REAL SCSETRHSSTOP
      EXTERNAL PSCNRM2
C
C***********************************************************************
      LOCLEN = IPAR(4)
      MAXIT = IPAR(10)
      EPSILON=SPAR(1)
C
C Compute RHSSTOP = EPSILON*|B|
C
      RHSSTOP = SCSETRHSSTOP(B,CDUMMY,EPSILON,IPAR,CDUMMY,PSCNRM2)
C
C Compute conjg(A')*B
C
      CALL CMATVEC(B,WRK(1,IACE),IDUMMY)
      BNORM = 0.
      DO 170 L = 1,LOCLEN
         BNORM = BNORM + REAL(B(L)*CONJG(B(L)))
         WRK(L,IGI) = WRK(L,IACE)
         WRK(L,IPI) = WRK(L,IGI)
  170 CONTINUE
C
C Compute |QI>=A|PI>
C
      CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY)
      QIQI = 0.
      GIGI = 0.
      DO 190 L = 1,LOCLEN
         GIGI = GIGI + REAL(CONJG(WRK(L,IGI))*WRK(l,IGI))
         QIQI = QIQI + REAL(CONJG(WRK(L,IQI))*WRK(L,IQI))
  190 CONTINUE
      ALPHAI = GIGI/QIQI
      DO 200 L = 1,LOCLEN
         X(L) = X(L) + ALPHAI*WRK(L,IPI)
  200 CONTINUE
C
C compute |AX1>:
C
      CALL MATVEC(X,WRK(1,IAXI),IDUMMY)
C
      DO 340 ITNO = 2,MAXIT
         IPAR(11) = ITNO
C
C Transfer <GI|GI> -> <GI-1|GI-1>:
C
         GI1GI1 = GIGI
C
C compute |GI>=AC|E>-AC|AXI>:
C
         CALL CMATVEC(WRK(1,IAXI), WRK(1,IGI), IDUMMY)
C
C Compute GIGI=<GI|GI>:
         GIGI = 0.
         DO 230 L = 1,LOCLEN
            WRK(L,IGI) = WRK(L,IACE) - WRK(L,IGI)
  230    CONTINUE
         DO 235 L = 1,LOCLEN
            GIGI = GIGI + REAL(CONJG(WRK(L,IGI))*WRK(L,IGI))
  235    CONTINUE
C
C Compute BETAI-1=<GI|GI>/<GI-1|GI-1>:
C
         BETAI1 = GIGI/GI1GI1
C
C Compute |PI>=|GI>+BETAI-1|PI-1>:
C
         DO 240 L = 1,LOCLEN
            WRK(L,IPI) = WRK(L,IGI) + BETAI1*WRK(L,IPI)
  240    CONTINUE
C
C Compute |QI>=A|PI>:
C
         CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY)
C
C Compute <QI|QI>:
C
         QIQI = 0.
         DO 250 L = 1,LOCLEN
            QIQI = QIQI + REAL(CONJG(WRK(L,IQI))*WRK(L,IQI))
  250    CONTINUE
C
C Compute ALPHAI = <GI|GI>/<QI|QI>:
C
         ALPHAI = GIGI/QIQI
C
C Compute |XI+1>=|XI>+ALPHAI*|PI>:
C
         DO 260 L = 1,LOCLEN
            X(L) = X(L) + ALPHAI*WRK(L,IPI)
  260    CONTINUE
C
C Except every 10TH iteration, compute |AXI+1>=|AXI>+ALPHAI*|QI>
C Draine: warning this part perhaps not checked 
C
         IF (ITNO.NE.10* (ITNO/10)) THEN
         DO 270 L = 1,LOCLEN
            WRK(L,IAXI) = WRK(L,IAXI) + ALPHAI*WRK(L,IQI)
  270    CONTINUE
         ELSE 
            CALL MATVEC(X,WRK(1,IAXI),IDUMMY)
         ENDIF
C
C Compute residual vector |RI>=A|XI>-|B>
C
         DO 300 L=1, LOCLEN
            WRK(L,IR) = WRK(L,IAXI)-B(L)
 300     CONTINUE
C
C Call STOPCRIT to check whether to terminate iteration.
C Criterion: terminate if <RI|RI> < TOL**2 * <B|B>
C Return with STATUS=0 if this condition is satisfied.
C
         CALL STOPCRIT(B,WRK(1,IR),CDUMMY,CDUMMY,CDUMMY,WRK,
     +                 RHSSTOP,IDUMMY,EXITNORM,
     +                 STATUS,IPAR,DUMMY,DUMMY,DUMMY,DUMMY,PSCNRM2)
C
C Call PROGRESS to report "error" = sqrt(<RI|RI>)/sqrt(<B|B>)
C
         CALL PROGRESS(LOCLEN,ITNO,EXITNORM,CDUMMY,CDUMMY,CDUMMY)
C      
      IF(STATUS.EQ.0)GOTO 500
 340  CONTINUE
 500  CONTINUE
      IPAR(12) = 0
      RETURN
      END



      SUBROUTINE COPYIT(CXA,CXB,NPTS)
      IMPLICIT NONE
C*** Arguments:
      INTEGER NPTS
      COMPLEX CXA(*), CXB(*)
C
C*** Local variables:
      INTEGER IPTS
C***********************************************************************
C Purpose: To copy a complex vector from one memory location to another.
C***********************************************************************
      DO 10 IPTS=1,NPTS
         CXB(IPTS)=CXA(IPTS)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE DIAGL(U,V,IPAR)
      IMPLICIT NONE

C Arguments:
      COMPLEX U(*), V(*)
      INTEGER IPAR(*)

C Common variables:

      CHARACTER*6 CMETHD
      INTEGER
     &   IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &   NAT,NAT0,NAT3,NX,NY,NZ

      COMPLEX CXADIA(1)
C-----------------------------------------------------------------------
      COMMON/M2/CXADIA
      COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3,NAT0,NX,NY,NZ,
     &          MXN3F,IDVOUT,CMETHD
C-----------------------------------------------------------------------

C cg package
C Left preconditioning subroutine - division by diagonal
C History:
C Originally written by P.J. Flatau, 1993
C 00.06.26 (BTD) Removed MXMEMF from COMMON/M6/
C end history
c-----------------------------------------------------------------------

      CALL CDIV(U,V,CXADIA,NAT3)
      RETURN
      END

      SUBROUTINE CDIV(U,V,CXA,N)
      IMPLICIT NONE

C Arguments:

      INTEGER N
      COMPLEX U(*), V(*), CXA(*)
C
C Local variables:
      INTEGER I

c-----------------------------------------------------------------------
c Part of conjugate gradient package
c needed by left preconditioning subroutine
c Copyright (C) 2003 B.T. Draine and P.J. Flatau
c This code is covered by the GNU General Public License
c-----------------------------------------------------------------------
      DO I=1,N
         V(I)=U(I)/CXA(I)
      ENDDO
      RETURN
      END

      SUBROUTINE MATVEC(CXX,CXY,IPAR)
      IMPLICIT NONE
C Arguments:
      COMPLEX CXX(*), CXY(*), IPAR(*)
C Common:
      INTEGER*2 IOCC(1)
      INTEGER
     &   IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &   NAT,NAT0,NAT3,NX,NY,NZ
      CHARACTER*6 CMETHD
      REAL AK(3),DX(3)
      COMPLEX CXADIA(1),CXAOFF(1),CXZC(1),CXZW(1)
C-----------------------------------------------------------------------
      COMMON/M1/AK,DX
      COMMON/M2/CXADIA
      COMMON/M3/CXZC
      COMMON/M4/CXZW
      COMMON/M5/IOCC
      COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3,NAT0, 
     &  NX,NY,NZ,MXN3F,IDVOUT,CMETHD
      COMMON/M7/CXAOFF
C-----------------------------------------------------------------------
C Subroutine MATVEC

C matvec  --> cxy = A cxx          'N'
C tmatvec --> cxy = A' cxx         'T'
C cmatvec --> cxy = conjg(A') cxx  'C'
C
C 
C Notice that DDSCAT has cases 'N' and 'C' implemented in CPROD
C Notice, however, that in DDSCAT 'T'='N'.
C Some CG methods use cases N, C (like Petravic).
C Some CG methods use cases N, T (like PIM)
C Some CG use all three products N, C, T  (parts of CCGPAK)
C
C Information about matrix A is transfered via commons /M1/-/M7/ with
C the main program. Arrays cxadia, cxzc, cxzwm,  iocc are properly
C dimensioned inside "cprod". You may get warning during the compilation
C about "misaligned common". Ignore these. Also notice mxnatf,
C mxn3f,mxnxf, mxnyf, mxnzf. These are defined by "parameter
C statement" in main code and this is the reason why they have "f" suffix.
C
C ipar - information array, not used, defined for compatibility
C
C History:
C 97.11.01 (BTD): modified to allow for noncubic lattice, introduced
C                 new vector DX(1-3) via COMMON/M1/
C                 added DX to argument list of CPROD
C 97.12.25 (BTD): added COMMON/M7/CXAOFF to communicate off-diagonal
C                 elements of inverse polarizability tensors.
C                 Added CXAOFF to argument list of CPROD
C 00.06.13 (BTD): Explicit declaration of all variables.
C 00.06.25 (BTD): Removed MXMEMF from COMMON/M6/
C end history
C Copyright (C) 1993,1997,2000
C                B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************

      CALL CPROD(AK,DX,CMETHD,'N',CXADIA,CXAOFF,
     &           CXX,CXY,CXZC,CXZW,
     &           IDVOUT,IOCC,MXN3f,MXNATf,MXNXf,MXNYf,MXNZf,
     &           NAT,NAT0,NAT3,NX,NY,NZ)
      RETURN
      END

      SUBROUTINE CMATVEC(CXX, CXY, IPAR)
      IMPLICIT NONE

C Arguments:

      COMPLEX CXX(*), CXY(*), IPAR(*)

C Common:

      INTEGER*2 IOCC(1)
      INTEGER
     &   IDVOUT,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &   NAT,NAT0,NAT3,NX,NY,NZ
      CHARACTER*6 CMETHD
      REAL AK(3),DX(3)
      COMPLEX CXADIA(1),CXAOFF(1),CXZC(1), CXZW(1)
C-----------------------------------------------------------------------
      COMMON/M1/AK,DX
      COMMON/M2/CXADIA
      COMMON/M3/CXZC
      COMMON/M4/CXZW
      COMMON/M5/IOCC
      COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3,NAT0,NX,NY,NZ,
     &          MXN3F,IDVOUT,CMETHD
      COMMON/M7/CXAOFF
C-----------------------------------------------------------------------
c Subroutine CMATVEC
c
C matvec  --> cxy = A cxx   
C tmatvec --> cxy = A' cxx   
C cmatvec --> cxy = conjg(A') cxx   
C History:
C 97.11.01 (BTD): modified to allow for noncubic lattice, introduced
C                 new vector DX(1-3) via COMMON/M1/
C                 added DX to argument list of CPROD
C 97.12.25 (BTD): added COMMON/M7/CXAOFF to communicate off-diagonal
C                 elements of inverse polarizability tensors.
C                 Added CXAOFF to argument list of CPROD
C 00.06.13 (BTD): Explicit declaration of all variables.
C 00.06.26 (BTD): Removed MXMEMF from COMMON/M6/
C end history
C Copyright (C) 1993,1997,2000
C               B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************

      CALL CPROD(AK,DX,CMETHD,'C',CXADIA,CXAOFF,
     &           CXX,CXY,CXZC,CXZW,
     &           IDVOUT,IOCC,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF,
     &           NAT,NAT0,NAT3,NX,NY,NZ)
      RETURN
      END

      SUBROUTINE CPROD(AK,DX,CMETHD,CWHAT,CXADIA,CXAOFF,
     &                 CXX,CXY,CXZC,CXZW,
     &                 IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ,
     &                 NAT,NAT0,NAT3,NX,NY,NZ)
      IMPLICIT NONE

c Arguments:

      INTEGER IDVOUT,MXN3,MXNAT,MXNX,MXNY,MXNZ,
     &        NAT,NAT0,NAT3,NX,NY,NZ
      CHARACTER CMETHD*6,CWHAT*1
      INTEGER*2 IOCC(MXNAT)
      REAL AK(3),DX(3)
      COMPLEX CXADIA(MXN3),CXAOFF(MXN3),
     &        CXX(MXN3),CXY(MXN3),CXZC(NX+1,NY+1,NZ+1,6),
     &        CXZW(MXNX*MXNY*MXNZ,*)

c Local variables:

      INTEGER J1,J2,J3
      REAL AKD

C Intrinsic functions:

      INTRINSIC CONJG

C***********************************************************************
c Subroutine CPROD

C Given:
C   AK(1-3)= k(1-3)*d , where k=k vector in vacuo
C                       d=effective lattice spacing=(dx*dy*dz)**(1/3)
C   DX(1-3)=(dx/d, dy/d, dz/d) 
C                 where dx,dy,dz=lattice spacing in x,y,z directions.
C                 By definition of effective lattice spacing d, must 
C                 have DX(1)*DX(2)*DX(3)=1
C   CMETHD = character variable specifying method used for FFT evaluation
C            (see ESELF)
C   CWHAT  = 'N' or 'C' determines what is to be done (see below)
C   CXADIA(1-NAT3) = diagonal elements of A(j,k)
C   CXAOFF(1-NAT3) = off-diagonal elements of 3x3 diagonal blocks of 
C                 A(j,k): CXAOFF(J,1-3)=a_{23},a_{31},a_{12} from 3x3
C                 matrix a_{ij} for dipole J.
C                 Recall that a_{ij} is inverse of 3x3 polarizability
C                 tensor alpha_{ij} for dipole J.
C   CXX(1-NAT3) = complex vector
C   CXZC        = (NX+1)*(NY+1)*(NZ+1)*6 array of Green function 
C                 coefficients used internally by ESELF, and generated
C                 by ESELF each time called with new k vector
C                 *****Not to be overwritten between calls to ESELF*****
C   CXZW        = complex workspace required by ESELF
C   IDVOUT      = device number for output
C   IOCC(1-NAT) = 0 or 1 if site is unoccupied or occupied
C   MXN3,MXNX,MXNY,MXNZ = dimensioning information (see DDSCAT)
C   NAT         = number of sites in extended target
C   NAT0        = number of occupied sites in extended target
C   NAT3        = 3*NAT
C   NX,NY,NZ    = extended target size is (NX*DX(1)) by (NY*DX(2)) by
C                 (NZ*DX(3))
C
C If CWHAT = 'N', this returns CXY(j)=A(j,k)*CXX(k) (sum over k)
C          = 'C', this returns CXY(j)=conjg(A(k,j))*CXX(k) (sum over k)
C
C It is assumed that matrix A(j,k) is symmetric
C Diagonal elements of A(j,k) are in vector CXADIA(j)
C Off-diagonal terms coupling dipole component x with direction y, etc
C are stored in vector CXAOFF(j), CXAOFF(j+NAT), CXAOFF(j+2*NAT)
C Other off-diagonal elements of A(j,k) are produced internally by
C   subroutine ESELF, which uses 3d-FFT to accomplish convolution
C   of CXX with A-matrix terms coupling dipole j with other dipoles
C   (i.e., E field at j due to other dipoles)
C
C Original version of CPROD created by P.J.Flatau, Colorado State Univ.
C Subsequently modified by B.T.Draine, Princeton Univ. Obs.
C History:
C 90.11.04 (BTD): Major modification to incorporate FFT method (FFT
C                 computations are handled by ESELF).
C 90.11.08 (BTD): Remove DO...ENDDO structures; renumber DO...CONTINUEs
C 90.11.21 (BTD): Remove "CCGMVM" option: only FFT stuff retained.
C 90.12.03 (BTD): No longer necessary to change ordering of CXX and CXY
C 90.12.05 (BTD): Changed call to ESELF to eliminate MXNX,MXNY,MXNZ
C 97.11.01 (BTD): Added comments and modified to accomodate noncubic
C                 lattice. Added DX to argument list for CPROD and to
C                 argument list in call to ESELF
C 97.12.25 (BTD): Added CXAOFF to argument list.
C 97.12.28 (BTD): Modified to use CXAOFF in calculation of CXY
C 98.01.01 (BTD): Modified coding to use indices J1,J2,J3 in loops
C                 computing with CXAOFF
C 00.06.13 (BTD): Explicit declaration of all variables.
C end history
C Copyright (C) 1993,1997,1998,2000
C               B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
c-----------------------------------------------------------------------

      AKD=SQRT(AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3))

C We assume that input vector CXX has zeroes for elements corresponding
C to vacuum sites.

      IF(CWHAT.EQ.'N')THEN
C***********************************************************************
C
         CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,DX,CXZC,CXZW,CXY)
C
C***********************************************************************
C ESELF omitted contribution from diagonal terms: add these.
C Also remember that ESELF computes -A_jk*x_k
         DO J1=1,NAT3
            CXY(J1)=CXADIA(J1)*CXX(J1)-CXY(J1)
         ENDDO

C Add contribution from off-diagonal elements of 3x3 diagonal blocks
C Version below assumes that off-diagonal elements a_ij are stored as
C CXAOFF(J,1-3)=(a_23,a_31,a_12) for dipole J.
C 97.12.28 (BTD):

         DO J1=1,NAT
            J2=J1+NAT
            J3=J2+NAT
            CXY(J1)=CXY(J1)+CXAOFF(J2)*CXX(J3)+CXAOFF(J3)*CXX(J2)
            CXY(J2)=CXY(J2)+CXAOFF(J3)*CXX(J1)+CXAOFF(J1)*CXX(J3)
            CXY(J3)=CXY(J3)+CXAOFF(J1)*CXX(J2)+CXAOFF(J2)*CXX(J1)
         ENDDO

      ELSEIF(CWHAT.EQ.'C')THEN

C Need to temporarily replace CXX by conjg(CXX)

         DO J1=1,NAT3
            CXX(J1)=CONJG(CXX(J1))
         ENDDO

C**********************************************************************
C
         CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,DX,CXZC,CXZW,CXY)
C
C************************************************************************
C
C ESELF omitted contribution from diagonal terms: add these.
C Also remember that ESELF computed -A_jk*conjg(x_k), while
C we seek to compute conjg(A_kj)* x_k.  Since A_kj=A_jk, we have
C conjg(A_kj)*x_k=conjg(A_jk*conjg(x_k))
C
         DO J1=1,NAT3
            CXY(J1)=CXADIA(J1)*CXX(J1)-CXY(J1)
         ENDDO
C
C Add contribution from off-diagonal elements of 3x3 diagonal blocks
C It is assumed that we have stored off-diagonal elements a_ij as
C CXAOFF(J,1-3)=(a_23,a_31,a_12) for element J
C Note that we could *probably* speed up the code by
C (1) "reducing" product vector CXY immediately after calling ESELF
C (2) restricting CXADIA to occupied sites
C (3) restricting CXAOFF to occupied sites
C but this would require modifications outside this routine, and
C therefore requires further study...
C 97.12.28 (BTD)

         DO J1=1,NAT
            J2=J1+NAT
            J3=J2+NAT
            CXY(J1)=CXY(J1)+CXAOFF(J2)*CXX(J3)+CXAOFF(J3)*CXX(J2)
            CXY(J2)=CXY(J2)+CXAOFF(J3)*CXX(J1)+CXAOFF(J1)*CXX(J3)
            CXY(J3)=CXY(J3)+CXAOFF(J1)*CXX(J2)+CXAOFF(J2)*CXX(J1)
         ENDDO

C Now compute conjugate of CXY, and restore CXX:

         DO J1=1,NAT3
            CXY(J1)=CONJG(CXY(J1))
            CXX(J1)=CONJG(CXX(J1))
         ENDDO
      ELSE
         WRITE(IDVOUT,FMT=*)' Error in CPROD: CWHAT= ',CWHAT
      ENDIF

C*** Now must "clean" product vector CXY:
C    (zero out elements corresponding to vacuum sites)

      IF(NAT0.LT.NAT)CALL NULLER(CXY,IOCC,MXNAT,MXN3,NAT)
      RETURN
      END

      SUBROUTINE CXC3DFFT(CZ,I1,I2,I3,I4,I5,I6,I7)
      IMPLICIT NONE
      COMPLEX CZ(1)
      INTEGER I1,I2,I3,I4,I5,I6,I7
C***********************************************************************
C    Note: This module should NOT be linked into DDSCAT when
C    running on a Convex with the vector library -- in that case
C    link with the -lveclib flag and use the native Convex C3DFFT
C    routine.
C    On other systems, however, you need to link in this routine
C    to avoid an error message from the compiler.
C
C History:
C 93.xx.xx (BTD) created to substitute for convex C3DFFT library routine
C 00.07.05 (BTD) renamed for greater Makefile consistency
C end history
C Copyright (C) 1993,2000 
C               B.T. Draine and P.J. Flatau
C
C This code is covered by the GNU General Public License.
C***********************************************************************
      CALL ERRMSG('FATAL','C3DFFT',' CONVEX option but dummy routine!')
      RETURN
      END

      SUBROUTINE CXFFTW(CX,MX,MY,MZ,ISIGN)
      IMPLICIT NONE
C
C Arguments:
C
      INTEGER ISIGN,MX,MY,MZ
      COMPLEX CX(MX,MY,MZ)
C***********************************************************************
C
C Purpose: This is a dummy routine to substitute for CXFFTW to allow
C          DDSCAT to be used on systems where the FFTW package has not
C          been installed.  If called, it will generate a fatal error
C          with an error message reporting DDSCAT has been compiled
C          with this dummy routine and option FFTW is unavailable.
C          
C History:
C 00.07.05 (BTD): created
C 03.07.13 (BTD): changed option FFTWFJ to option FFTW21
C end history
C Copyright (C) 2000, 2003
C               B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C
C***********************************************************************
C
      WRITE(0,*)'FATAL ERROR: DDSCAT has been compiled with dummy',
     &          ' version of CXFFTW'
      WRITE(0,*)' *** option FFTW21 cannot be used with dummy routine'
      WRITE(0,*)' *** to enable option FFTW21 it is necessary to:'
      WRITE(0,*)' --- have fftw 2.1.x library installed on system'

      WRITE(0,*)' --- in Makefile, variable LIBFFTW should include',
     &          ' correct path to libsfftw.a'
      WRITE(0,*)' --- in Makefile, define cxfftw = cxfftw',
     &          ' (rather than dummycxfftw.f)'
      WRITE(0,*)' --- make ddscat'
      STOP
      END



















      SUBROUTINE 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,DX,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***********************************************************************
C Dummy version of subroutine WRITENET
C (genuine subroutine WRITENET writes NetCDF portable binary file
C to file cnetfile 'dd.cdf')
C
C This dummy version is provided in order that use of DDSCAT is
C possible even if NetCDF software is not installed on system.
C
C History:
C 96.11.15 (PJF): First version of NetCDF  (Wolff's request)
C                 5a6 ---> 5a7 Converted to subroutine.
C 96.11.20 (BTD): dummy version created from real WRITENET
C 98.01.11 (BTD): added DX(3) to argument list for compatibility with
C                 new version of WRITENET
C end history
C Copyright (C) 1996,1998 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 IBETA,ICTHM,IDNC,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,CMDFFT*6,CMDTRQ*6,CNETFLAG*6,CSHAPE*6,
     &          CNETFILE*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),DX(3),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),
     &     SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA),
     &     TIMERS(MXTIMERS),WAVEA(MXWAV)
      INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3)
C
C This subroutine should never be called, so if it is, bail out...
C
      CALL ERRMSG('FATAL','DUMMYWRITENET',
     &   ' NetCDF is not installed. Change to NOTCDF in ddscat.par!')
      RETURN
      END

      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_null
C
C This version of timeit is a dummy which does not provide any
C timing information.
C 94.06.20 (PJF) add dtime to the formal parameters
C 99.05.12 (BTD) explicitly declare DTIME
C
C Copyright (C) 1993,1999 B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
C Arguments:
      CHARACTER CMSGTM*(*)
      REAL DTIME
C Local variables:
      CHARACTER CSTA*3,CMSGNM*70
C***********************************************************************
      DATA CSTA/'ONE'/
      IF(CSTA.EQ.'ONE')THEN
         CSTA='TWO'
      ELSEIF(CSTA.EQ.'TWO')THEN
         CSTA='ONE'
         WRITE(CMSGNM,FMT='(A,A)')'Timing results for: ',CMSGTM
         CALL WRIMSG('TIMEIT',CMSGNM)
         IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN
            WRITE(CMSGNM,FMT='(A)')
     &         'Dummy timing routine: no timings available'
            CALL WRIMSG('TIMEIT',CMSGNM)
         ENDIF
      ENDIF
      RETURN
      END
