PROGRAM DDSCAT C********************************************************************** C C DDSCAT is a program to use the Discrete Dipole Approximation (DDA) C to calculate the absorption and scattering properties of targets C of arbitrary geometry and dielectric properties. The only C approximation is that the target may be approximated by an array C of point dipoles on a cubic lattice. C DDA theory and validity criteria are presented by B. T. Draine C (1988; Astrophys. J. 333, pp. 848-872). C The original version of DDSCAT was developed by B. T. Draine, C Princeton University Observatory. C C The current version of DDSCAT (Version 4c) was jointly developed by C B. T. Draine (Princeton University Observatory) and C P. J. Flatau (University of California, San Diego, Scripps Inst. C of Oceanography, La Jolla) C C DDSCAT (Version 4c) includes the following features: C 1. Optional generation of target arrays for to approximate several C standard target shapes, including C a. ellipsoid (sphere is special case) C b. cylindrical prism C c. rectangular prism C d. hexagonal prism C e. regular tetrahedron C f. cylindrical prism with anisotropic diel. tensor ("UNICYL") C (uniaxial material with crystal axis = cylinder axis) C g. ellipsoid with anisotropic diel. tensor ("ANIELL") C h. two touching ellipsoids with either isotropic or aniostropic C dielectric tensors ("TWOELL" or "TWOAEL") C i. three thouching ellipsoids with either isotropic or aniostropic C dielectric constants ("THRELL" or "THRAEL") C j. two concentric ellipsoids ("CONELL") C 2. Alternatively, program may read in target array properties from C ascii file. C 3. User may easily specify target orientation with respect to C incident plane wave. C 4. Solution is found iteratively using Complex Conjugate Gradient C methods. C 5. User may select FFT implementation: C a. Brenner's implementation; C b. Temperton's implementation; C c. Convex vector library implementation (if on Convex); C d. GPFA code of Temperton (recommended for general use). C 6. User may select prescription for determining dipole C polarizabilities: C a. Lattice Dispersion Relation ("LATTDR"; this is recommended) C b. Isotropized Lattice Dispersion Relation ("LDRISO"); C c. Prescription of Goedecke & O'Brien (1988) ("GOBR88"); C d. Clausius-Mossotti with Radiative Reaction, as used by C Draine (1988) ("DRAI88"). C 7. Automatic computation of absorption, extinction, and scattering C cross sections. C 8. User may easily specify scattering directions for which C scattering matrix is to be computed. C 9. Automatic computation of orientational averages, with C orientation weights provided by simple subroutine. C 10. Multicomponent targets are allowed. C 11. Anisotropic dielectric tensors are allowed, although with C the restriction that the principal axes must be aligned with C the xyz directions in the "target frame". C 12. Wavelength-dependent dielectric functions are read from user- C provided files, with automatic interpolation if necessary. C 13. Automatic looping over radii and wavelengths if desired. C C History: C Note: Overall history of significant changes to the DDSCAT package may C be found in the comments in subroutine version.f . C Here we limit comments to changes to this program module, DDSCAT.f: C 90.11.30 (BTD): Removed LSC3 portion of CXSC complex scratch space. C 90.12.03 (BTD): Reordered vectors. C Removed LSC2 and LSC1 portions of CXSC scratch space. C 90.12.21 (BTD): Added code to create output files C 'qtable' (summary of Q values) C 'mtable' (summary of diel. func. for material 1) C 91.01.03 (BTD): Added MXBETA,MXPHI,MXTHET to argument list for REAPAR C 91.01.05 (BTD): Changed I4 -> I7 in format statements 9010,9011. C Set IOSHP=-1 to suppress printing of target.out file. C 91.05.08 (BTD): Added CALPHA to argument lists for REAPAR and ALPHA C and added CALPHA to various WRITE statements. C 91.05.23 (BTD): Added IWRKSC to arg. list for REAPAR and use it to C control creation of wAArBBkCCC.sca files C Move lines calculating G. C Move code writing to qtable and mtable. C Added IDVOUT to argument list for TARGET C 91.08.14 (BTD): Add QBKSCA to arg. list for GETFML C Modify code to print out Q_bk C 91.08.15 (BTD): Provide different headings for qtable depending on C whether IORTH=1 or 2 (add statement 9043). C 91.09.17 (BTD): Remove calls to ALPHA and EVALA (needed to move these C to subroutine GETFML since alpha from Lattice C Dispersion Relation depends on propagation direction C and polarization state) C Add CALPHA,CXEPS,ICOMP,MXCOMP to argument list for C GETFML (information required by ALPHA) C 91.11.12 (BTD): Added variable IDVSHP to argument list for TARGET C (device number for REASHP to use in reading shape file) C 93.01.15 (BTD): Divided output file qtable into 2 files: C qtable (containing Q_ext,Q_abs,Q_sca,g,Q_bk) and C qtable2 (containing Q_pha, Q_pol, Q_cpol) C 93.01.16 (BTD): Modify so that qtable, qtable2, and mtable are closed C after writing, and reopened for each new write C 93.01.20 (BTD): Add MXNX, MXNY, MXNZ to argument list for subroutine C EXTEND to permit checking of target size against C maximum allowed dimensions C 93.03.11 (BTD): Deleted all code associated with unused variable C CMACHN (originally included to identify machine/OS) C 93.03.12 (BTD): Changed CDESCR*60 -> CDESCR*67 C Moved WRITE(IDVOUT,9010) to subr. TARGET C 93.06.02 (BTD): Corrected error in FORMAT statement 9045 C 93.06.03 (BTD): Corrected error in computation of angle-averaged C backscattering cross section QBKSUM(JO) (had neglected C to reset sum to zero for each new target). C 93.09.28 (BTD): Add "fix" to work around Sun compiler/OS bug (see C description below. C 93.12.15 (BTD): Corrected calculation of PPOL. Added comments on C correspondence between our notation and elements C of 2x2 amplitude scattering matrix and 4x4 Mueller C scattering matrix. C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) to C allow up to 6 shape parameters C 94.06.20 (PJF): Comment about NEWTMP FFT and corresponding changes C in reapar, eself, extend. Add FFT timing code C TSTFFT to distribution and changes to cxfft3 and C timeit routines. C 94.12.20 (BTD): Introduce subroutine VERSION to set value of C variable CSTAMP. version.f will now serve as C location to maintain log of significant code C changes. C C Copyright (C) 1993,1994 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C Adjustable Parameters: C C MXNX = max. extent of target in x direction C MXNY = y C MXNZ = z C MXMEM = 0 to use Brenner FFT only, 1 to allow use of Temperton FFT C MXCOMP = max. number of different dielectric functions C MXTHET = max. number of target rotation angles THETA C MXBETA = max. number of target rotation angles BETA C MXPHI = max. number of target rotation angles PHI C MXSCA = max. number of scattered directions C MXRAD = max. number of radii C MXWAV = max. number of wavelengths C MXWAVT = max. number of wavelengths in dielectric tables C C*** Set parameter MXMEM *********************************************** INTEGER MXMEM C PARAMETER(MXMEM=0) PARAMETER(MXMEM=1) C C*** Set parameters MXNX,MXNY,MXNZ ************************************* INTEGER MXNX,MXNY,MXNZ PARAMETER(MXNX=8,MXNY=8,MXNZ=8) C PARAMETER(MXNX=16,MXNY=16,MXNZ=16) C PARAMETER(MXNX=12,MXNY=18,MXNZ=18) C PARAMETER(MXNX=16,MXNY=24,MXNZ=24) C PARAMETER(MXNX=20,MXNY=20,MXNZ=20) C PARAMETER(MXNX=25,MXNY=25,MXNZ=25) C PARAMETER(MXNX=32,MXNY=32,MXNZ=32) C PARAMETER(MXNX=24,MXNY=48,MXNZ=42) C PARAMETER(MXNX=40,MXNY=40,MXNZ=40) C PARAMETER(MXNX=48,MXNY=16,MXNZ=16) C PARAMETER(MXNX=48,MXNY=48,MXNZ=48) C PARAMETER(MXNX=96,MXNY=32,MXNZ=32) C C*** Set parameters MXCOMP,MXTHET,MXBETA,MXPHI,MXRAD,MXWAV,MXWAVT,MXSCA INTEGER MXCOMP PARAMETER(MXCOMP=9) INTEGER MXTHET,MXBETA,MXPHI PARAMETER(MXTHET=10,MXBETA=10,MXPHI=100) INTEGER MXRAD,MXWAV,MXWAVT PARAMETER(MXRAD=100,MXWAV=100,MXWAVT=500) INTEGER MXSCA PARAMETER(MXSCA=1000) C C*** Derived Parameters: ********************************************** INTEGER MXNAT,MXN3 PARAMETER(MXNAT=MXNX*MXNY*MXNZ,MXN3=3*MXNAT) INTEGER MXCXSC C MXCXSC specifies required complex scratch array CXSC PARAMETER(MXCXSC=7*MXN3) C C Local variables: C COMPLEX CXEN,CXVAR1 C REAL & AEFF,AK1,AK3,BETAD,BETAMI,BETAMX,BETMID,BETMXD, & DEGRAD,DENOM,DEPOLR,FREQ, & PHID,PHIMAX,PHIMID,PHIMIN,PHIMXD,PHIND,PI,PIA2,PPOL, & RVAR1,RVAR2,RVAR3,RVAR4,SUMWGT, & THETAD,THETMI,THETMX,THETND,THTMID,THTMXD,WAVE,TOL, & WG,XX C C Target Properties: C AEFF=aeff=effective radius of target (physical units) C PIA2=pi*(aeff/d)**2=pi*(3*NAT/4*pi)**(2/3) C C Incident Wave Properties: C AK1=(k*d), where k=2*pi/wave=propagation vector and d=dipole spacing C AK3=(k*d)**3 C AKE2=(k*d \dot E0), where E0 is incident polarization vector C WAVE=wavelength (in vacuo) in physical units C XX=k*aeff=2*pi*aeff/wave C C Scattering Properties: C G= for scattered radiation C GSUM(1-2)=weighted sum of G*QSCA over target orientations, for each C incident polarization C QPHA(1-2)=(phaseshift cross section)/(pi*aeff**2) C QEXSUM(1-2)=running sum of QEXT over target orientations C QABSUM(1-2)=running sum of QABS over target orientations C QBKSUM(1-2)=running sum of QBKSCA over target orientations C QSCSUM(1-2)=running sum of QSCA over target orientations C THETND=scattering angle theta in degrees C PHIND=scattering angle phi in degrees C SINPHI=sin(phi) C DEPOLR=depolarization ratio for scattered of initially polarized C light into direction theta,phi C PPOL=percentage polarization for scattering of initially unpolarized C light into direction theta,phi C INTEGER & IBETA,ICTHM,IDIR,IDIRWR,IDVERR,IDVOUT,INIT,IOPAR,IORTH, & IDVSHP,ILIN10,ILIN12,IOSHP,IPHI,IPHIM,IRAD,ITASK,ITHETA, & IWAV,IWRKSC, & J,JO, & LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0, & NAT,NAT0,NAT3,NBETA,NCOMP,ND,NORI,NPHI,NRAD,NSCAT,NTHETA, & NWAV,NX,NY,NZ C C NAT=number of dipoles in target C NAT3=3*NAT C NBETA=number of different beta values for target orientation C NTHETA=number of different theta values for target orientation C NPHI=number of different phi values for target orientation C NORI=NBETA*NTHETA*NPHI=total number of target orientations C NSCAT=number of different scattering directions for each orientation C NWAV=number of different wavelengths C NRAD=number of different target radii C INIT=0,1,2 for choice of initial vector |x0> in complex conjugate gradient C method C IORTH=1 or 2 to do 1 or 2 incident polarizations C IWRKSC=0 or 1 to suppress or generate "wAArBBkCCC.sca" files for each C target orientation C NRXY=total number of elements of CXAOFF C NX=x-length of extended target C NY=y-length of extended target C NZ=z-length of extended target C CHARACTER & CALPHA*6,CDIEL*6,CMETHD*6,CSHAPE*6,CSTAMP*11, & CFLAVG*13,CFLPAR*13,CFLSCA*14,CFLSHP*13,CDESCR*67 CHARACTER*40 CFLEPS(MXCOMP) C .. C .. Local Arrays .. COMPLEX & CX1112(MXSCA),CX1121(MXSCA),CX1122(MXSCA), & CX1221(MXSCA),CX1222(MXSCA),CX2122(MXSCA), & CXADIA(MXN3),CXALPH(MXN3), & CXE(MXN3),CXE01(3),CXE02(3), & CXE01R(3),CXE02R(3),CXEPS(MXCOMP), & CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),CXF22(MXSCA), & CXRFR(MXCOMP),CXSC(MXCXSC), & CXXI(MXN3), & CXZC(MXNX+1,MXNY+1,MXNZ+1,6), & CXZW(2*MXNX,2*MXNY,2*MXNZ,3+MXMEM) C C Arrays describing target properties: C CXADIA(1-3*NAT)=diagonal elements of A matrix C CXEPS(1-MXCOMP)=dielectric const. for each of MXCOMP materials C CXREFR(1-3)=complex refractive index at current wavelength C CXALPH(1-3*NAT)=diagonal elements of dipole polarizability tensor alpha C C Arrays describing incident radiation: C CXE01(1-3)=incident polarization vector e01 prior to target rotation C CXE02(1-3)=incident polarization vector e02 prior to target rotation C CXE01R(1-3)=rotated incident polarization vector e01 in Target Frame C CXE02R(1-3)=rotated incident polarization vector e02 in Target Frame C CXE0R(1-3)=incident polarization vector (=CXE01R or CXE02R) being used C CXE(1-3*NAT)=incident electric field vector in Target Frame C C Array describing polarization of target for one incident wave: C CXXI(1-3*NAT)=polarization vector at each dipole C C Arrays describing scattering properties of target: C CXFF11(1-NSCAT)=scattering matrix element f11 for each of NSCAT directions C CXFF12(1-NSCAT)= f12 C CXFF21(1-NSCAT)= f21 C CXFF22(1-NSCAT)= f22 C CX1112(1-NSCAT)=sum of f11* \times f12 over orientations C CX1121(1-NSCAT)= f11* \times f21 C CX1122(1-NSCAT)= f11* \times f22 C CX1221(1-NSCAT)= f12* \times f21 C CX1222(1-NSCAT)= f12* \times f22 C CX2122(1-NSCAT)= f21* \times f22 C C Scratch arrays: C CXSC(1-MXCXSC)=scratch array used by DDACCG C INTEGER*2 IOCC(MXNAT),ICOMP(MXN3),ISCR1(MXN3),IXYZ0(MXNAT,3) C C IX0(1-NAT0)=x/d for dipoles 1-NAT in real target C IY0(1-NAT0)=y/d C IZ0(1-NAT0)=z/d C C Scratch arrays: C ISCR1,ISCR2,ISCR3,ISCR4 (used by EXTEND) C C REAL & A1(3),A2(3),AEFFA(MXRAD),AKSR(3,MXSCA),AKR(3), & BETA(MXBETA), & E1A(MXWAVT),E2A(MXWAVT),EM1(3,MXSCA),EM1R(3,MXSCA), & EM2(3,MXSCA),EM2R(3,MXSCA),EN0R(3),ENSC(3,MXSCA), & G(2),GSUM(2),PHI(MXPHI),PHIN(MXSCA), & QABS(2),QABSUM(2),QAV(8),QBKSCA(2),QBKSUM(2),QEXT(2), & QEXSUM(2),QPHA(2),QPHSUM(2),QSCAT(2),QSCATG(2),QSCSUM(2), & S1111(MXSCA),S1212(MXSCA),S2121(MXSCA),S2222(MXSCA), & SCRRS1(MXN3),SHPAR(6),THETA(MXTHET),THETAN(MXSCA), & WAVEA(MXWAV),WGTA(MXTHET,MXPHI),WGTB(MXBETA),WVA(MXWAVT) C C Arrays describing target properties: C A1(1-3)=initial direction of target axis a1 in Lab Frame C A2(1-3)=initial direction of target axis a2 in Lab Frame C AEFFA(1-NRAD)=values of effective radius of target (physical units) C SHPAR(1-6)=up to 6 parameters for defining target geometry C C Arrays describing target rotation: C THETA(1-NTHETA)=desired values of target rotation angle theta C BETA(1-NBETA)=desired values of target rotation angle beta C PHI(1-NPHI)=desired values of target rotation angle phi C WGTA(1-NTHETA,1-NPHI)=weight factor for each orientation of target C axis a1 in Lab Frame C WGTB(1-NBETA)=weight factor for rotation of target around a1 C Arrays describing incident wave: C WAVEA(1-NWAV)=values of wavelength (physical units) C EN0R(3)=unit vector in direction of incidence in (rotated) Target C Frame C AKR(3)=k*d vector in (rotated) Target Frame C C Arrays describing scattering properties: C ENSC(1-3,1-NSCAT)=scattering vectors in Lab Frame C AKSR(1-3,1-NSCAT)=scattering k-vectors in Target Frame C EM1(1-3,1-NSCAT)=scattering polarization state 1 in Lab Frame C EM2(1-3,1-NSCAT)=scattering polarization state 2 in Lab Frame C EM1R(1-3,1-NSCAT)=scattering polarization state 1 in Target Frame C EM1R(1-3,1-NSCAT)=scattering polarization state 2 in Target Frame C THETAN(1-NSCAT)=values of theta for NSCAT scattering directions nhat C at which scattering matrix fml is to be calculated C PHIN(1-NSCAT)=values of phi for scattering NSCAT directions nhat at C which scattering matrix fml is to be calculated C S1111(1-NSCAT)=sum of |f11|**2 over target orientations C S1212(1-NSCAT)= |f12|**2 C S2121(1-NSCAT)= |f21|**2 C S2222(1-NSCAT)= |f22|**2 C QABS(1-2)=(absorption cross section)/(pi*aeff**2) for 2 polarizations C QBKSCA(1-2)=4*pi*(diff.scatt.cross section for theta=pi)/(pi*aeff**2) C for 2 polarizations C QEXT(1-2)=(extinction cross section)/(pi*aeff**2) C QSCA(1-2)=(scattering cross section)/(pi*aeff**2)=QEXT-QABS C QSCAT(1-2)=(scattering cross section)/(pi*aeff**2) estimated by SCAT C C Scratch array: C SCRRS1(1-3*NAT)=scratch array used by DDACCG and SCAT C C External subroutines: EXTERNAL DIELEC,EVALA,GETSET,NAMER,REAPAR C C Intrinsic functions: INTRINSIC ATAN,SQRT,CONJG,REAL C C Set device numbers: C IDVOUT = device number for running I/O C on Sun4OS, IDVOUT=0 produces unbuffered output to standard C output C IDVOUT=6 produces buffered output to standard C output C IDVERR = device number for error messages C IOPAR = device number for reading ddscat.par file C IOSHP = device number for writing target.out (IOSHP<0 to suppress) C IDVSHP = device number for reading shape file DATA IDVOUT/0/ DATA IDVERR/0/ DATA IDVSHP/10/ DATA IOPAR/1/ DATA IOSHP/-1/ C VMS version: C DATA CFLPAR,CFLSHP/'DDSCAT.PAR','SHAPE.DAT'/ C unix version: DATA CFLPAR,CFLSHP/'ddscat.par','shape.dat'/ C*** BTD 93.09.28: Workaround ****************************************** C Problem with Sun Fortran V1.4 (1 Mar 1991) running under SunOS4.1.1 C on Sparc IPX with tmpfs (cohabitation of /tmp and swap on same disk C partition). C Symptom: without this fix, the value of IDVOUT was "stochastic" C from one compilation of DDSCAT.f to another, instead of being C initialized to 0 by DATA statement above. It has not been possible C to understand why/where this problem is occuring, but it appears to C be either a compiler or operating system bug. Work around this by C adding line explicitly (re-)initializing IDVOUT: IDVOUT=0 C*********************************************************************** CALL VERSION(CSTAMP) WRITE(IDVOUT,FMT=9000)CSTAMP,MXNX,MXNY,MXNZ,MXMEM C .. C Pseudo names of scratch array positions: LACE=1 LAXI=1*MXN3+1 LCLM=2*MXN3+1 LGI=3*MXN3+1 LPI=4*MXN3+1 LQI=5*MXN3+1 LSC0=6*MXN3+1 C PI=4.*ATAN(1.) DEGRAD=180./PI C*** temporary hack -- do we really need getset? CALL GETSET('SET','IOOUT',IDVOUT) CALL GETSET('SET','IOERR',IDVERR) C C**** Read program control parameters C CMETHD=choice of solution methods C CSHAPE=description of shape C SHPAR(1-6)= up to 6 parameters for defining target geometry C INIT=selection of |x0> in CCG method C TOL=error tolerance in CCG method C ICTHM=number of theta values for computation of CSCA and G by SCAT C IPHIM=number of phi values for computation of CSCA and G by SCAT C MXWAV=dimensioning information for array WAVEA C NWAV=number of wavelengths to be used C WAVEA(1-NWAV)=wavelengths (physical units) C MXRAD=dimensioning information for array AEFFA C NRAD=number of radii to be used C AEFFA(1-NRAD)=effective radii (physical units) C IORTH=1 or 2 for 1 or 2 incident polarization states to be used C NBETA=number of beta values to be used for target orientation C BETAMI=minimum value of beta (radians) C BETAMX=maximum value of beta (radians) C NTHETA=number of theta values to be used for target orientation C THETMI=minimum value of theta (radians) C THETMX=maximum value of theta (radians) C NPHI=number of phi values to be used for target orientation C PHIMIN=minimum value of phi (radians) C PHIMAX=maximum value of phi (radians) C DEGRAD=degrees/radian (conversion factor) C MXSCA=dimensioning information for arrays PHIN, THETAN C NSCAT=number of scattered directions for computation of fml C THETAN(1-NSCAT)=theta values (radians) for scattered directions C PHIN(1-NSCAT)=phi values (radians) for scattered directions C CALL REAPAR(AEFFA, & BETAMI,BETAMX, & CALPHA,CDIEL,CFLEPS,CFLPAR,CMETHD,CSHAPE, & CXE01,CXE02, & DEGRAD, & ICTHM,IDVERR,IDVOUT,INIT,IOPAR,IORTH,IPHIM,IWRKSC, & MXBETA,MXCOMP,MXMEM,MXPHI,MXRAD,MXSCA,MXTHET,MXWAV, & NBETA,NCOMP,NPHI,NRAD,NSCAT,NTHETA,NWAV, & PHIMIN,PHIMAX,PHIN, & THETAN,THETMI,THETMX,TOL, & WAVEA,SHPAR) C NORI=NTHETA*NBETA*NPHI C Compute angles in degrees for printout BETMID=DEGRAD*BETAMI BETMXD=DEGRAD*BETAMX THTMID=DEGRAD*THETMI THTMXD=DEGRAD*THETMX PHIMID=DEGRAD*PHIMIN PHIMXD=DEGRAD*PHIMAX C C*** Obtain scattering vectors and scattering pol. vectors in Lab Frame C CALL SCAVEC(MXSCA,NSCAT,THETAN,PHIN,CXE01,ENSC,EM1,EM2) C C*** Specify target properties C C A1(1-3)=target axis 1 in Target Frame C A2(1-3)=target axis 2 in Target Frame C CALL TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,NAT0,IXYZ0(1,1),IXYZ0(1,2), & IXYZ0(1,3),ICOMP,IDVOUT,MXN3,NAT3) C C NAT0=number of dipoles in "real" target C*** Extend target to rectangular volume suitable for FFT C CALL EXTEND(CMETHD,ICOMP,ISCR1,IDVERR,IDVOUT,IOCC, & IXYZ0(1,1),IXYZ0(1,2),IXYZ0(1,3), & MXNX,MXNY,MXNZ,MXNAT,MXN3,NAT,NAT0,NAT3,NX,NY,NZ) C C Note: IXYZ0(1-NAT0,1-3) contain coordinates of occupied lattice sites. C WRITE(IDVOUT,FMT=9011)NAT,NX,NY,NZ C C**** Specify target orientations C C MXBETA=dimensioning information for array BETA of beta values C MXTHET=dimensioning information for array THETA of theta values C MXPHI=dimensioning information for array PHI of phi values C BETAMI,BETAMX=minimum,maximum beta values (radians) C THETMI,THETMX=minimum,maximum theta values (radians) C PHIMIN,PHIMAX=minimum,maximum phi values (radians) C BETA(1-NBETA)=beta values for target orientations (radians) C THETA(1-NTHETA)=theta values for target orientations (radians) C PHI(1-NPHI)=phi values for target orientations (radians) C WGTA(1-NTHETA,1-NPHI)=weight for each orientation of axis a1 in LF C WGTB(1-NBETA)=weight for each rotation of target around a1 C CALL ORIENT(BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX, & MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI, & BETA,THETA,PHI,WGTA,WGTB) C C*** Open file 'qtable' for summary of Q values and C 'mtable' for summary of refractive indices C OPEN(UNIT=10,FILE='qtable',STATUS='UNKNOWN') OPEN(UNIT=11,FILE='mtable',STATUS='UNKNOWN') OPEN(UNIT=12,FILE='qtable2',STATUS='UNKNOWN') C C*** Write out header for 'qtable', 'qtable2', and 'mtable' files WRITE(10,FMT=9020)CSTAMP,CDESCR,CMETHD,CALPHA,CSHAPE,NAT0,CDIEL ILIN10=7 WRITE(10,FMT=9021)(CFLEPS(J),J=1,NCOMP) ILIN10=ILIN10+NCOMP WRITE(12,FMT=9020)CSTAMP,CDESCR,CMETHD,CALPHA,CSHAPE,NAT0,CDIEL ILIN12=7 WRITE(12,FMT=9021)(CFLEPS(J),J=1,NCOMP) ILIN12=ILIN12+NCOMP IF(IORTH.EQ.1)THEN WRITE(10,FMT=9043)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN10=ILIN10+6 WRITE(12,FMT=9045)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN12=ILIN12+6 ELSE WRITE(10,FMT=9044)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN10=ILIN10+6 WRITE(12,FMT=9046)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH ILIN12=ILIN12+6 ENDIF CLOSE(10) CLOSE(12) WRITE(11,FMT=9020)CSTAMP,CDESCR,CMETHD,CALPHA,CSHAPE,NAT0,CDIEL WRITE(11,FMT=9021)(CFLEPS(J),J=1,NCOMP) WRITE(11,FMT=9048) C C*** Loop over wavelengths: C DO 150 IWAV=1,NWAV WAVE=WAVEA(IWAV) C C*** Determine dielectric properties of target C CALL DIELEC(CDIEL,WAVE,IDVOUT,CFLEPS,CXEPS, & MXCOMP,MXWAVT,NCOMP,E1A,E2A,WVA) C C*** Obtain complex polarizabilities C DO 20 J=1,NCOMP CXRFR(J)=SQRT(CXEPS(J)) WRITE(IDVOUT,FMT=9031)CXRFR(J),CXEPS(J),J 20 CONTINUE C C**** Loop over grain sizes: C DO 140 IRAD=1,NRAD AEFF=AEFFA(IRAD) XX=2.*PI*AEFF/WAVE WRITE(IDVOUT,FMT=9032)AEFF,WAVE,XX C C Compute AK1=length of k vector in "natural" units=k*d C (natural unit of length = d = lattice spacing) C Remember that NAT = number of dipoles including "vacuum" sites in C extended target C NAT0= number of dipoles in "real" target AK1=XX*(4.*PI/(3.*REAL(NAT0)))**(1./3.) AK3=AK1**3 C PIA2=pi*(aeff/d)**2 PIA2=PI*(.75*NAT0/PI)**(2./3.) C C*** Initialize various sums over target orientation. DO 27 JO=1,IORTH QSCSUM(JO)=0. QABSUM(JO)=0. QEXSUM(JO)=0. QBKSUM(JO)=0. QPHSUM(JO)=0. GSUM(JO)=0. 27 CONTINUE DO 30 IDIR=1,NSCAT S1111(IDIR)=0. S1212(IDIR)=0. S2121(IDIR)=0. S2222(IDIR)=0. CX1112(IDIR)=(0.,0.) CX1121(IDIR)=(0.,0.) CX1122(IDIR)=(0.,0.) CX1221(IDIR)=(0.,0.) CX1222(IDIR)=(0.,0.) CX2122(IDIR)=(0.,0.) 30 CONTINUE C C Loop over target orientations (target rotations in Lab Frame) C (actually, loop over lab rotations in Target Frame) C IDIR=0 IDIRWR=0 DO 110 ITHETA=1,NTHETA THETAD=THETA(ITHETA)*DEGRAD DO 100 IBETA=1,NBETA BETAD=BETA(IBETA)*DEGRAD DO 97 IPHI=1,NPHI PHID=PHI(IPHI)*DEGRAD IDIR=IDIR+1 C C**** Obtain propagation vectors, polarization vectors, and scattering C vectors in Target Frame for desired target rotation in Lab Frame C CALL ROTATE(A1,A2,AK1,CXE01,CXE02,BETA(IBETA), & THETA(ITHETA),PHI(IPHI),EN0R, & CXE01R,CXE02R, & MXSCA,NSCAT,ENSC,EM1,EM2, & AKSR,EM1R,EM2R) C DO 40 J=1,3 AKR(J)=AK1*EN0R(J) 40 CONTINUE C CALL NAMER(IWAV,IRAD,IDIR,CFLSCA,CFLAVG) C CALL GETFML(AKR,AK3,AKSR, & CALPHA,CMETHD,CXADIA,CXALPH,CXE,CXE01R,CXE02R,CXEPS, & CXF11,CXF12,CXF21,CXF22,CXXI,CXSC,CXZC,CXZW, & EM1R,EM2R, & ICOMP,ICTHM,IDVOUT,INIT,IOCC,IORTH,IPHI,IPHIM,ITASK,IXYZ0, & LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0, & MXCOMP,MXCXSC,MXMEM,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXSCA, & NAT,NAT0,NAT3,NSCAT,NX,NY,NZ, & PHI,PHIN,PIA2, & QABS,QBKSCA,QEXT,QPHA,QSCAT,QSCATG, & SCRRS1,THETAN,TOL) C C*** Compute g= DO 45 J=1,IORTH G(J)=QSCATG(J)/QSCAT(J) 45 CONTINUE C C********* Sum scattering properties over orientations **************** C C Correspondence to notation of Bohren & Huffman (1983): C Our f_jk correspond to notation of Draine (1988). Relationship between C f_jk and elements S_j of the amplitude scattering matrix as defined C by Bohren & Huffman (1983) is as follows: C S_1 = -i*f_22 C S_2 = -i*f_11 C S_3 = i*f_12 C S_4 = i*f_21 C The elements of the 4x4 Mueller matrix (see Bohren & Huffman 1983) are C related to our f_jk as follows: C S_11 = [|f_11|^2 + |f_22|^2 + |f_12|^2 + |f_21|^2]/2 C S_12 = [|f_11|^2 + |f_21|^2 - |f_22|^2 - |f_12|^2]/2 C S_13 = -Re[f_11*conjg(f_12) + f_22*conjg(f_21)] C S_14 = Im[f_22*conjg(f_21) - f_11*conjg(f_12)] C S_21 = [|f_11|^2 + |f_12|^2 - |f_22|^2 - |f_21|^2]/2 C S_22 = [|f_11|^2 + |f_22|^2 - |f_21|^2 - |f_12|^2]/2 C S_23 = Re[f_22*conjg(f_21) - f_11*conjg(f_12)] C S_24 = -Im[f_11*conjg(f_12) + f_22*conjg(f_21)] C S_31 = -Re[f_11*conjg(f_21) + f_22*conjg(f_22)] C S_32 = Re[f_22*conjg(f_12) - f_22*conjg(f_21)] C S_33 = Re[f_22*conjg(f_11) + f_12*conjg(f_21)] C S_34 = Im[f_11*conjg(f_22) + f_21*conjg(f_22)] C S_41 = -Im[f_21*conjg(f_11) + f_22*conjg(f_12)] C S_42 = Im[f_22*conjg(f_12) - f_22*conjg(f_11)] C S_43 = Im[f_22*conjg(f_11) - f_12*conjg(f_21)] C S_44 = Re[f_22*conjg(f_11) - f_12*conjg(f_21)] C C Notation internal to this program: C ij.ne.kl: CX_ijkl = \f_ij* \times f_kl (summed over orientations) C S_ijij = \f_ij* \times f_ij (summed over orientations) C 4 pure real elements: S1111,S1212,S2121,S2222 C 8 complex elements: CX1112,CX1121,CX1122,CX1221,CX1222,CX2122 C 4 redundant elts.: CX1211 = Conjg(CX1112) C CX2111 = Conjg(CX1121) C CX2112 = Conjg(CX1221) C CX2211 = Conjg(CX1122) C WG=WGTA(ITHETA,IPHI)*WGTB(IBETA) SUMWGT=0. DO 50 ND=1,NSCAT S1111(ND)=S1111(ND)+ & WG*CONJG(CXF11(ND))*CXF11(ND) S1212(ND)=S1212(ND)+ & WG*CONJG(CXF12(ND))*CXF12(ND) S2121(ND)=S2121(ND)+ & WG*CONJG(CXF21(ND))*CXF21(ND) S2222(ND)=S2222(ND)+ & WG*CONJG(CXF22(ND))*CXF22(ND) CX1112(ND)=CX1112(ND)+ & WG*CONJG(CXF11(ND))*CXF12(ND) CX1121(ND)=CX1121(ND)+ & WG*CONJG(CXF11(ND))*CXF21(ND) CX1122(ND)=CX1122(ND)+ & WG*CONJG(CXF11(ND))*CXF22(ND) CX1221(ND)=CX1221(ND)+ & WG*CONJG(CXF12(ND))*CXF21(ND) CX1222(ND)=CX1222(ND)+ & WG*CONJG(CXF12(ND))*CXF22(ND) CX2122(ND)=CX2122(ND)+ & WG*CONJG(CXF21(ND))*CXF22(ND) SUMWGT=SUMWGT+WG 50 CONTINUE C C**** Now augment sums of QABS,QBKSCA,QEXT,QSCA,G*QSCA over incident C directions and polarizations. DO 60 JO=1,IORTH QEXSUM(JO)=QEXSUM(JO)+QEXT(JO)*WG QABSUM(JO)=QABSUM(JO)+QABS(JO)*WG QBKSUM(JO)=QBKSUM(JO)+QBKSCA(JO)*WG QPHSUM(JO)=QPHSUM(JO)+QPHA(JO)*WG QSCSUM(JO)=QSCSUM(JO)+QSCAT(JO)*WG GSUM(JO)=GSUM(JO)+QSCATG(JO)*WG 60 CONTINUE C**** Choose whether or not to write out scattering properties for C current orientation; conditional may be replaced if desired. C C**** If IWRKSC=1, write scattering properties of current orientation: IF(IWRKSC.EQ.1)THEN IDIRWR=IDIRWR+1 CALL NAMER(IWAV,IRAD,IDIRWR,CFLSCA,CFLAVG) OPEN(UNIT=8,FILE=CFLSCA,STATUS='UNKNOWN') WRITE(8,FMT=9030)CSTAMP,CDESCR, & CMETHD,CALPHA,CSHAPE,NAT0 WRITE(8,FMT=9032)AEFF,WAVE,XX WRITE(8,FMT=9031)(CXRFR(ND),CXEPS(ND),ND, & ND=1,NCOMP) WRITE(8,FMT=9033)TOL,ICTHM,IPHIM,A1,A2 WRITE(8,FMT=9035)AKR,CXE01R,CXE02R WRITE(8,FMT=9040)BETAD,THETAD,PHID WRITE(8,FMT=9050)QEXT(1),QABS(1),QSCAT(1),G(1), & QBKSCA(1),QPHA(1) IF(IORTH.EQ.2)THEN QAV(1)=.5*(QEXT(1)+QEXT(2)) QAV(2)=.5*(QABS(1)+QABS(2)) QAV(3)=.5*(QSCAT(1)+QSCAT(2)) QAV(4)=(QSCATG(1)+QSCATG(2))/ & (QSCAT(1)+QSCAT(2)) QAV(5)=.5*(QBKSCA(1)+QBKSCA(2)) QAV(6)=.5*(QPHA(1)+QPHA(2)) QAV(7)=QEXT(1)-QEXT(2) QAV(8)=QPHA(1)-QPHA(2) WRITE(8,FMT=9055)QEXT(2),QABS(2),QSCAT(2), & G(2),QBKSCA(2),QPHA(2),QAV ENDIF C convert orientation angles to degrees BETAD=BETA(IBETA)*DEGRAD THETAD=THETA(ITHETA)*DEGRAD PHID=PHI(IPHI)*DEGRAD IF(IORTH.EQ.1)THEN WRITE(8,FMT=9060) DO 67 ND=1,NSCAT C convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) RVAR1=CONJG(CXF11(ND))*CXF11(ND) RVAR2=CONJG(CXF21(ND))*CXF21(ND) CXVAR1=CONJG(CXF11(ND))*CXF21(ND) WRITE(8,FMT=9070)ND,THETND,PHIND,RVAR1, & RVAR2,CXVAR1 67 CONTINUE ELSEIF(IORTH.EQ.2)THEN WRITE(8,FMT=9080) DO 73 ND=1,NSCAT C convert scAttering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) RVAR1=CONJG(CXF11(ND))*CXF11(ND) RVAR2=CONJG(CXF21(ND))*CXF21(ND) RVAR3=CONJG(CXF12(ND))*CXF12(ND) RVAR4=CONJG(CXF22(ND))*CXF22(ND) C Compute DEPOLR = depolarization ratio (defined by Bohren & Huffman, C 1983, p. 403). DENOM=RVAR1+RVAR2+RVAR3+RVAR4 DEPOLR=2.*(RVAR2+RVAR3)/DENOM C Compute PPOL = degree of linear polarization of scattered light C for incident unpolarized light (defined by Bohren & C Huffman, 1983, p. 382). Note: prior to 93.12.15 this C was NOT calculated correctly. Incorrect formula C was PPOL=(RVAR2-RVAR1+RVAR4-RVAR3)/DENOM C Corrected by BTD 93.12.15: PPOL=(RVAR3-RVAR1+RVAR4-RVAR2)/DENOM WRITE(8,FMT=9090)ND,THETND,PHIND,RVAR1, & RVAR2,RVAR3,RVAR4, & DEPOLR,PPOL 73 CONTINUE ENDIF CLOSE(8) ENDIF 97 CONTINUE C C*** Have completed loop over target orientation IPHI C 100 CONTINUE C C*** Have completed loop over target orientation IBETA C 110 CONTINUE C C*** Have completed loop over target orientation ITHETA C*** Have completed all loops over target orientations C C*** Write dielectric information to 'mtable' FREQ=1.E4/WAVE CXEN=SQRT(CXEPS(1)) WRITE(11,FMT=9058)WAVE,FREQ,CXEN,CXEPS(1) C C******** Now calculate and print out orientational averages ********** C OPEN(UNIT=8,FILE=CFLAVG,STATUS='UNKNOWN') WRITE(8,FMT=9030)CSTAMP,CDESCR,CMETHD,CALPHA,CSHAPE,NAT0 WRITE(8,FMT=9032)AEFF,WAVE,XX WRITE(8,FMT=9031)(CXRFR(ND),CXEPS(ND),ND,ND=1,NCOMP) WRITE(8,FMT=9033)TOL,ICTHM,IPHIM,A1,A2 RVAR1=AK1 RVAR2=0. RVAR3=0. WRITE(8,FMT=9037)RVAR1,RVAR2,RVAR3,CXE01,CXE02 WRITE(8,FMT=9042)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH C RVAR1=GSUM(1)/QSCSUM(1) WRITE(8,FMT=9050)QEXSUM(1),QABSUM(1),QSCSUM(1), & RVAR1,QBKSUM(1),QPHSUM(1) IF(IORTH.EQ.1)THEN WRITE(8,FMT=9060) DO 120 ND=1,NSCAT C Convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) WRITE(8,FMT=9070)ND,THETND,PHIND,S1111(ND),S2121(ND), & CX1121(ND) 120 CONTINUE C C*** Write orientationally-averaged Q values (except Qpha) to 'qtable': OPEN(UNIT=10,FILE='qtable',STATUS='OLD') DO 121 ND=1,ILIN10 READ(10,*) 121 CONTINUE C C 93/01/15 (BTD): Replace following line: C WRITE(10,9056)AEFF,WAVE,QEXT(1),QABS(1),QSCAT(1),G(1), C & QBKSCA(1) C RVAR1=GSUM(1)/QSCSUM(1) WRITE(10,9056)AEFF,WAVE,QEXSUM(1),QABSUM(1),QSCSUM(1), & RVAR1,QBKSUM(1) ILIN10=ILIN10+1 CLOSE(10) OPEN(UNIT=12,FILE='qtable2',STATUS='OLD') DO 122 ND=1,ILIN12 READ(12,*) 122 CONTINUE WRITE(12,9057)AEFF,WAVE,QPHSUM(1) ILIN12=ILIN12+1 CLOSE(12) ELSEIF(IORTH.EQ.2)THEN RVAR1=GSUM(2)/QSCSUM(2) QAV(1)=.5*(QEXSUM(1)+QEXSUM(2)) QAV(2)=.5*(QABSUM(1)+QABSUM(2)) QAV(3)=.5*(QSCSUM(1)+QSCSUM(2)) QAV(4)=(GSUM(1)+GSUM(2))/(QSCSUM(1)+QSCSUM(2)) QAV(5)=.5*(QBKSUM(1)+QBKSUM(2)) QAV(6)=.5*(QPHSUM(1)+QPHSUM(2)) QAV(7)=QEXSUM(1)-QEXSUM(2) QAV(8)=QPHSUM(1)-QPHSUM(2) WRITE(8,FMT=9055)QEXSUM(2),QABSUM(2),QSCSUM(2),RVAR1, & QBKSUM(2),QPHSUM(2),QAV C C*** Write orientationally-averaged Q values (except Qpha) to 'qtable': OPEN(UNIT=10,FILE='qtable',STATUS='OLD') DO 123 ND=1,ILIN10 READ(10,*) 123 CONTINUE WRITE(10,9056)AEFF,WAVE,QAV(1),QAV(2),QAV(3),QAV(4), & QAV(5) ILIN10=ILIN10+1 CLOSE(10) OPEN(UNIT=12,FILE='qtable2',STATUS='OLD') DO 124 ND=1,ILIN12 READ(12,*) 124 CONTINUE WRITE(12,9057)AEFF,WAVE,QAV(6),QAV(7),QAV(8) ILIN12=ILIN12+1 CLOSE(12) C WRITE(8,FMT=9080) DO 130 ND=1,NSCAT C Convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) C Compute depolarization ratio (Bohren & Huffman p. 403) DENOM=S1111(ND)+S1212(ND)+S2121(ND)+S2222(ND) DEPOLR=2.*(S1212(ND)+S2121(ND))/DENOM C Compute PPOL = degree of linear polarization of scattered light C for incident unpolarized light (Bohren & Huffman p. 382) C Note: prioir to 93.12.15 this was NOT calculated correctly. Incorrect C formula was PPOL=(S2121(ND)+S2222(ND)-S1111(ND)-S1212(ND))/DENOM C Corrected by BTD 93.12.15: PPOL=(S1212(ND)-S1111(ND)+S2222(ND)-S2121(ND))/DENOM WRITE(8,FMT=9090)ND,THETND,PHIND,S1111(ND),S2121(ND), & S1212(ND),S2222(ND),DEPOLR,PPOL 130 CONTINUE ENDIF CLOSE(8) 140 CONTINUE C C*** Have completed loop over target size C 150 CONTINUE C C*** Have completed loop over wavelengths C STOP 'Program DDSCAT terminates normally.' C 9000 FORMAT(' >DDSCAT --- ',A,' compiled with MXNX=',I3,' MXNY=',I3, & ' MXNZ=',I3,' MXMEM=',I1) 9011 FORMAT(7X,I7,' = NAT = number of dipoles in extended target',/, & 7X,3I4,' = x,y,z length of extended target (Targ. Frame)') 9020 FORMAT(' DDSCAT --- ',A,/,' TARGET ---',A,/, & ' ',A,' --- method of solution ',/, & ' ',A,' --- prescription for polarizabilities',/, & ' ',A,' --- shape ',/, & ' NAT0 = ',I7,' = number of dipoles',/, & ' ',A,' --- source of dielectric function') 9021 FORMAT(' ',A) 9030 FORMAT(' DDSCAT --- ',A,/,' TARGET ---',A,/, & ' ',A,' --- method of solution ',/, & ' ',A,' --- prescription for polarizabilies',/, & ' ',A,' --- shape ',/, & ' NAT0 = ',I7,' = number of dipoles') 9031 FORMAT(' n= (',F7.4,' , ',F7.4,'), epsilon= (',F8.4,' , ',F7.4, & ') for material',I2) 9032 FORMAT(' AEFF=',F10.5,' =',' effective radius (physical units)', & /,' WAVE=',F10.5,' = wavelength (physical units)',/, & ' K*A=',F10.5,' = 2*pi*aeff/lambda') 9033 FORMAT(' TOL=',1P,E10.3,' = error tolerance for CCG method',/, & ' ICTHM=',I3,' = theta values used in comp. of Qsca,g',/, & ' IPHIM=',I3,' = phi values used in comp. of Qsca,g',/, & 1X,0P,'(',3F7.4,') = target axis A1 in Target Frame',/, & 1X,'(',3F7.4,') = target axis A2 in Target Frame') 9035 FORMAT(1X,'(',3F7.4,') = k vector (latt. units) in TF',/, & 1X,3('(',F7.4,',',F7.4,')'),' inc. pol. vec. 1 in TF',/, & 1X,3('(',F7.4,',',F7.4,')'),' inc. pol. vec. 2 in TF') 9037 FORMAT(1X,'(',3F7.4,') = k vector (latt. units) in Lab Frame',/, & 1X,3('(',F7.4,',',F7.4,')'),' inc. pol. vec. 1 in LF',/, & 1X,3('(',F7.4,',',F7.4,')'),' inc. pol. vec. 2 in LF') 9040 FORMAT(' BETA =',F7.3,' = rotation of target around A1',/, & ' THETA=',F7.3,' = angle between A1 and k',/, & ' PHI =',F7.3,' = rotation of A1 around k') 9042 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations') 9043 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff wave Q_ext Q_abs Q_sca ', & 'g= Q_bk') 9044 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff wave Q_ext Q_abs Q_sca ', & 'g= Q_bk') 9045 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff wave Q_pha') 9046 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/, & 2F8.3,' = theta_min, theta_max; NTHETA=',I2,/, & 2F8.3,' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations',/, & ' aeff wave Q_pha Q_pol Q_cpol') 9048 FORMAT(' wave(um) f(cm-1) Re(m) Im(m) Re(eps) ', & 'Im(eps)') 9050 FORMAT(' >DDSCAT',3X,'Qext',6X,'Qabs',6X,'Qsca',5X,'g=',4X, & 'Qbk',7X,'Qpha',/,3X,'JO=1:',1P,6E10.3) 9055 FORMAT(3X,'JO=2:',1P,6E10.3,/,3X,'mean:',6E10.3,/, & 3X,'Qpol=',E10.3,34X,'Qcpol=',E10.3) 9056 FORMAT(1PE10.4,E11.4,3E10.3,E11.3,E10.3) 9057 FORMAT(1PE10.4,E11.4,3E11.3) 9058 FORMAT(1PE10.4,E11.4,4E10.3) 9060 FORMAT('**** Selected scattering directions ', & ' [note: incident pol state 1 is FIXED!]',/, & ' ND THETA PHI <|f11|^2> <|f21|^2> Re', & ' Im') 9070 FORMAT(I3,2F6.1,1P,2E10.3,2E11.3) 9080 FORMAT(' ***************** Selected scattering directions *****', & '************',/, & ' ND THETA PHI <|f11|^2> <|f21|^2> <|f12|^2> ', & '<|f22|^2> DEPOLR PPOL') 9090 FORMAT(I3,2F6.1,1P,4E10.3,0P,F8.5,F9.5) C END