      SUBROUTINE INTERP(X,Y,Z,XA,YA,ZA,IDVOUT,MXTAB,NTAB)
C Arguments:
      INTEGER IDVOUT,MXTAB,NTAB
      REAL X,Y,Z
      REAL XA(MXTAB),YA(MXTAB),ZA(MXTAB)
C Local variables:
      INTEGER I1,I2,I3,I4,INC
      REAL SGN,X1,X2,X3,X4,Y1,Y2,Y3,Y4
C***********************************************************************
C Given array of tabulated values XA(1-NTAB), YA(1-NTAB), ZA(1-NTAB),
C and given independent variable X, this routine interpolates for two
C dependent variables: Y and Z.
C Uses only parabolic interpolation for safety.
C When called with X outside range of tabulated values XA, returns
C extrapolated values Y and Z but issues warning statement.
C Note: XA values must be either monotonically increasing OR
C       monotonically decreasing.
C B.T.Draine, Princeton Univ. Observatory
C 89/12/19 (BTD): Revised.
C 91/05/02 (BTD): Added IDVOUT to argument list.
C                 Changed WRITE(0 -> WRITE(IDVOUT
C 92/03/24 (BTD): Added test IF(I2.GE.NTAB) in case INTERP is used for
C                 data sets with different lengths.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      DATA I2/1/
      IF(I2.GE.NTAB)I2=NTAB-1
      INC=0
      SGN=1.
C*** Check whether X is increasing or decreasing:
      IF(XA(1).GE.XA(NTAB))SGN=-1.
C*** Check whether outside table limits
      IF((SGN*(X-XA(1))).LT.0.)THEN
         I2=2
         WRITE(IDVOUT,6990)X
         GOTO 4600
      ENDIF
      IF(SGN*(X-XA(NTAB)).GT.0.)THEN
         I2=NTAB-1
         WRITE(IDVOUT,6990)X
         GOTO 4600
      ENDIF
C*** X is within table limits.  Find I2
 1000 IF(SGN*(XA(I2)-X))2000,3000,4000
 2000 IF(INC)2100,2200,2200
 2100 IF(I2+2-NTAB)4700,4700,2500
 2200 INC=1
      I2=I2+1
      IF(I2+1.LE.NTAB)GO TO 1000
 2500 I2=NTAB-1
      GO TO 4600
 3000 Y=YA(I2)
      Z=ZA(I2)
      RETURN
 4000 IF(INC)4200,4200,4100
 4100 I2=I2-1
      IF(I2-2)4500,4700,4700
 4200 INC=-1
      I2=I2-1
      IF(I2.GE.2)GOTO 1000
 4500 I2=2
 4600 I1=I2-1
      I3=I2+1
      X1=XA(I1)
      X2=XA(I2)
      X3=XA(I3)
      Y1=YA(I1)
      Y2=YA(I2)
      Y3=YA(I3)
      CALL PARAB3(X,Y,X1,X2,X3,Y1,Y2,Y3)
      Y1=ZA(I1)
      Y2=ZA(I2)
      Y3=ZA(I3)
      CALL PARAB3(X,Z,X1,X2,X3,Y1,Y2,Y3)
      RETURN
 4700 I1=I2-1
      I3=I2+1
      I4=I2+2
      X1=XA(I1)
      X2=XA(I2)
      X3=XA(I3)
      X4=XA(I4)
      Y1=YA(I1)
      Y2=YA(I2)
      Y3=YA(I3)
      Y4=YA(I4)
      CALL PARAB4(X,Y,X1,X2,X3,X4,Y1,Y2,Y3,Y4)
      Y1=ZA(I1)
      Y2=ZA(I2)
      Y3=ZA(I3)
      Y4=ZA(I4)
      CALL PARAB4(X,Z,X1,X2,X3,X4,Y1,Y2,Y3,Y4)
      RETURN
 6990 FORMAT('Warning from INTERP: outside table limits for X=',1PE12.5)
      END

      SUBROUTINE PARAB3(X,Y,X1,X2,X3,Y1,Y2,Y3)
      REAL A,B,X,X1,X2,X3,Y,Y1,Y2,Y3
C***********************************************************************
C Subroutine PARAB3 does parabolic interpolation, with parabola
C constrained to fit (x1,y1),(x2,y2),(x3,y3) exactly.
C B.T.Draine, Institute for Advanced Study, March 1980.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
C***********************************************************************
      A=(Y3-Y2-(X3-X2)*(Y1-Y2)/(X1-X2))/((X3-X2)*(X3-X1))
      B=(Y1-Y2)/(X1-X2)-A*(X1-X2)
      Y=(A*(X-X2)+B)*(X-X2)+Y2
      RETURN
      END

      SUBROUTINE PARAB4(X,Y,X1,X2,X3,X4,Y1,Y2,Y3,Y4)
      REAL A,B,X,X1,X2,X3,X4,Y,Y1,Y2,Y3,Y4
C***********************************************************************
C Subroutine PARAB4 is designed to do parabolic interpolation, with
C parabola constrained to match (x2,y2) and (x3,y3) exactly, and to
C minimize sum of squared deviations from (x1,y1) and (x4,y4).
C It is assumed that x1.lt.x2.le.x.lt.x3.lt.x4
C Fit parabola Y=A*(X-X2)**2+B*(X-X2)+Y2
C B.T.Draine, Institute for Advanced Study, March 1980.
C
C Copyright (C) 1993, B.T. Draine and P.J. Flatau
C This code is covered by the GNU General Public License.
***********************************************************************
      A=((X1-X2)*(Y3-Y2)/(X3-X2)+Y2-Y1)*(X1-X2)*(X1-X3)
     &  +((X4-X2)*(Y3-Y2)/(X3-X2)+Y2-Y4)*(X4-X2)*(X4-X3)
      A=-A/(((X1-X2)*(X1-X3))**2+((X4-X2)*(X4-X3))**2)
      B=(Y3-Y2)/(X3-X2)-A*(X3-X2)
      Y=(A*(X-X2)+B)*(X-X2)+Y2
      RETURN
      END
      SUBROUTINE NAMER(IWAV,IRAD,IDIR,CFLSCA,CFLAVG)
C***********************************************************************
C Purpose: to generate file names for specific (wavelength,
C                                               size,
C                                               incident direction)
C Present version allows up to 100 wavelengths, (00-99)
C                              100 sizes        (00-99)
C                             1000 directions   (000-999)
C Given:
C         IWAV=integer (0-99) identifying wavelength
C         IRAD=integer (0-99) identifying size
C         IDIR=integer (0-999) identifying direction
C Returns:
C          CFLSCA=name for output file containing Qext,Qabs,Qpha,Qsca,g
C                 and fml values for selected directions
C          CFLAVG=name for output file for given wavelength and size,
C                 containing scattering etc. averaged over inc.dir.
C
C B.T.Draine, Princeton Univ. Observatory, 1988
C History:
C 90/11/21 (BTD): Modified to allow 3 digits to identify direction
C                 and to change indices to begin from w00r00k000.sca
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:
      INTEGER IDIR,IDIR0,IRAD,IRAD0,IWAV,IWAV0
      CHARACTER CFLSCA*14,CFLAVG*13
C
C Local variables:
      CHARACTER A1*1,A2*1,A3*1,EXT*4,ZSCA*14,ZAVG*13
      CHARACTER N(10)*1,V1(2)*1,V2(2)*1,V3(3)*1,W(14)*1
C
C Equivalences:
      EQUIVALENCE (A1,W(1)),(V1,W(2)),(A2,W(4)),(V2,W(5)),(A3,W(7)),
     $            (V3,W(8)),(EXT,W(11))
      EQUIVALENCE (W,ZSCA)
      EQUIVALENCE (W,ZAVG)
C
C DATA statements:
      DATA N/'0','1','2','3','4','5','6','7','8','9'/
C***********************************************************************
C*** Filename CFLSCA will be of the form w01r01k001
C     (for IWAV=2,IRAD=2,IDIR=2)
C
      A1='w'
      A2='r'
      A3='k'
C*** Set IWAV0=IWAV-1 to run over range 0-99
C    Set IRAD0=IRAD-1 to run over range 0-99
C    Set IDIR0=IDIR-1 to run over range 0-999
      IWAV0=IWAV-1
      IRAD0=IRAD-1
      IDIR0=IDIR-1
      V1(1)=N(IWAV0/10+1)
      V1(2)=N(IWAV0-10*(IWAV0/10)+1)
      V2(1)=N(IRAD0/10+1)
      V2(2)=N(IRAD0-10*(IRAD0/10)+1)
      V3(1)=N(IDIR0/100+1)
      V3(2)=N((IDIR0-100*(IDIR0/100))/10+1)
      V3(3)=N(IDIR0-10*(IDIR0/10)+1)
      EXT='.sca'
      CFLSCA=ZSCA
C*** CFLAVG is of slightly different form,
C    e.g., w00r00ori.avg
      A3='o'
      V3(1)='r'
      V3(2)='i'
      V3(3)='.'
      EXT='avg '
      CFLAVG=ZAVG
      RETURN
C
      END
      SUBROUTINE NULLER(CXVEC,IOCC,MXNAT,MXN3,NAT)
C Arguments:
      INTEGER MXNAT,MXN3,NAT
      INTEGER*2 IOCC(MXNAT)
      COMPLEX*8 CXVEC(NAT,3)
C Local variables:
      INTEGER J,L
C***********************************************************************
C This routine takes as input a complex vector with NAT3 elements
C and sets to zero those elements corresponding to "vacuum" sites.
C This is accomplished by multiplying by vector IOCC(J), whose
C elements are either 1 or 0 depending on whether site is occupied
C or "vacuum".
C Note: CXVEC should be given dimension CXVEC(NAT,3) here
C       so that first 3*NAT elements of CXVEC are employed.
C
C B.T.Draine, Princeton Univ. Obs., 90/11/1
C History:
C 90/12/15 (BTD): Corrected error in dimensioning of CXVEC.
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 110 L=1,3
         DO 100 J=1,NAT
            CXVEC(J,L)=IOCC(J)*CXVEC(J,L)
  100    CONTINUE
  110 CONTINUE
      RETURN
      END
      SUBROUTINE ORIENT(BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX,
     &                  MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI,
     &                  BETA,THETA,PHI,WGTA,WGTB)
C Arguments:
      INTEGER MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI
      REAL BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX
      REAL BETA(MXBETA),THETA(MXTHET),PHI(MXPHI),
     &     WGTA(MXTHET,MXPHI),WGTB(MXBETA)
C Local variables:
      INTEGER I,J
      REAL DELTA,WG
C***********************************************************************
C Given:
C        BETAMI=minimum value of beta (radians)
C        BETAMX=maximum value of beta (radians)
C        THETMI=minimum value of theta (radians)
C        THETMX=maximum value of theta (radians)
C        PHIMIN=minimum value of phi (radians)
C        PHIMAX=maximum value of phi (radians)
C        MXBETA,MXTHET,MXPHI=parameters for dimensioning of arrays
C                             BETA,THETA,PHI
C        NBETA=desired number of values of beta
C        NTHETA=desired number of values of theta
C        NPHI=desired number of values of PHI
C Returns:     
C        BETA(1-NBETA)=beta values (radians)
C        THETA(1-NTHETA)=theta values (radians)
C        PHI(1-NPHI)=phi values (radians)
C        WGTA(1-NTHETA,1-NPHI)=weighting of each orientation of target
C                              axis a1 in Lab Frame
C                              (sum of weights = 1)
C        WGTB(1-NBETA)=weighting of each rotation of target around a1
C                              (sum of weights = 1)
C            Note: it is assumed that target orientation weight function
C                  can be factored into WGTA*WGTB -- i.e., that rotations
C                  around a1 are decoupled from orientation of a1.
C
C Purpose: to generate list of desired target orientations
C
C Definitions: beta=angle of rotation of target around target axis a1
C                   (beta=0 is defined to be such that a2 lies in a1-x
C                    plane, with a2 in direction of increasing theta)
C              theta=angle between target axis a1 and x axis
C                   (0.le.theta.le.pi)
C              phi=angle of rotation of axis a1 around axis
C              (phi=0 is defined to be such that a1 lies in xy plane)
C Present version assumes:
C        beta to be uniformly distributed between BETAMI and BETAMX
C        cos(theta) to be uniformly distributed between cos(THETMI) and
C                   cos(THETMX)
C        phi to be uniformly distributed between PHIMIN and PHIMAX
C
C Values assigned to BETA and PHI are midpoints of uniform intervals in
C beta and phi
C Note that if NBETA=1, BETA is set to midpoint of range in phi
C Likewise, if NPHI=1, PHI is set to midpoint of range in phi
C
C If NTHETA=1: specify a single value of THETA
C        cos(theta)=0.5*(cos(thetami)+cos(thetamx))
C
C If NTHETA>1:
C    If NTHETA is even: 
C       divide range of cos(theta) into NTHETA equal intervals
C       set THETA to midpoints of these intervals
C       give intervals equal weights
C    If NTHETA is odd:
C       divide range of cos(theta) into NTHETA-1 equal intervals
C       set THETA to endpoints of these intervals; first value of
C       THETA is THETAMI and last value of THETA is THETAMX
C
C Note: values chosen for beta, phi are at midpoints of
C        uniform intervals
C
C B.T.Draine, Princeton Univ. Obs., 89.11.20
C History:
C 92.04.01 (BTD) Modify to allow Simpson's rule or trapezoidal rule
C                for integration over cos(theta).
C 99.03.01 (BTD) Modified to properly assign THETA values when NPHI>1
C                Code up to this time was not consistent with
C                description in UserGuide.  This problem was brought
C                to my attention by Miroslav Kocifaj
C                (Miroslav.Kocifaj@swh.sk).
C 99.03.05 (BTD) Corrected errors in assignment of THETA values.
C                Thanks to Miroslav Kocifaj for detecting these.
C 99.04.30 (BTD) Corrected errors in assignment of BETA and
C                PHI values when BETAMI.ne.0 or PHIMIN.ne.0
C                These errors (inadvertent omission of BETAMI and
C                PHIMIN from expressions assigning BETA(J) and PHI(J))
C                were introduced in 99.03.05 revision.
C                Thanks to Henriette Lemke for calling attention to
C                these errors.
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
C Assign values to BETA:
C
      DELTA=(BETAMX-BETAMI)/REAL(NBETA)
      DO J=1,NBETA
         BETA(J)=BETAMI+DELTA*(REAL(J)-0.5)
      ENDDO
C
C Assign values to THETA:
C
      IF(NTHETA.EQ.2*(NTHETA/2))THEN
C
C THETA is even:
C
        DELTA=(COS(THETMX)-COS(THETMI))/REAL(NTHETA)
        DO J=1,NTHETA
           THETA(J)=ACOS(COS(THETMI)+DELTA*(REAL(J)-0.5))
        ENDDO
      ELSE
C
C THETA is odd:
C
         IF(NTHETA.EQ.1)THEN
            THETA(1)=ACOS(0.5*(COS(THETMI)+COS(THETMX)))
         ELSE
            DELTA=(COS(THETMX)-COS(THETMI))/REAL(NTHETA-1)
            THETA(1)=THETMI
            THETA(NTHETA)=THETMX
            DO J=2,NTHETA-1
               THETA(J)=ACOS(COS(THETMI)+DELTA*REAL(J-1))
            ENDDO
         ENDIF
      ENDIF
C
C Assign values to PHI:
C
      DELTA=(PHIMAX-PHIMIN)/REAL(NPHI)
      DO J=1,NPHI
         PHI(J)=PHIMIN+DELTA*(REAL(J)-0.5)
      ENDDO
C
C** Specify weight factors WGTA, WGTB
C   (Note: weight function WGTA = 4*pi*P/(NTHETA*NPHI), where
C    P=(probability/solid angle) of orientation in direction THETA,PHI .
C    The orientational averaging program automatically samples uniformly
C    in cos(theta) to allow for d(solid angle)=sin(theta)d(theta)d(phi).
C
C   Present version assumes random orientations.
C
C   When NTHETA is >1 and even, we use trapezoidal integration
C   When NTHETA is >1 and odd, we use Simpson's rule integration.
C
      IF(NTHETA.EQ.1.OR.NTHETA.EQ.2*(NTHETA/2))THEN
C
C either 1 or even number of theta values: midpoints of intervals
C
         WG=1./REAL(NTHETA*NPHI)
         DO J=1,NPHI
            DO I=1,NTHETA
               WGTA(I,J)=WG
            ENDDO
         ENDDO
      ELSE
C
C odd number >1 of theta values: use Simpson's rule weighting
C
         WG=1./REAL(3*(NTHETA-1)*NPHI)
         DO J=1,NPHI
            WGTA(1,J)=WG
            WGTA(NTHETA,J)=WG
         ENDDO
         WG=4./REAL(3*(NTHETA-1)*NPHI)
         DO I=2,NTHETA,2
            DO J=1,NPHI
               WGTA(I,J)=WG
            ENDDO
         ENDDO
         IF(NTHETA.GE.5)THEN
            WG=2./REAL(3*(NTHETA-1)*NPHI)
            DO I=3,NTHETA-1,2
               DO J=1,NPHI
                  WGTA(I,J)=WG
               ENDDO
            ENDDO
         ENDIF
      ENDIF
C
C Assign weights for rotation angle beta:
C
      WG=1./REAL(NBETA)
      DO I=1,NBETA
         WGTB(I)=WG
      ENDDO
C
      RETURN
      END
      SUBROUTINE REAPAR(cbinflag, cnetflag, AEFFA,
     &                  BETAMI,BETAMX,
     &                  CALPHA,CDIEL,CFLEPS,CFLPAR,CMDSOL,CMDFFT,CSHAPE,
     &                  CMDTRQ,CXE01,CXE02,
     &                  DEGRAD,
     &                  ICTHM,IDVERR,IDVOUT,INIT,IOPAR,IORTH,IPHM,
     &                  IWRKSC,
     &                  MXBETA,MXCOMP,MXMEM,MXPHI,MXRAD,MXSCA,MXTHET,
     &                  MXWAV,
     &                  NBETA,NCOMP,NPHI,NRAD,NSCA,NTHETA,NWAV,
     &                  PHIMIN,PHIMAX,PHIN,
     &                  THETAN,THETMI,THETMX,TOL,
     &                  WAVEA,SHPAR)
C***********************************************************************
C Subroutine REAPAR handles the reading of input parameters from the
C "ddscat.par" file, as well as elementary processing with those input
C parameters to generate arrays.
C Original version created by P.J.Flatau, Colorado State Univ.
C Modified by B.T.Draine, Princeton Univ. Obs.
C History:
C 90.11.09 (BTD): Modified so as to remove sensitivity to errors
C                 in sign of entered value of DTHETA.
C 90.11.30 (BTD): Added test to compare CMDFFT and MXMEM
C 90.12.14 (BTD): Changed option MKRECT to RCTNGL
C 91.01.03 (BTD): Added MXBETA,MXPHI,MXTHET to argument list
C                 and added checking of NBETA,NPHI,NTHETA
C 91.05.08 (BTD): Added CALPHA to allow choice of prescription
C                 for polarizabilities
C 91.05.08 (BTD): Some rewriting done.
C 91.05.14 (BTD): Added separate options LDRXYZ and LDRAVG
C 91.05.15 (BTD): Added new option LDR000
C 91.05.23 (BTD): Added variable IWRKSC to argument list and
C                 added code to read IWRKSC
C 91.05.30 (BTD): Changed LDRXYZ->LDR100
C                 Added option LDR111
C 91.11.13 (BTD): Now three options: DRAI88, GOBR88, LATTDR
C                 (subroutine ALPHA now allows for directional
C                 dependence for LATTDR option)
C 92.05.14 (BTD): Allow for experimental option LDRISO
C 92.09.21 (BTD): Allow for option ANIELL (ellipsoid of anisotropic
C                 material); change option GRPHIT -> UNICYL (cylinder
C                 of uniaxial material)
C 93.01.07 (BTD): Allow for options TWOELL and TWOAEL (two touching
C                 ellipsoids, each either isotropic or anisotropic)
C 93.03.11 (BTD): Deleted all code associated with variable CMACHN
C                 which previously was used to identify machine/OS
C                 This information was not used or needed.
C 93.12.15 (BTD): Added code to check for compatibility of NCOMP
C                 when CSHAPE=TWOELL,THRELL,TWOAEL,or THRAEL
C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6)
C                 Added CONELL (two concentric ellipsoids) as target
C                 option.
C 94.06.20 (PJF): Added code to handle NEWTMP
C 94.12.02 (BTD): Added options 
C                    LB1991 (Dungey-Bohren 1991 polarizability)
C                    LR1994 (Lumme-Rahola 1994 polarizability)
C 95.05.26 (BTD): Corrected construction of orthogonal polarization
C                 state for case where E01 is complex (e.g.,
C                 circular polarization)
C 95.06.14 (BTD): Added code to check that user is not requesting
C                 f_ml when illuminating target with other then
C                 linearly polarized light
C 95.06.19 (BTD): Added argument CMDSOL to allow selection of
C                 different CCG methods
C                 Added argument CMDTRQ to allow torque calculations
C                 to be skipped.
C 96.01.05 (BTD): Added shape options BLOCKS and DW1996
C 96.01.25 (BTD): Added shape option TWOSPH
C 96.01.25 (BTD): Disable code added 95.06.14 so that user can
C                 calculate f_ml even when illuminating target with
C                 elliptically polarized light.
C 96.10.18 (BTD): Changed option NEWTMP to GPFAFT (Generalized Prime
C                 Factor Algorithm for Fourier Transform).
C 96.11.06 (BTD): Added provision to 
C                 (1) produce fatal error if CXE01(1).NE.(0.,0.)
C                 (2) normalize CXE01 if necessary
C 96.11.14 (PJF): cbinflag for binary file
C 96.11.15 (PJF): netCDF netflag
C 96.11.21 (BTD): modified handling of CBINFLAG and CNETFLAG to
C                 be consistent with other option flags
C 98.12.21 (BTD): changed dimension of CFLPAR from CHARACTER*40 to
C                 CHARACTER*60 to allow longer file names
C                 (also changed in DDSCAT.f and dielec.f)
C end history
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***********************************************************************
C
C Arguments:
      CHARACTER
     &   CALPHA*(*),CDIEL*(*),CFLPAR*(*),CMDSOL*(*),CMDFFT*(*),
     &   CMDTRQ*(*),CSHAPE*(*), cbinflag*(*), cnetflag*(*)
      INTEGER
     &   ICTHM,IDVERR,IDVOUT,INIT,IOPAR,IORTH,IPHM,IWRKSC,
     &   MXBETA,MXCOMP,MXMEM,MXPHI,MXRAD,MXSCA,MXTHET,MXWAV,
     &   NBETA,NCOMP,NPHI,NRAD,NSCA,NWAV,NTHETA
      CHARACTER*60 CFLEPS(MXCOMP)
      REAL
     &   BETAMI,BETAMX,DEGRAD,PHIMIN,PHIMAX,THETMI,THETMX,TOL
      REAL
     &   AEFFA(MXRAD),PHIN(MXSCA),SHPAR(6),THETAN(MXSCA),WAVEA(MXWAV)
      COMPLEX CXE01(3),CXE02(3)
C
C Local variables:
      CHARACTER CDIVID*3,CLINE*70,CMSGNM*70
      INTEGER J,NSCA0
      REAL
     &   AEFEND,AEFINI,DELTA,DTHETA,E1,PHI1,THETA1,THETA2,
     &   WAVEND,WAVINI
C
C External routines:
      EXTERNAL DIVIDE,ERRMSG,WRIMSG
C
C***********************************************************************
      OPEN (UNIT=IOPAR,FILE=CFLPAR,STATUS='OLD')
C***********************************************************************
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG('REAPAR',CLINE)
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
C
C    Specify whether torques are to be calculated:
C
C    CMDTRQ*6 = 'DOTORQ' -  calculate torque on grain
C             = 'NOTORQ' -  do not calculate torque on grain
C
      READ(IOPAR,FMT=*,ERR=40)CMDTRQ
      IF(CMDTRQ.EQ.'dotorq')CMDTRQ='DOTORQ'
      IF(CMDTRQ.EQ.'notorq')CMDTRQ='NOTORQ'
      IF(CMDTRQ.EQ.'DOTORQ')THEN
         WRITE(CMSGNM,FMT='(A,A)')CMDTRQ,' - compute torques '
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSEIF(CMDTRQ.EQ.'NOTORQ')THEN
         WRITE(CMSGNM,FMT='(A,A)')CMDTRQ,' - do not compute torques '
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMDTRQ')
      ENDIF
C
C    Define method used for iterative solution of complex linear equations
C
C    CMDSOL*6 = 'PETRKP' -  Petravic & Kuo-Petravic method
C             = 'PBCGST' -  PIM BiConjugate Gradient with Stabilization
C
      READ(IOPAR,FMT=*,ERR=40)CMDSOL
      IF(CMDSOL.EQ.'petrkp')CMDSOL='PETRKP'
      IF(CMDSOL.EQ.'pbcgst')CMDSOL='PBCGST'
      IF(CMDSOL.EQ.'PETRKP'.OR.
     &   CMDSOL.EQ.'PBCGST')THEN
         WRITE(CMSGNM,FMT='(A,A)')CMDSOL,' - CCG Method  '
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMDSOL')
      ENDIF
C
C    Define FFT method:
C
C    CMDFFT*6 = 'BRENNR' -  Use Brenner's FOURX FFT routine.
C             = 'TMPRTN' -  Use Temperton's FFT routine.
C             = 'GPFAFT' -  GPFA code of Temperton
C             = 'CONVEX' -  Use native Convex veclib routine C3DFFT.
C
      READ(IOPAR,FMT=*,ERR=40)CMDFFT
      IF(CMDFFT.EQ.'BRENNR'.OR.
     &   CMDFFT.EQ.'TMPRTN'.OR.
     &   CMDFFT.EQ.'CONVEX'.OR.
     &   CMDFFT.EQ.'GPFAFT')THEN
         WRITE(CMSGNM,FMT='(A,A)')CMDFFT,' - Solution Method  '
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMDFFT')
      ENDIF
C*** Check whether compiled with enough memory for TMPRTN option:
      IF(CMDFFT.EQ.'TMPRTN'.AND.MXMEM.EQ.0)THEN
         CALL ERRMSG('FATAL','REAPAR',
     &   ' Must recompile with MXMEM=1 to use TMPRTN option')
      ENDIF
      IF(CMDFFT.NE.'TMPRTN'.AND.MXMEM.EQ.1)THEN
         CALL ERRMSG('WARNING','REAPAR',
     &   ' You do not need MXMEM=1 with BRENNR,CONVEX, or GPFAFT')
      ENDIF
C***********************************************************************
C
C    Define prescription for computing polarizabilities:
C
C    CALPHA*6 = 'DRAI88' - Draine (1988) prescription
C             = 'GOBR88' - Goedecke & Obrien (1988) prescription
C             = 'LATTDR' - Lattice Dispersion Relation
C             = 'LDRISO' - Isotropized Lattice Dispersion Relation
C             = 'DB1991' - Dungey & Bohren (1991) prescription
C             = 'LR1994' - Lumme & Rahola (1994) prescription
C
      READ(IOPAR,FMT=*,ERR=40)CALPHA
      IF(CALPHA.EQ.'DRAI88')THEN
         CALL WRIMSG('REAPAR',
     &               'DRAI88 - Draine (1988) prescription for alpha')
      ELSEIF(CALPHA.EQ.'GOBR88')THEN
         CALL WRIMSG('REAPAR',
     &      'GOBR88 - Goedecke & OBrien (1988) prescription for alpha')
      ELSEIF(CALPHA.EQ.'LATTDR')THEN
         CALL WRIMSG('REAPAR',
     &      'LATTDR - Lattice Dispersion Relation for alpha')
      ELSEIF(CALPHA.EQ.'LDRISO')THEN
         CALL WRIMSG('REAPAR',
     &     'LATTDR - Isotropized Lattice Dispersion Relation for alpha')
      ELSEIF(CALPHA.EQ.'DB1991')THEN
         CALL WRIMSG('REAPAR',
     &     'DB1991 - Dungey & Bohren (1991) prescription for alpha')
      ELSEIF(CALPHA.EQ.'LR1994')THEN
         CALL WRIMSG('REAPAR',
     &     'LR1994 - Lumme & Rahola (1994) prescription for alpha')
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' wrong definition of CALPHA')
      ENDIF
C
C binary file flag
      READ(IOPAR,FMT=*,ERR=40)CBINFLAG
      IF(CBINFLAG.EQ.'allbin')CBINFLAG='ALLBIN'
      IF(CBINFLAG.EQ.'oribin')CBINFLAG='ORIBIN'
      IF(CBINFLAG.EQ.'notbin')CBINFLAG='NOTBIN'
      IF(CBINFLAG.EQ.'ALLBIN'
     &   .OR.CBINFLAG.EQ.'ORIBIN'
     &   .OR.CBINFLAG.EQ.'NOTBIN')THEN
         WRITE(CMSGNM,FMT='(A,A)')CBINFLAG,
     &                            ' - Unformatted binary dump option'
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' Wrong definition of CBINFLAG')
      ENDIF
C
C binary file flag
      READ(IOPAR,FMT=*,ERR=40)CNETFLAG
      IF(CNETFLAG.EQ.'allcdf')CNETFLAG='ALLCDF'
      IF(CNETFLAG.EQ.'oricdf')CNETFLAG='ORICDF'
      IF(CNETFLAG.EQ.'notcdf')CNETFLAG='NOTCDF'
      IF(CNETFLAG.EQ.'ALLCDF'  
     &   .OR.CNETFLAG.EQ.'ORICDF' 
     &   .OR.CNETFLAG.EQ.'NOTCDF')THEN
         WRITE(CMSGNM,FMT='(A,A)')CNETFLAG,
     &                            ' - netCDF file option'
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' Wrong definition of CNETFLAG')
      ENDIF
C
C***********************************************************************
C
C    Define shape:
C
C    CSHAPE*6 = 'FRMFIL'     read shape data from file
C             = 'RCTNGL'     make rectangular solid
C             = 'ELLIPS'     make ellipsoidal solid
C             = 'CYLNDR'     make cylindrical solid
C             = 'HEXGON'     make hexagonal prism
C             = 'TETRAH'     make regular tetrahedron
C             = 'UNICYL'     make cylinder of unixaxial material
C             = 'ANIELL'     make ellipsoid of anisotropic material
C             = 'TWOELL'     make two touching ellipsoids of materials 1,2
C             = 'TWOAEL'     make two touching anisotropic ellipsoids of
C                            materials 1-6
C             = 'THRELL'     make three touching isotropic ellipsoids of
C                            materials 1,2,3
C             = 'THRAEL'     make three touching anisotropic ellipsoids of
C                            materials 1-9
C             = 'CONELL'     make two concentric ellipsoids of
C                            materials 1,2
C             = 'BLOCKS'     make target out of cubic blocks defined by
C                            data in file 'blocks.par'
C             = 'DW1996'     make 13-cube target used by Draine &
C                            Weingartner 1996
C             = 'TWOSPH'     make two touching spheroids with symmetry axes at
C                            specified angle
C
      READ(IOPAR,FMT=*,ERR=40)CSHAPE
      IF(    CSHAPE.EQ.'FRMFIL'
     &   .OR.CSHAPE.EQ.'RCTNGL'
     &   .OR.CSHAPE.EQ.'ELLIPS'
     &   .OR.CSHAPE.EQ.'CYLNDR'
     &   .OR.CSHAPE.EQ.'HEXGON'
     &   .OR.CSHAPE.EQ.'TETRAH'
     &   .OR.CSHAPE.EQ.'UNICYL'
     &   .OR.CSHAPE.EQ.'ANIELL'
     &   .OR.CSHAPE.EQ.'TWOELL'
     &   .OR.CSHAPE.EQ.'TWOAEL'
     &   .OR.CSHAPE.EQ.'THRELL'
     &   .OR.CSHAPE.EQ.'THRAEL'
     &   .OR.CSHAPE.EQ.'CONELL'
     &   .OR.CSHAPE.EQ.'BLOCKS'
     &   .OR.CSHAPE.EQ.'DW1996'
     &   .OR.CSHAPE.EQ.'TWOSPH')THEN
         WRITE(CMSGNM,FMT='(A,A)')CSHAPE,' - Shape definition '
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' Wrong shape directive ')
      ENDIF
C
C Read shape parameters
C Need 3 shape parameters for all options except 
C      FRMFIL (here skip a line, later read target from file)
C      TETRAH, DW1996 (1 shape parameter)
C      CYLNDR, HEXGON, UNICYL (2 shape parameters)
C      TWOSPH (5 shape parameters)
C      CONELL (6 shape parameters)
C
      IF(CSHAPE.EQ.'FRMFIL')THEN
         READ(IOPAR,FMT=*,ERR=40)
      ELSEIF(
     &   CSHAPE.EQ.'TETRAH'.OR.
     &   CSHAPE.EQ.'DW1996')THEN
         READ(IOPAR,FMT=*,ERR=40)SHPAR(1)
      ELSEIF(
     &   CSHAPE.EQ.'CYLNDR'.OR.
     &   CSHAPE.EQ.'HEXGON'.OR.
     &   CSHAPE.EQ.'UNICYL')THEN
         READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2)
      ELSEIF(CSHAPE.EQ.'TWOSPH')THEN
         READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2),SHPAR(3),
     &                           SHPAR(4),SHPAR(5)
      ELSEIF(CSHAPE.EQ.'CONELL')THEN
         READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2),SHPAR(3),
     &                           SHPAR(4),SHPAR(5),SHPAR(6)
      ELSE
         READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2),SHPAR(3)
      ENDIF
C***********************************************************************
C
C   Obtain names of file(s) containing dielectric function(s)
C   NCOMP = number of different dielectric functions
C   CFLEPS(1-NCOMP) = names of files containing dielectric functions
C
      READ(IOPAR,FMT=*,ERR=40)NCOMP
C
C*** Check that MXCOMP is large enough *********************************
      IF(NCOMP.GT.MXCOMP)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP > MXCOMP: must increase MXCOMP in DDSCAT')
C
C*** Check that NCOMP=2 if CSHAPE=UNICYL *******************************
C                     3           ANIELL
C                     2           TWOELL
C                     3           THRELL
C                     6           TWOAEL
C                     9           THRAEL
C                     2           CONELL
C                     2           TWOSPH
      IF(CSHAPE.EQ.'UNICYL'.AND.NCOMP.NE.2)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 2 for option UNICYL')
      IF(CSHAPE.EQ.'ANIELL'.AND.NCOMP.NE.3)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 3 for option ANIELL')
      IF(CSHAPE.EQ.'TWOELL'.AND.NCOMP.NE.2)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 2 for option TWOELL')
      IF(CSHAPE.EQ.'THRELL'.AND.NCOMP.NE.3)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 3 for option THRELL')
      IF(CSHAPE.EQ.'TWOAEL'.AND.NCOMP.NE.6)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 6 for option TWOAEL')
      IF(CSHAPE.EQ.'THRAEL'.AND.NCOMP.NE.9)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 9 for option THRAEL')
      IF(CSHAPE.EQ.'CONELL'.AND.NCOMP.NE.2)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 2 for option CONELL')
      IF(CSHAPE.EQ.'TWOSPH'.AND.NCOMP.NE.2)
     &   CALL ERRMSG('FATAL','REAPAR',
     &   'NCOMP must be 2 for option TWOSPH')
C***********************************************************************
      READ(IOPAR,FMT=*,ERR=40)CDIEL
      IF(CDIEL.EQ.'TABLES')THEN
         DO 1000 J=1,NCOMP
             READ(IOPAR,FMT=*,ERR=40)CFLEPS(J)
 1000    CONTINUE
      ELSEIF(CDIEL.EQ.'H2OICE'.OR.CDIEL.EQ.'H2OLIQ')THEN
         IF(NCOMP.GT.1)CALL ERRMSG('FATAL','REAPAR',
     &                              'H2OICE/LIQ with NCOMP>1?')
      ELSE
         CALL ERRMSG('FATAL','REAPAR','Unrecognized option for CDIEL')
      ENDIF
C***********************************************************************
C
C    Define INITIALIZATION:
C
C    INIT = 0 to start with |X0>=0
C           1 to obtain |X0> from 4th-order expansion in polarizability
C           2 to read |X0> from file solvp.in
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
      READ(IOPAR,FMT=*,ERR=40)INIT
      IF(INIT.EQ.0)THEN
         CALL WRIMSG('REAPAR',
     &               'INIT=0 to start with |X0>=0 (CCG method)')
      ELSEIF(INIT.EQ.1)THEN
         CALL WRIMSG('REAPAR',
     &         'INIT=1 to obtain |X0> from 4th-order expansion (CCG)')
      ELSEIF(INIT.EQ.2)THEN
         CALL WRIMSG('REAPAR',
     &        'INIT=2 to read |X0> from file solvp.in')
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' Wrong value of INIT')
      ENDIF
C***********************************************************************
C
C    Define error tolerance:
C
C    TOL= maximum acceptable value of |Ax-E|/|E|
C
      READ(IOPAR,FMT=*,ERR=40)TOL
      IF(TOL.GT.1.E-15.AND.TOL.LT.1.E-1)THEN
         WRITE(CMSGNM,FMT='(1P,E10.3,A)')TOL,
     &         ' = TOL = max. acceptable normalized residue |Ax-E|/|E|'
         CALL WRIMSG('REAPAR',CMSGNM)
      ELSE
         CALL ERRMSG('WARNING','REAPAR',' Strange value of TOL ')
      ENDIF
C***********************************************************************
C
C    Define NUMBER OF THETA AND PHI VALUES FOR EVALUATION OF CSCA AND G
C
C    ICTHM = number of thetan values to use in evaluation of CSCA, G
C            (must be odd; recommend ICTHM=33)
C    IPHM  = number of phin values to use in evaluation of CSCA, G
C            (even or odd; recommend IPHM=12 or more)
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
      READ(IOPAR,FMT=*)ICTHM
      IF((ICTHM+1)/2*2 .NE. (ICTHM+1))THEN
         CALL ERRMSG('FATAL','REAPAR',' ICTHM SHOULD BE ODD ')
      ELSE
         WRITE(CMSGNM,FMT='(I4,A)')ICTHM,
     &' - ICTHM; number of THETAN values for calc.of g and Qsca by SCAT'
         CALL WRIMSG('REAPAR',CMSGNM)
      ENDIF
      READ(IOPAR,FMT=*)IPHM
      WRITE(CMSGNM,FMT='(I4,A)')IPHM,
     &' - IPHM; number of PHIN values for calc. of g and Qsca by SCAT'
      CALL WRIMSG('REAPAR',CMSGNM)
C***********************************************************************
C
C    Define WAVELENGTH :
C
C    WAVINI = first wavelength (physical units)
C    WAVEND = last wavelength (physical units)
C    NWAV   = number of wavelengths
C    CDIVID = 'LIN', 'LOG', or 'INV' for equal increments in
C             wavelength, log(wavelength), or frequency
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
      READ(IOPAR,FMT=*,ERR=40)WAVINI,WAVEND,NWAV,CDIVID
      WRITE(CMSGNM,FMT='(I3,A,F7.4,A,F7.4)')
     &      NWAV,' wavelengths from ',WAVINI,' to ',WAVEND
      CALL WRIMSG('REAPAR',CMSGNM)
      CALL DIVIDE(CDIVID,WAVINI,WAVEND,NWAV,MXWAV,WAVEA)
C***********************************************************************
C
C    Define AEFF = radius of equal-volume sphere (physical units)
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
      READ(IOPAR,FMT=*)AEFINI,AEFEND,NRAD,CDIVID
      WRITE(CMSGNM,FMT='(I3,A,F7.4,A,F7.4)')
     &      NRAD,' eff. radii from ',AEFINI,' to ',AEFEND
      CALL WRIMSG('REAPAR',CMSGNM)
      CALL DIVIDE(CDIVID,AEFINI,AEFEND,NRAD,MXRAD,AEFFA)
C***********************************************************************
C
C    Define incident polarizations (in Lab Frame)
C    It is assumed that incident radiation is along x axis in Lab Frame
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
C Read complex polarization vector CXE01=e01 (normalize if necessary).
      READ(IOPAR,FMT=*,ERR=40)CXE01
      E1=SQRT(CXE01(1)*CONJG(CXE01(1)))
      IF(E1.NE.0.)THEN
         CALL ERRMSG('FATAL','REAPAR',' CXE01(1) must be zero!')
      ENDIF
C Normalize:
      E1=SQRT(CXE01(2)*CONJG(CXE01(2))+CXE01(3)*CONJG(CXE01(3)))
      CXE01(2)=CXE01(2)/E1
      CXE01(3)=CXE01(3)/E1
C Construct orthogonal normalized polarization vector CXE02=e02
C using xhat cross e01* = e02
C 
      CXE02(1)=(0.,0.)
      CXE02(2)=-CONJG(CXE01(3))
      CXE02(3)=CONJG(CXE01(2))
C
C IORTH=1 to calculate only for single polarization
C      =2 to also calculate for orthogonal polarization
      READ(IOPAR,FMT=*,ERR=40)IORTH
C
C***********************************************************************
C
C    Specify whether or not to write ".sca" files
C
C IWRKSC=0 to NOT write ".sca" file for each target orientation
C       =1 to write ".sca" file for each target orientation
      READ(IOPAR,FMT=*,ERR=40)IWRKSC
C
C***********************************************************************
C
C    Read information determining target rotations
C
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
      READ(IOPAR,FMT=*)BETAMI,BETAMX,NBETA
      READ(IOPAR,FMT=*)THETMI,THETMX,NTHETA
      READ(IOPAR,FMT=*)PHIMIN,PHIMAX,NPHI
C
C*** Check that values of NBETA,NTHETA,NPHI are OK:
      IF(NBETA.GT.MXBETA)
     &   CALL ERRMSG('FATAL','REAPAR',' NBETA>MXBETA: must recompile')
      IF(NPHI.GT.MXPHI)
     &   CALL ERRMSG('FATAL','REAPAR',' NPHI>MXPHI: must recompile')
      IF(NTHETA.GT.MXTHET)
     &   CALL ERRMSG('FATAL','REAPAR',' NTHETA>MXTHET: must recompile')
C
C If NPHI>1, then set IORTH=2 regardless of value input.
      IF(IORTH.EQ.1.AND.NPHI.GT.1)THEN
         IORTH=2
         CALL WRIMSG('REAPAR',
     &               'Set IORTH=2 since NPHI>1 ')
      ENDIF
C
      IF(IORTH.EQ.1)THEN
         CALL WRIMSG('REAPAR',
     &               'Calculate only for single polarization ')
      ELSEIF(IORTH.EQ.2)THEN
         CALL WRIMSG('REAPAR',
     &               'Do orthogonal polarization for each case ')
      ELSE
         CALL ERRMSG('FATAL','REAPAR',' WRONG VALUE OF IORTH ')
      ENDIF
      WRITE(CMSGNM,FMT='(2F7.2,A,I4)')BETAMI,BETAMX,
     &  ' Range of BETA values ; NBETA =',NBETA
      CALL WRIMSG('REAPAR',CMSGNM)
      WRITE(CMSGNM,FMT='(2F7.2,A,I4)')THETMI,THETMX,
     &  ' Range of THETA values; NTHETA=',NTHETA
      CALL WRIMSG('REAPAR',CMSGNM)
      WRITE(CMSGNM,FMT='(2F7.2,A,I4)')PHIMIN,PHIMAX,
     &  ' Range of PHI values ;   NPHI =',NPHI
      CALL WRIMSG('REAPAR',CMSGNM)
C Convert from degrees to radians
      BETAMI=BETAMI/DEGRAD
      BETAMX=BETAMX/DEGRAD
      THETMI=THETMI/DEGRAD
      THETMX=THETMX/DEGRAD
      PHIMIN=PHIMIN/DEGRAD
      PHIMAX=PHIMAX/DEGRAD
C***********************************************************************
C
C    Specify scattering directions to be calculated
C    THETAN,PHIN = Direction of n from n0 in calculation of the
C                  scattering matrix; cos(THETAN) is n0 \dot n
C                  and PHIN is azimuthal angle of n from Lab xy plane
C                  Note: THETA1, THETA2, and PHI for each scattering
C                  plane are entered in degrees, and immediately 
C                  converted to radians.  Arbitrary number
C                  of scattering planes may be considered.
C
      NSCA=0
      READ(IOPAR,FMT=*,ERR=40)CLINE
      CALL WRIMSG(' ',CLINE)
 22   READ(IOPAR,FMT=*,ERR=40,END=25)PHI1,THETA1,THETA2,DTHETA
      WRITE(CMSGNM,FMT='(3F7.1,A)')PHI1,THETA1,THETA2,
     &  ' = phi, theta_min, theta_max for scattering plane'
      CALL WRIMSG('REAPAR',CMSGNM)
C Convert to radians:
      PHI1=PHI1/DEGRAD
      THETA1=THETA1/DEGRAD
      THETA2=THETA2/DEGRAD
      DTHETA=DTHETA/DEGRAD
C Allow for possibility that user did not enter correct sign for
C DTHETA (if I did it, others will too...)
      IF(THETA2.LT.THETA1)DTHETA=-ABS(DTHETA)
C compute theta values for this scattering plane
      NSCA0=1
      DELTA=THETA2-THETA1
      IF(DELTA.NE.0)NSCA0=1+INT(DELTA/DTHETA+0.5)
C Check that there is enough storage space:
      IF((NSCA+NSCA0).GT.MXSCA)THEN
         WRITE(IDVOUT,FMT='(I6,I6,A)')MXSCA,NSCA0,
     &         ' =MXSCA,NSCA0 (dirs in this scatt. plane)'
         CALL ERRMSG('FATAL','REAPAR',
     &               'MXSCA too small in DDSCAT')
      ENDIF
      THETAN(NSCA+1)=THETA1
      IF(NSCA0.GT.2)THEN
         DO 23 J=2,NSCA0-1
            THETAN(NSCA+J)=THETA1+(J-1)*DTHETA
 23      CONTINUE
      ENDIF
      THETAN(NSCA+NSCA0)=THETA2
      DO 24 J=1,NSCA0
         PHIN(NSCA+J)=PHI1
 24   CONTINUE
      NSCA=NSCA+NSCA0
      WRITE(CMSGNM,FMT='(I4,A)')NSCA0,
     &  ' = number of scattering angles in this scattering plane'
      CALL WRIMSG('REAPAR',CMSGNM)
      GOTO 22
 25   CONTINUE
C
C Check that if user has requested computation of f_ml, user has
C also requested linearly-polarized incident radiation.
C 96.01.25 (BTD) We have disabled this "protection".  So long as user
C is aware that when using elliptical incident polarization the
C f_ml are not the "usual" ones with index l referring to linearly
C polarized incident modes, there is not reason not to allow user
C to calculate f_ml for the the scattered radiation.
C
C      IF(NSCA.GT.0)THEN
C         DO 26 J=1,3
C            IF(AIMAG(CXE01(J)).NE.0.)CALL ERRMSG('FATAL','REAPAR',
C     &      'pol. vector E01 is complex: DDSCAT cannot compute f_ml, so 
C     &do not ask it to!')
C 26      CONTINUE
C      ENDIF
C
      CLOSE (IOPAR)
      RETURN
   40 CONTINUE
      CALL ERRMSG('FATAL','REAPAR',' Error reading ddscat.par file')
      RETURN
      END
      SUBROUTINE REASHP(CFLSHP,IDVOUT,IDVSHP,A1,A2,CDESCR,MXNAT,NAT,
     &                  IXYZ,ICOMP,MXN3,NAT3)
C Arguments:
      INTEGER IDVOUT,IDVSHP,MXNAT,NAT,MXN3,NAT3
      CHARACTER CDESCR*67,CFLSHP*(*)
      INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3)
      REAL A1(3),A2(3)
C Local variables:
      INTEGER J,JJ
C***********************************************************************
C Given:
C    CFLSHP = name of (unformatted) file containing target description
C    IDVSHP  = device number to use in reading this file
C    IDVOUT = device number to use for "running" output
C    MXNAT,MXN3 = dimensioning information
C Returns:
C    A1(3)  = target axis 1
C    A2(3)  = target axis 2
C    CDESCR = string describing target
C    NAT    = number of dipoles in target
C    NAT3   = 3*NAT
C    IXYZ(1-NAT,1-3) = lattice coordinates of dipoles
C    ICOMP(1-NAT,1-3)= composition for directions x,y,z (in Target
C                      Frame) of dipoles 1-NAT.
C
C Note: It is expected that many users will wish to modify this routine
C       to be consistent with whatever target generation software they
C       use.  However, as an example (and for test purposes) this routine
C       is supplied with code compatible with the "target.out" files which
C       are generated by the "calltarget" program when called for
C       homogenous isotropic targets (shape.dat files do not include
C       information on target composition).
C       Also provided below (commented out) are
C       * a module of code appropriate for reading formatted target.out
C         file generated by the routine tarcel.f (concentric ellipsoids)
C         as an example of an inhomogeneous target.
C       * a module of code appropriate for reading an unformatted 
C         "shape.dat" file for a general target, assuming user's
C         target generation routine writes an unformatted file with
C         the same structure.
C
C
C***********************************************************************
C History:
C 91.05.02 (BTD): Added IDVOUT to argument list.
C 91.05.23 (BTD): Added A1, A2 to argument list.
C 91.11.11 (BTD): Modified for compatibility with WRITE statements in
C                 program CALLTARGET; Changed dimensioning of ICOMP from
C                 ICOMP(MXN3) to ICOMP(MXNAT,3)
C 93.09.24 (BTD): Modified to include examples of READ statements
C                 for both unformatted and formatted shape.dat files.
C                 Leave the formatted example (compatible with the
C                 "target.out" files generated by "calltarget") in effect.
C 95.04.07 (BTD): Changed CDESCR*60 -> CDESCR*67 for compatibility with
C                 calling routine TARGET
C 95.05.15 (BTD): Changed  READ(DVSHP,FMT='(A60)') to ..(A67)..
C                 Added module for inhomogeneous/anisotropic targets
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***********************************************************************
C
C Enable following module of code if "shape.dat" file has been written
C in formatted form consistent with the "target.out" files generated by
C routines tar2el.f, tar3el.f, tarcyl.f, tarell.f, tarhex.f, tarrec.f,
C and tartet.f (in particular, this assumes homogeneous target):
C
      OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='FORMATTED')
      READ(IDVSHP,FMT='(A67)')CDESCR
      READ(IDVSHP,*)NAT
      READ(IDVSHP,*)A1
      READ(IDVSHP,*)A2
      READ(IDVSHP,*)
      READ(IDVSHP,*)(JJ,IXYZ(J,1),IXYZ(J,2),IXYZ(J,3),J=1,NAT)
      DO 0200 JJ=1,3
         DO 0100 J=1,NAT
            ICOMP(J,JJ)=1
 0100    CONTINUE
 0200 CONTINUE
      CLOSE(IDVSHP)
C
C***********************************************************************
C
C Enable following module of code if "shape.dat" file has been written
C in formatted form consistent with the "target.out" files generated by
C routine tarcel.f (inhomogenous target):
C
C      OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='FORMATTED')
C      READ(IDVSHP,FMT='(A67)')CDESCR
C      READ(IDVSHP,*)NAT
C      READ(IDVSHP,*)A1
C      READ(IDVSHP,*)A2
C      READ(IDVSHP,*)
C      READ(IDVSHP,*)(JJ,IXYZ(J,1),IXYZ(J,2),IXYZ(J,3),
C     &               ICOMP(J,1),ICOMP(J,2),ICOMP(J,3),J=1,NAT)
C      CLOSE(IDVSHP)
C
C***********************************************************************
C
C Enable following module of code if "shape.dat" file has been written
C in unformatted form consistent with these read statements:
C
C      INTEGER JX,JY,JZ
C      OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='UNFORMATTED')
C      READ(IDVSHP)CDESCR,NAT,A1,A2
C      READ(IDVSHP)((IXYZ(JX,JY),JX=1,NAT),JY=1,3),
C     &            ((ICOMP(J,JJ),J=1,NAT),JJ=1,3)
C      CLOSE(IDVSHP)
C
C***********************************************************************
      NAT3=3*NAT
      RETURN
      END

      SUBROUTINE REDUCE(CXV,IOCC,MXN3,MXNAT,NAT,NAT0)
C***********************************************************************
C Given:
C      CXV(1-3*NAT) defined for NAT lattice sites
C      IOCC(1-NAT) = 0 or 1 depending on whether lattice site
C                      is vacant or occupied
C      MXN3,MXNAT = dimensioning information
C      NAT = number of lattice sites
C      NAT0 = number of occupied lattice sites
C Returns
C      CXV(1-3*NAT0) defined for NAT0 occupied lattice sites
C
C Arguments:
      INTEGER MXN3,MXNAT,NAT,NAT0
      INTEGER*2 IOCC(MXNAT)
      COMPLEX CXV(MXN3)
C Local variables:
      INTEGER J,JOC,M
C
C B.T.Draine, Princeton Univ. Observatory, 90/11/29
C History:
C 90/12/03 (BTD): Modified for new ordering of vectors.
C 90/12/05 (BTD): Corrected errors
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 M=0,2
         JOC=0
         DO 10 J=1,NAT
            IF(IOCC(J).EQ.1)THEN
               JOC=JOC+1
               CXV(JOC+M*NAT0)=CXV(J+M*NAT)
            ENDIF
  10     CONTINUE
  20  CONTINUE
      RETURN
      END
