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 Subsequent versions of DDSCAT were 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 5a) 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 cross sections for: C a. absorption; C b. extinction; C c. scattering; C d. vector radiation pressure; C e. vector torque. C 8. User may easily specify scattering directions for which C scattering matrixs 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 95.06.15 (BTD): Changed CSTAMP from character*11 to character*22 C to include a date. C 95.06.14 (BTD): Create version 5x.1 for internal use. C Version 5x.1 computes vector force and torque C on grain due to incident radiation. C Changed QSCATG(2) to QSCAG(3,2) to allow for C asymmetric scattering by target. C Added variable QTRQSC(1-3,1-2) and added to C argument list of GETFML C Replaced variable GSUM(1-2) to QSCGSUM(1-3,1-2) C Added variables QTRQAB(1-3,1-2) and C QTRQABSUM(1-3,1-2) C 95.06.20 (PJF+: Introduce changes into DDSCAT in order to support C BTD): use of modular iterative solvers (e.g., CCG). C : add CMDSOL to argument list to specify method C of solution C : add CMDTRQ to argument list of REAPAR and GETFML C 95.07.10 (BTD): deleted redundant declarations of NAT,NAT0,NAT3, C NX,NY,NZ (already declared in COMMON/M6/ C 95.07.20 (BTD): added new scratch array SCRRS2(MXNAT) for use in C torque calculations in SCAT C 95.07.27 (BTD): Added new scratch variable CXSCR1 for GETFML to C use for storing first polarization solution while C solving for second. C 96.01.26 (BTD): Replaced IX0,IY0,IZ0 by IXYZ0 C (with simultaneous change in TARGET) C 96.02.23 (BTD): Modified format statements to get additional digits C for a_1, a_2, k vector, and pol. vectors. C 96.11.05 (BTD): Modified to remove various computations to C subroutines GETMUELLER and WRITESCA C 96.11.14 (PJF): Add TIMERS, MXTIMERS, NTIMERS C Add cbinflag, cbinfile, IOBIN (binary file) C Add cbinflag to formal parameters of REAPAR routine C Remove GETSET routine calls C 96.11.15 (PJF): Add cnetflag, cnetfile to reapar call. Add netCDF close C 96.11.21 (BTD): Declared IDNC,NTIMERS,etc. C Removed call to NCCLOS C Added call to WRITESCA in order to call WRITENET in C order to call NCCLOS in order to close netCDF file C 96.12.02 (BTD): initialized smori to zero before looping over C orientations C eliminated variables S1212,S222,CX1112,CX1122, C CX1221,CX1222,CX2122 C 97.07.24 (BTD): Changed argument list for GETMUELLER, replacing C PHI by PHIN, in order to correct error in computation C of Mueller matrix elements and amplitude scattering C matrix elements (see GETMUELLER). 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 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 PARAMETER(MXMEM=0) C PARAMETER(MXMEM=1) C C*** Set parameters MXNX,MXNY,MXNZ ************************************* INTEGER MXNX,MXNY,MXNZ PARAMETER(MXNX=8,MXNY=8,MXNZ=8) C PARAMETER(MXNX=12,MXNY=12,MXNZ=12) 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=24,MXNY=24,MXNZ=24) C PARAMETER(MXNX=25,MXNY=25,MXNZ=25) C PARAMETER(MXNX=32,MXNY=32,MXNZ=32) C PARAMETER(MXNX=36,MXNY=36,MXNZ=36) 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=72,MXNY=72,MXNZ=72) 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=100,MXBETA=100,MXPHI=100) INTEGER MXRAD,MXWAV,MXWAVT PARAMETER(MXRAD=100,MXWAV=100,MXWAVT=1500) INTEGER MXSCA PARAMETER(MXSCA=1000) INTEGER MXTIMERS,NTIMERS PARAMETER (MXTIMERS=20, NTIMERS=12) 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 C increase from 7*MXN3 required by PETR C to 10*MXN3 required by PIMCBICGSTAB method C (PIM Complex BiConjugate Gradient with Stabilization) PARAMETER(MXCXSC=10*MXN3) C C C Common variables (required for communication with the C conjugate gradient package): C CHARACTER CMDFFT*6 INTEGER IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT3,NAT0,NX,NY,NZ INTEGER*2 IOCC(MXNAT) REAL AKR(3) COMPLEX & CXADIA(MXN3), & CXZC(MXNX+1,MXNY+1,MXNZ+1,6), & CXZW(2*MXNX,2*MXNY,2*MXNZ,3+MXMEM) C----------------------------------------------------------------------- COMMON/M1/AKR COMMON/M2/CXADIA COMMON/M3/CXZC COMMON/M4/CXZW COMMON/M5/IOCC COMMON/M6/MXNATF,MXNXF,MXNYF,MXNZF,NAT,NAT3, & NAT0,NX,NY,NZ,MXMEMF,MXN3F,IDVOUT,CMDFFT C----------------------------------------------------------------------- C C Local variables: C CHARACTER CMDSOL*6,CMDTRQ*6 C COMPLEX CXEN C REAL & AEFF,AK1,AK3,BETAD,BETAMI,BETAMX,BETMID,BETMXD, & DEGRAD,FREQ, & PHID,PHIMAX,PHIMID,PHIMIN,PHIMXD,PI,PIA2, & THETAD,THETMI,THETMX,THTMID,THTMXD,TOL, & WAVE,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 QSCGSUM(1-3,1-2)=running (weighted) sum of g(1-3)*Q_sca C over target orientations, for incident polarizations 1-2, C where C g(1)= for scattered radiation C g(2)= for scattered radiation C g(3)= for scattered radiation C theta is measured relative to incident direction C phi is measured relative to x,y plane (Lab Frame) 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 QTRQSCSUM(1-3,1-2)=running sum of QTRQSC(1-3,1-2) 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,IDNC,IDVERR,IDVSHP,INIT,IOPAR,IORI,IORTH, & ILIN10,ILIN12,IOBIN,IOSHP,IPHI,IPHIM,IRAD,ITASK,ITHETA, & IWAV,IWRKSC, & J,JO, & LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0, & NBETA,NCOMP,NORI,NPHI,NRAD,NSCAT,NTHETA, & NWAV 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,CSHAPE*6,CSTAMP*22, & CFLAVG*13,CFLPAR*13,CFLSCA*14,CFLSHP*13,CDESCR*67, & cbinflag*6, cbinfile*14, cnetflag*6, cnetfile*14 CHARACTER*40 CFLEPS(MXCOMP) C .. C .. Local Arrays .. COMPLEX & CX1121(MXSCA), & CXALPH(MXN3), & CXE(MXN3),CXE01(3),CXE02(3), & CXE01R(3),CXE02R(3),CXEPS(MXCOMP), & CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),CXF22(MXSCA), & CXRFR(MXCOMP),CXS1(MXSCA),CXS2(MXSCA),CXS3(MXSCA),CXS4(MXSCA), & CXSC(MXCXSC),CXSCR1(MXN3),CXXI(MXN3) 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 CX1121(1-NSCAT)=sum of f11* \times f21 over orientations C C Scratch arrays: C CXSC(1-MXCXSC)=scratch array used by DDACCG C INTEGER*2 ICOMP(MXN3),ISCR1(MXN3),IXYZ0(MXNAT,3) C C IXYZ0(1-NAT0,1-3)=x/d,y/d,z/d for dipoles 1-NAT0 in real target C C Scratch arrays: C ISCR1,ISCR2,ISCR3,ISCR4 (used by EXTEND) C C REAL & A1(3),A2(3),AEFFA(MXRAD),AKSR(3,MXSCA), & BETA(MXBETA), & E1A(MXWAVT),E2A(MXWAVT),EM1(3,MXSCA),EM1R(3,MXSCA), & EM2(3,MXSCA),EM2R(3,MXSCA),EN0R(3),ENSC(3,MXSCA), & PHI(MXPHI),PHIN(MXSCA), & QABS(2),QABSUM(2),QBKSCA(2),QBKSUM(2),QEXT(2), & QEXSUM(2),QPHA(2),QPHSUM(2),QSCAT(2),QSCAG(3,2),QSCGSUM(3,2), & QSCSUM(2),QTRQAB(3,2),QTRQSC(3,2), & QTRQABSUM(3,2),QTRQSCSUM(3,2), & S1111(MXSCA),S2121(MXSCA), & SCRRS1(MXN3),SCRRS2(MXNAT),SHPAR(6),SM(MXSCA,4,4), & SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA), & WAVEA(MXWAV),WGTA(MXTHET,MXPHI),WGTB(MXBETA),WVA(MXWAVT), & TIMERS(MXTIMERS) 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 S2121(1-NSCAT)= |f21|**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 SCRRS2(1-NAT)=scratch array used by SCAT C C External subroutines: EXTERNAL DIELEC,EVALA,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 Initialize various variables in common blocks which are derived C from parameters: C MXNATF=MXNAT MXNXF=MXNX MXNYF=MXNY MXNZF=MXNZ MXMEMF=MXMEM MXN3F=MXN3 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 C C**** Read program control parameters C CMDFFT=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(cbinflag, cnetflag, AEFFA, & BETAMI,BETAMX, & CALPHA,CDIEL,CFLEPS,CFLPAR,CMDSOL,CMDFFT,CSHAPE, & CMDTRQ,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,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(CMDFFT,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 why not use cbinfile ? and open it in writesca ? IOBIN=51 OPEN(UNIT=IOBIN,file='dd.bin',form='unformatted') C C*** Write out header for 'qtable', 'qtable2', and 'mtable' files WRITE(10,FMT=9020)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0,CDIEL ILIN10=7 WRITE(10,FMT=9021)(CFLEPS(J),J=1,NCOMP) ILIN10=ILIN10+NCOMP WRITE(12,FMT=9020)CSTAMP,CDESCR,CMDFFT,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,CMDFFT,CALPHA,CSHAPE,NAT0,CDIEL WRITE(11,FMT=9021)(CFLEPS(J),J=1,NCOMP) WRITE(11,FMT=9048) C C*** Loop over wavelengths: C DO 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 J=1,NCOMP CXRFR(J)=SQRT(CXEPS(J)) WRITE(IDVOUT,FMT=9031)CXRFR(J),CXEPS(J),J ENDDO C C**** Loop over grain sizes: C DO 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 JO=1,IORTH QSCSUM(JO)=0. QABSUM(JO)=0. QEXSUM(JO)=0. QBKSUM(JO)=0. QPHSUM(JO)=0. DO J=1,3 QSCGSUM(J,JO)=0. QTRQABSUM(J,JO)=0 QTRQSCSUM(J,JO)=0. ENDDO ENDDO DO IDIR=1,NSCAT S1111(IDIR)=0. S2121(IDIR)=0. CX1121(IDIR)=(0.,0.) ENDDO DO J=1,4 DO JO=1,4 DO IDIR=1,NSCAT SMORI(IDIR,JO,J)=0. ENDDO ENDDO ENDDO C C Loop over target orientations (target rotations in Lab Frame) C (actually, loop over lab rotations in Target Frame) C IORI=0 DO ITHETA=1,NTHETA THETAD=THETA(ITHETA)*DEGRAD DO IBETA=1,NBETA BETAD=BETA(IBETA)*DEGRAD DO IPHI=1,NPHI PHID=PHI(IPHI)*DEGRAD IORI=IORI+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 J=1,3 AKR(J)=AK1*EN0R(J) ENDDO C CALL NAMER(IWAV,IRAD,IORI,CFLSCA,CFLAVG) C CALL GETFML(AKR,AK3,AKSR, & CALPHA,CMDSOL,CMDFFT,CMDTRQ, & CXADIA,CXALPH,CXE,CXE01R,CXE02R,CXEPS, & CXF11,CXF12,CXF21,CXF22,CXXI,CXSC,CXSCR1,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,QSCAG,QTRQAB,QTRQSC, & SCRRS1,SCRRS2,THETAN,TOL,TIMERS, MXTIMERS, NTIMERS) C C C Compute Mueller scattering matrix C CALL GETMUELLER(IBETA,IORTH,IPHI,ITHETA,MXBETA, & MXSCA,MXPHI,MXTHET,NSCAT, & CMDTRQ,PHIN,SM,SMORI, & S1111,S2121, & CX1121,CXE01,CXE02, & CXF11,CXF12,CXF21,CXF22, & CXS1,CXS2,CXS3,CXS4, & QABS,QABSUM,QBKSCA,QBKSUM,QEXSUM, & QEXT,QPHA,QPHSUM,QSCAG,QSCAT, & QSCGSUM,QSCSUM,QTRQAB,QTRQABSUM, & QTRQSC,QTRQSCSUM,WGTA,WGTB) C C---------------------------------------------------------------------- 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 CALL WRITESCA(MXNX,MXNY,MXNZ, & NX,NY,NZ,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD,ITHETA, & IBETA,IPHI,THETA,MXTHET,NTHETA,BETA,MXBETA, & NBETA,PHI,MXPHI,NPHI, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMORI,THETAN, & CX1121,CXE01,CXE01R,CXE02, & CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) ENDIF ENDDO C C*** Have completed loop over target orientation IPHI C ENDDO C C*** Have completed loop over target orientation IBETA C ENDDO 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 Set IORI=0 prior to calling WRITESCA in order to print orientational C averages C IORI=0 C C******** Now calculate and print out orientational averages ********** C CALL WRITESCA(MXNX,MXNY,MXNZ, & NX,NY,NZ,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD,ITHETA, & IBETA,IPHI,THETA,MXTHET,NTHETA,BETA,MXBETA, & NBETA,PHI,MXPHI,NPHI, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMORI,THETAN, & CX1121,CXE01,CXE01R,CXE02, & CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) ENDDO C C*** Have completed loop over target size C ENDDO C C*** Have completed loop over wavelengths C C close binary file (this perhaps should be moved to writesca) CLOSE(IOBIN) IF(CNETFLAG.NE.'NOTCDF')THEN CALL WRIMSG('DDSCAT',' close NetCDF') CNETFLAG='NCCLOS' CALL WRITESCA(MXNX,MXNY,MXNZ, & NX,NY,NZ,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD,ITHETA, & IBETA,IPHI,THETA,MXTHET,NTHETA,BETA,MXBETA, & NBETA,PHI,MXPHI,NPHI, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMORI,THETAN, & CX1121,CXE01,CXE01R,CXE02, & CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) ENDIF 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) 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') 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)') 9058 FORMAT(1PE10.4,E11.4,4E10.3) C END