text:
# This document can be cited as: Draine, B.T., and Flatau, P.J. 2000, “User Guide for the Discrete Dipole Approximation Code DDSCAT (Version 5a10)”, http://arxiv.org/abs/astro-ph/0008151v4 User Guide for the Discrete Dipole Approximation Code DDSCAT (Version 5a10)

## 1 Introduction

DDSCAT.5a is a Fortran software package to calculate scattering and absorption of electromagnetic waves by targets with arbitrary geometries using the “discrete dipole approximation” (DDA). In this approximation the target is replaced by an array of point dipoles (or, more precisely, polarizable points); the electromagnetic scattering problem for an incident periodic wave interacting with this array of point dipoles is then solved essentially exactly. The DDA (sometimes referred to as the “coupled dipole approximation”) was apparently first proposed by Purcell & Pennypacker (1973). DDA theory was reviewed and developed further by Draine (1988), Draine & Goodman (1993), and recently reviewed by Draine & Flatau (1994) and Draine (2000).

DDSCAT.5a is a Fortran implementation of the DDA developed by the authors. The previous version, DDSCAT.5a9, was released 1998 December 15. The current version, DDSCAT.5a10, released 2000 June 15, adds a new “multisphere” target option. DDSCAT.5a is intended to be a versatile tool, suitable for a wide variety of applications ranging from interstellar dust to atmospheric aerosols. As provided, DDSCAT.5a10 should be usable for many applications without modification, but the program is written in a modular form, so that modifications, if required, should be fairly straightforward.

The authors make this code openly available to others, in the hope that it will prove a useful tool. We ask only that:

* If you publish results obtained using DDSCAT, please consider acknowledging the source of the code.
* If you discover any errors in the code or documentation, please promptly communicate them to the authors.
* You comply with the “copyleft” agreement (more formally, the GNU General Public License) of the Free Software Foundation: you may copy, distribute, and/or modify the software identified as coming under this agreement. If you distribute copies of this software, you must give the recipients all the rights which you have. See the file doc/copyleft distributed with the DDSCAT software.

We also strongly encourage you to send email to the authors identifying yourself as a user of DDSCAT; this will enable the authors to notify you of any bugs, corrections, or improvements in DDSCAT.

The current version, DDSCAT.5a10 , uses the DDA formulae from Draine (1988). The code incorporates Fast Fourier Transform (FFT) methods (Goodman, Draine, & Flatau 1991). The “lattice dispersion relation” (LDR) prescription (Draine & Goodman 1993) is used for determining dipole polarizabilities.

This User Guide assumes that you have already obtained the Fortran source code for DDSCAT.5a10 either via the World Wide Web (http://www.astro.princeton.edu/$``$draine) or via anonymous ftp following the instructions in the README file.111To obtain the README file: (1) anonymous ftp to astro.princeton.edu , (2) cd draine/scat/DDA/ver5a, and (3) get README.  We refer you to the list of references at the end of this document for discussions of the theory and accuracy of the DDA \[first see the recent reviews by Draine and Flatau (1994) and Draine (2000)\]. In §4 we describe the principal changes between DDSCAT.5a10 and the previous releases.222 The first “official release” was DDSCAT.4b, although DDSCAT.4c – while never announced – was made available to a number of interested users. DDSCAT.5a8 was released in 1997 April. DDSCAT.5a9 was released in 1998 December. DDSCAT.5a10 was released 2000 June 15.  The succeeding sections contain instructions for:

* compiling and linking the code;
* running a sample calculation;
* understanding the output from the sample calculation;
* modifying the parameter file to do your desired calculations;
* specifying target orientation;
* changing the DIMENSIONing of the source code to accommodate your desired calculations.

The instructions for compiling, linking, and running will be appropriate for a UNIX system; slight changes will be necessary for non-UNIX sites, but they are quite minor and should present no difficulty.

Finally, this User Guide can be obtained by any of the following means:

* http://xxx.lanl.gov/abs/astro-ph/0008151 – you will be offered the options of downloading

+ Latex source
    + Postscript
    + Other formats – click on this to obtain the UserGuide as a PDF file.
* anonymous ftp from astro.princeton.edu, subdirectory draine/scat/ddscat/ver5a10.

+ For postscript, get either the uncompressed postscript file UserGuide.ps or the gzipped file UserGuide.ps.gz (remember to set the binary option in ftp before using the get command).
    + Alternatively, you can get the LaTeX source UserGuide.tex.gz and the figures f1.eps.gz, f2.eps.gz, f3.eps.gz, f4.eps.gz, and f5.eps.gz.

## 2 Applicability of the DDA

The principal advantage of the DDA is that it is completely flexible regarding the geometry of the target, being limited only by the need to use an interdipole separation $`d`$ small compared to (1) any structural lengths in the target, and (2) the wavelength $`\lambda `$. Numerical studies (Draine & Goodman 1993; Draine & Flatau 1994; Draine 2000) indicate that the second criterion is adequately satisfied for calculations of total cross sections if

$$|m|kd<1,$$
(1)
where $`m`$ is the complex refractive index of the target material, and $`k2\pi /\lambda `$, where $`\lambda `$ is the wavelength in the surrounding medium (normally taken to be vacuum, but see §3.2). However, if accurate calculations of the scattering phase function (e.g., radar or lidar cross sections) are desired, a more conservative criterion

$$|m|kd<0.5$$
(2)
will ensure that differential scattering cross sections $`dC_{\mathrm{sca}}/d\mathrm{\Omega }`$ are accurate to within a few percent of the average differential scattering cross section $`C_{\mathrm{sca}}/4\pi `$ (see Draine 2000).

Let $`V`$ be the target volume. If the target is represented by an array of $`N`$ dipoles, located on a cubic lattice with lattice spacing $`d`$, then

$$V=Nd^3.$$
(3)
We characterize the size of the target by the “effective radius”

$$a_{\mathrm{eff}}(3V/4\pi )^{1/3},$$
(4)
the radius of an equal volume sphere. A given scattering problem is then characterized by the dimensionless “size parameter”

$$xka_{\mathrm{eff}}=\frac{2\pi a_{\mathrm{eff}}}{\lambda }.$$
(5)
The size parameter can be related to $`N`$ and $`|m|kd`$:

$$x\frac{2\pi a_{\mathrm{eff}}}{\lambda }=\frac{62.04}{|m|}\left(\frac{N}{10^6}\right)^{1/3}|m|kd.$$
(6)
Equivalently, the target size can be written

$$a_{\mathrm{eff}}=9.873\frac{\lambda }{|m|}\left(\frac{N}{10^6}\right)^{1/3}|m|kd.$$
(7)
Practical considerations of CPU speed and computer memory currently available on scientific workstations typically limit the number of dipoles employed to $`N<10^6`$ (see §16 for limitations on $`N`$ due to available RAM); for a given $`N`$, the limitations on $`|m|kd`$ translate into limitations on the ratio of target size to wavelength.

For calculations of total cross sections $`C_{\mathrm{abs}}`$ and $`C_{\mathrm{sca}}`$, we require $`|m|kd<1`$:

$$a_{\mathrm{eff}}<9.88\frac{\lambda }{|m|}\left(\frac{N}{10^6}\right)^{1/3}\mathrm{or}x<\frac{62.04}{|m|}\left(\frac{N}{10^6}\right)^{1/3}.$$
(8)
For scattering phase function calculations, we require $`|m|kd<0.5`$:

$$a_{\mathrm{eff}}<4.94\frac{\lambda }{|m|}\left(\frac{N}{10^6}\right)^{1/3}\mathrm{or}x<\frac{31.02}{|m|}\left(\frac{N}{10^6}\right)^{1/3}.$$
(9)
It is therefore clear that the DDA is not suitable for very large values of the size parameter $`x`$, or very large values of the refractive index $`m`$. The primary utility of the DDA is for scattering by dielectric targets with sizes comparable to the wavelength. As discussed by Draine & Goodman (1993), Draine & Flatau (1994), and Draine (2000), total cross sections calculated with the DDA are accurate to a few percent provided $`N>10^4`$ dipoles are used, criterion (1) is satisfied, and $`|m1|<2`$.

Examples illustrating the accuracy of the DDA are shown in Figs. 1-2 which show overall scattering and absorption efficiencies as a function of wavelength for different discrete dipole approximations to a sphere, with $`N`$ ranging from 304 to 59728. The DDA calculations assumed radiation incident along the (1,1,1) direction in the “target frame”. Figs. 3-4 show the scattering properties calculated with the DDA for $`x=ka=7`$. Additional examples can be found in Draine & Flatau (1994) and Draine (2000).

## 3 DDSCAT.5a

### 3.1 What Does It Calculate?

DDSCAT.5a solves the problem of scattering and absorption by an array of polarizable point dipoles interacting with a monochromatic plane wave incident from infinity. DDSCAT.5a has the capability of automatically generating dipole array representations for a variety of target geometries (see §19) and can also accept dipole array representations of targets supplied by the user (although the dipoles must be located on a cubic lattice). The incident plane wave can have arbitrary elliptical polarization (see §21), and the target can be arbitrarily oriented relative to the incident radiation (see §17). The following quantities are calculated by DDSCAT.5a :

* Absorption efficiency factor $`Q_{\mathrm{abs}}C_{\mathrm{abs}}/\pi a_{\mathrm{eff}}^2`$, where $`C_{\mathrm{abs}}`$ is the absorption cross section;
* Scattering efficiency factor $`Q_{\mathrm{sca}}C_{\mathrm{sca}}/\pi a_{\mathrm{eff}}^2`$, where $`C_{\mathrm{sca}}`$ is the scattering cross section;
* Extinction efficiency factor $`Q_{\mathrm{ext}}Q_{\mathrm{sca}}+Q_{\mathrm{abs}}`$;
* Phase lag efficiency factor $`Q_{\mathrm{pha}}`$, defined so that the phase-lag (in radians) of a plane wave after propagating a distance $`L`$ is just $`n_tQ_{\mathrm{pha}}\pi a_{\mathrm{eff}}^2L`$, where $`n_t`$ is the number density of targets.
* The 4$`\times `$4 Mueller scattering intensity matrix $`S_{ij}`$ describing the complete scattering properties of the target for scattering directions specified by the user.
* The radiation force efficiency vector $`𝐐_{pr}`$ (see §15).
* The radiation torque efficiency vector $`𝐐_\mathrm{\Gamma }`$ (see §15).

### 3.2 Application to Targets in Dielectric Media

Let $`\omega `$ be the angular frequency of the incident radiation. For many applications, the target is essentially in vacuo, in which case the dielectric function $`ϵ`$ which the user should supply to DDSCAT is the actual complex dielectric function $`ϵ_{\mathrm{target}}(\omega )`$, or complex refractive index $`m_{\mathrm{target}}(\omega )=\sqrt{ϵ_{\mathrm{target}}}`$ of the target material.

However, for many applications of interest (e.g., marine optics, or biological optics) the “target” body is embedded in a (nonabsorbing) dielectric medium, with (real) dielectric function $`ϵ_{\mathrm{medium}}(\omega )`$, or (real) refractive index $`m_{\mathrm{medium}}(\omega )=\sqrt{ϵ_{\mathrm{medium}}}`$. DDSCAT.5a is fully applicable to these scattering problems, except that:

* The “dielectric function” or “refractive index” supplied to DDSCAT.5a should be the effective dielectric function

$$ϵ(\omega )=\frac{ϵ_{\mathrm{target}}(\omega )}{ϵ_{\mathrm{medium}}(\omega )}$$
(10)
or effective refractive index:

$$m(\omega )=\frac{m_{\mathrm{target}}(\omega )}{m_{\mathrm{medium}}(\omega )}.$$
(11)
* The wavelength $`\lambda `$ specified in ddscat.par should be the wavelength in the medium:

$$\lambda =\frac{\lambda _{vac}}{m_{\mathrm{medium}}},$$
(12)
where $`\lambda _{vac}=2\pi c/\omega `$ is the wavelength in vacuo.

The absorption, scattering, extinction, and phase lag efficiency factors $`Q_{\mathrm{abs}}`$, $`Q_{\mathrm{sca}}`$, and $`Q_{\mathrm{ext}}`$ calculated by DDSCAT will then be equal to the physical cross sections for absorption, scattering, and extinction divided by $`\pi a_{\mathrm{eff}}^2`$ – e.g., the attenuation coefficient for radiation propagating through a medium with a density $`n_t`$ of scatterers will be just $`\alpha =n_tQ_{\mathrm{ext}}\pi a_{\mathrm{eff}}^2`$. Similarly, the phase lag (in radians) after propagating a distance $`L`$ will be $`n_tQ_{\mathrm{pha}}\pi a_{\mathrm{eff}}^2`$.

The elements $`S_{ij}`$ of the 4$`\times `$4 Mueller scattering matrix $`𝐒`$ calculated by DDSCAT will be correct for scattering in the medium:

$$𝐈_{\mathrm{sca}}=\left(\frac{\lambda }{2\pi r}\right)^2𝐒𝐈_{\mathrm{in}},$$
(13)
where $`𝐈_{\mathrm{in}}`$ and $`𝐈_{\mathrm{sca}}`$ are the Stokes vectors for the incident and scattered light (in the medium), $`r`$ is the distance from the target, and $`\lambda `$ is the wavelength in the medium (eq. 12). See §23 for a detailed discussion of the Mueller scattering matrix.

The time-averaged radiative force and torque (see §15) on a target in a dielectric medium are

$$𝐅_{\mathrm{rad}}=𝐐_{pr}\pi a_{\mathrm{eff}}^2u_{\mathrm{rad}},$$
(14)
$$𝚪_{\mathrm{rad}}=𝐐_\mathrm{\Gamma }\pi a_{\mathrm{eff}}^2u_{\mathrm{rad}}\frac{\lambda }{2\pi },$$
(15)
where the time-averaged energy density is

$$u_{\mathrm{rad}}=ϵ_{\mathrm{medium}}\frac{E_0^2}{8\pi },$$
(16)
where $`E_0\mathrm{cos}(\omega t+\varphi )`$ is the electric field of the incident plane wave in the medium.

## 4 What’s New?

### 4.1 DDSCAT.5a

DDSCAT.5a differs from previous versions in five major respects:

1. Use of the new Generalized Prime Factor Algorithm (GPFA) developed by Clive Temperton (1992) for FFT calculations. The GPFA algorithm is generally faster than the previous algorithms, yet requires no more memory than the algorithm of Brenner (1969). See §13.
2. Availability of a new algorithm for iterative solution of the system of complex linear equations. This algorithm is often faster than the algorithm of Petravic and Kuo-Petravic (1979) which was used through DDSCAT.4b (and which remains available as an option in DDSCAT.5a). See §14.
3. Automatic calculation of the transverse components of the radiative force on the target. See §15.
4. Capability to compute the electromagnetic torque on the target, due to absorption and scattering of light from the incident beam, as described by Draine and Weingartner (1996). See §15
5. Automatic calculation of the 4$`\times `$4 Mueller scattering matrix. See §23.
6. Improved output handling, including FORTRAN unformatted binary write option as well as a netCDF interface and IDL postprocessing code.
7. Finally, this much-improved and expanded User Guide!

We also call users’ attention to a minor but possibly confusing change: DDSCAT.5a uses a different convention for specifying the target axes $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ for certain target choices (RCTNGL, ELLIPS, …). We think the new convention is more straightforward. See the discussion in §19 below for details.

### 4.2 DSCAT.5a9

DDSCAT.5a9 differs from DDSCAT.5a8 in one respect:

* DDSCAT.5a8 contained an error in the calculation of the elements of the scattering matrix (see §23) for scattering planes with $`\varphi 0`$. While DDSCAT.5a8 computed the scattering matrix correctly for $`\varphi _s=0`$, for irregular targets only the $`S_{11}`$ element was computed correctly for scattering planes with $`\varphi _s0`$. This bug has been corrected in DDSCAT.5a9.

### 4.3 DDSCAT.5a10

DDSCAT.5a10 differs from DDSCAT.5a9 in two respects:

* The public-domain LAPACK code has been replaced by the latest routines obtained from NETLIB (http://www.netlib.org). Subroutine PRINAXIS did not function properly when the older LAPACK routines were compiled using the g77 compiler, most likely due to some variable not being SAVE-d. The new LAPACK code appears to execute properly.
* A new shape option – NSPHER – has been added, permitting a target to be defined as the union of the volumes of an arbitrary number of spheres. The locations and radii of the NSPH spheres are input via a file. Adding this option required the argument list of subroutine TARGET to be changed, and modifications were also made to subroutines REAPAR and REASHP.
* A new shape option – PRISM3 – has been added (on 2002.02.12), enabling construction of a triangular prism target. Concomitant modifications were made to subroutines TARGET and REAPAR.

## 5 Obtaining the Source Code

The easiest way to download the source code is from the DDSCAT home page:
http://www.astro.princeton.edu/$``$draine/DDSCAT.html
where you can obtain the file ddscat5a10.tar.gz – a “gzipped tarfile” containing the complete source code and documentation. This can also be obtained by anonymous ftp from astro.princeton.edu, sudirectory draine/scat/ddscat/ver5a10. Note that it is a compressed (binary) file.

The source code can be installed as follows. First, place the file ddscat5a10.tar.gz in the directory where you would like DDSCAT to reside. You should have at least 5 Mbytes of disk space available.

If you are on a Linux system, you should be able to type
tar xvzf ddscat5a10.tar.gz
which will “gunzip” the tarfile and then “extract” the files from the tarfile. If your version of “tar” doesn’t support the “z” option (e.g., you are running Solaris) then try
zcat ddscat5a10.tar.gz | tar xvf -
If neither of the above work on your system, try the two-stage procedure
gunzip ddscat5a10.tar.gz
tar xvf ddscat5a10.tar
(The disadvantage of the two-stage procedure is that it uses more disk space, since after the second step you will have the uncompressed tarfile ddscat5a10.tar – about 3.8 Mbytes – in addition to all the files you have extracted from the tarfile – another 4.6 Mbytes).

Any of the above approaches should create subdirectories src, doc, misc, and IDL. The source code will be in subdirectory src, and documentation in subdirectory doc.

If you are running Windows on a PC, you will need the “winzip” program, which can be downloaded from http://www.winzip.com. winzip should be able to “unzip” the gzipped tarfile ddscat5a10.tar.gz and “extract” the various files from it automatically.

## 6 Compiling and Linking

In the discussion below, it will be assumed that the source code for DDSCAT.5a has been installed in a directory DDA/src . To compile the code on a Unix system, position yourself in the directory DDA/src. The Makefile supplied has compiler options appropriate for Sun Fortran under Solaris 2.x. If you have a different Fortran compiler, you will probably need to edit Makefile to provide the desired compiler options. Makefile contains samples of compiler options for selected operating systems, including HP AUX, IBM AIX, and SGI IRIX operating systems.

DDSCAT can be compiled on Windows systems using standard Fortran compilers (e.g., Compaq or Microsoft Visual Fortran), and run successfully. Note, however, that many of the instructions below are specifically for Unix or Linux operating systems – some steps may have to be done differently on a non-Unix platform.

So far as we know, there are only two operating-system-dependent aspects of DDSCAT.5a: (1) the device number to use for “standard output”, and (2) the TIMEIT routine. There is, in addition, one installation-dependent aspect to the code: the procedure for linking to the netCDF library (see §10.3 for discussion of the possibility of writing binary files using the machine-independent netCDF standard).

### 6.1 Device Numbers IDVOUT and IDVERR

The variables IDVOUT and IDVERR specify device numbers for “running output” and “error messages”, respectively. Normally these would both be set to the device number for “standard output” (e.g., writing to the screen if running interactively). The variables IDVOUT and IDVERR are set by DATA statements in the “main” program DDSCAT.f and in the output routine WRIMSG (file wrimsg.f). Under Sun Fortran, DATA IDVOUT/0/ results in unbuffered output to “standard output”; unbuffered output (if available) is nice so that if you choose to direct your output to a file (e.g., using ddscat \>& ddscat.out &) the output file will contain up-to-date information. Other operating systems or compilers may not support this, and you may need to edit DDSCAT.f to change the two data statements to DATA IDVOUT/6/ and DATA IDVERR/6/, and edit wrimsg.f to change DATA IDVOUT/0/ to DATA IDVOUT/6/.

### 6.2 Subroutine TIMEIT

The only other operating system-dependent part of DDSCAT.5a is the single subroutine TIMEIT. Several versions of TIMEIT are provided:

* timeit\_sun.f uses the SunOS system call etime
* timeit\_convex.f uses the Convex OS system call etime
* timeit\_cray.f uses the system call second
* timeit\_hp.f uses the HP-AUX system calls sysconf and times
* timeit\_ibm6000.f uses the AIX system call mclock
* timeit\_osf.f uses the DEC OSF system call etime
* timeit\_sgi.f uses the IRIX system call etime
* timeit\_vms.f uses the VMS system calls LIB$INIT\_TIMER and LIB$SHOW\_TIMER
* timeit\_titan.f uses the system call cputim
* timeit\_null.f is a dummy routine which provides no timing information, but can be used under any operating system.

You must compile and link one of the timeit\_xxx.f routines, possibly after modifying it to work on your system; timeit\_null.f is the easiest choice, but it will return no timing information.333 Non-Unix sites: The VMS-compatible version of TIMEIT is included in the file SRC9.FOR. For non-VMS sites, you will need to replace this version of TIMEIT with one of the other versions, which are to be found in the file MISC.FOR. If you are having trouble getting this to work, choose the “dummy” version of TIMEIT from MISC.FOR: this will return no timing information, but is platform-independent.

### 6.3 Leaving the netCDF Capability Disabled

The Makefile supplied with the distribution of DDSCAT.5a is set up to link to a “dummy” subroutine dummywritenet.f instead of subroutine writenet.f, in order to minimize problems during initial compilation and linking. The “dummy” routine has no functionality, other than bailing out with a fatal error message if the user makes the mistake of trying to specify one of the netCDF options (ALLCDF or ORICDF). First-time users of DDSCAT.5a should not try to use the netCDF option – simply skip this section, and specify option NOTCDF in ddscat.par. After successfully using DDSCAT.5a you can return to §6.4 to try to enable the netCDF capability.

### 6.4 Enabling the netCDF Capability

Subroutine WRITENET (file writenet.f) provides the capability to output binary data in the netCDF standard format (see §10.3). In order to use this routine (instead of dummywritenet.f), it is necessary to take the following steps:444 Non-UNIX sites: You need to replace the “dummy” version of SUBROUTINE WRITENET in SRC1.FOR with the version provided in MISC.FOR. You will also need to consult your systems administrator to verify that the netCDF library has already been installed on your system, and to find out how to link to the library routines.

1. Have netCDF already installed on your system (check with your system administrator).
2. Find out where netcdf.inc is located and edit the include statement in writenet.f to specify the correct path to netcdf.inc.
3. Find out where the libnetcdf.a library is located, and edit the Makefile so that the variable LIBNETCDF specifies the correct path to this library.
4. Edit the Makefile to “comment out” (with a # symbol in column 1) the line writenet = dummywritenet and “uncomment” (remove the # symbol) the line writenet = writenet so that writenet.f will be compiled instead of the dummy routine dummywritenet.f.

### 6.5 Compiling and Linking…

After suitably editing the Makefile (while still positioned in DDA/src) simply type555 Non-Unix sites: see §6.6  make ddscat which should create an executable file DDA/src/ddscat . The resulting executable will not have netCDF capability, and will contain timing instructions compatible with the Solaris 1.x or 2.x operating systems on Sun computers, as well as several other versions of Unix. To add netCDF capability, see §6.4. To replace the Sun-compatible timing routine with another, see §6.2.

### 6.6 Installation on Non-Unix Systems

DDSCAT.5a is written in standard Fortran-77 plus the DO ... ENDDO extension which appears to be supported by all current Fortran-77 compatible compilers. It is possible to run DDSCAT on non-Unix systems. If the Unix “make” utility is not available, here in brief is what needs to be accomplished:

All of the necessary Fortran code to compile and link DDSCAT.5a is included in the following files: SRC0.FOR, SRC1.FOR, SRC2.FOR, SRC3.FOR, SRC4.FOR, SRC5.FOR, SRC6.FOR, SRC7.FOR, SRC8.FOR, SRC9.FOR, CGCOMMON.FOR, GPFA.FOR, LAPACK.FOR, and PIM.FOR. There is an additional file MISC.FOR, but this is not needed for “basic” use of the code (see below).

The main program DDSCAT is found in SRC0.FOR. It calls a number of subroutines, which are included in the other \*.FOR files.

Three of the subroutines exist in more than one version. The “default” version of each is located in the file SRC1.FOR:

* Select the version of SUBROUTINE WRITENET from SRC1.FOR (do not use the version in MISC.FOR). The resulting code will not support netCDF capability. (If netCDF capability is required, you will need to install netCDF libraries on your system – see §6.4).
* Select the version of SUBROUTINE C3DFFT from SRC1.FOR (do not use the version in MISC.FOR)
* There is one version of SUBROUTINE TIMEIT included in SRC1.FOR, and a number of additional versions in MISC.FOR. Use the version from SRC1.FOR, which begins

```
      SUBROUTINE TIMEIT(CMSGTM,DTIME)
C
C     timeit_null
C
C This version of timeit is a dummy which does not provide any
C timing information.

```

This version does not use any system calls, and therefore the code should compile and link without problems; however, when you run the code, it will not report any timing information reporting how much time was spent on different parts of the calculation.

If you wish to obtain timing information, you will need to find out what system calls are supported by your operating system. You can look at the other versions of SUBROUTINE TIMEIT in MISC.FOR to see how this has been done under VMS and various version of Unix.

Once you have selected the appropriate versions of WRITENET, C3DFFT, and TIMEIT, you can simply compile and link just as you would with any other Fortran code with a number of modules. You should end up with an executable with a (system-dependent) name like DDSCAT.EXE.

In addition to program DDSCAT, we provide two other Fortran 77 programs which may be useful. Program CALLTARGET can be used to call the target generation routines which may be helpful for testing purposes. Program TSTFFT is useful for comparing speeds of different FFT options.

## 7 Moving the Executable

Now reposition yourself into the directory DDA (e.g., type cd ..), and copy the executable from src/ddscat to the DDA directory by typing

cp src/ddscat ddscat

This should copy the file DDA/src/ddscat to DDA/ddscat . Similarly, copy the sample parameter file ddscat.par and the file diel.tab to the DDA directory by typing

cp doc/ddscat.par ddscat.par cp doc/diel.tab diel.tab

## 8 The Parameter File ddscat.par

The directory DDA should now contain a sample file ddscat.par which provides parameters to the program ddscat. As provided (see AppendixA), the file ddscat.par is set up to calculate scattering by a 8$`\times `$6$`\times `$4 rectangular array of 192 dipoles, with an effective radius $`a_{\mathrm{eff}}=1\mu \mathrm{m}`$, at a wavelength of $`6.2832\mu \mathrm{m}`$ (for a “size parameter” $`2\pi a_{\mathrm{eff}}/\lambda =1`$).

The dielectric function of the target material is provided in the file diel.tab, which is a sample file in which the refractive index is set to $`m=1.33+0.01i`$ at all wavelengths; the name of this file is provided to ddscat by the parameter file ddscat.par.

The sample parameter file as supplied calls for the new GPFA FFT routine (GPFAFT) of Temperton (1992) to be employed and the PBCGST iterative method to be used for solving the system of linear equations. (See section §13 and §14 for discussion of choice of FFT algorithm and choice of equation-solving algorithm.)

The sample parameter file specifies (via option LATTDR) that the “Lattice Dispersion Relation” of Draine and Goodman (1993) be employed to determine the dipole polarizabilities. See §11 for discussion of other options.

The sample ddscat.par file specifies that the calculations be done for a single wavelength ($`6.2832\mu \mathrm{m}`$) and a single effective radius ($`1.00\mu \mathrm{m}`$). Note that in DDSCAT.5a the “effective radius” $`a_{\mathrm{eff}}`$ is the radius of a sphere of equal volume – i.e., a sphere of volume $`Nd^3`$ , where $`d`$ is the lattice spacing and $`N`$ is the number of occupied (i.e., non-vacuum) lattice sites in the target. Thus the effective radius $`a_{\mathrm{eff}}=(3N/4\pi )^{1/3}d`$ .

The incident radiation is always assumed to propagate along the $`x`$ axis in the “Lab Frame”. The sample ddscat.par file specifies incident polarization state $`\widehat{𝐞}_{01}`$ to be along the $`y`$ axis (and consequently polarization state $`\widehat{𝐞}_{02}`$ will automatically be taken to be along the $`z`$ axis). IORTH=2 in ddscat.par calls for calculations to be carried out for both incident polarization states ($`\widehat{𝐞}_{01}`$ and $`\widehat{𝐞}_{02}`$ – see 21).

The target is assumed to have two vectors $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ embedded in it; $`\widehat{𝐚}_2`$ is perpendicular to $`\widehat{𝐚}_1`$. In the case of the 8$`\times `$6$`\times `$4 rectangular array of the sample calculation, the vector $`\widehat{𝐚}_1`$ is along the “long” axis of the target, and the vector $`\widehat{𝐚}_2`$ is along the “intermediate” axis. The target orientation in the Lab Frame is set by three angles: $`\beta `$, $`\mathrm{\Theta }`$, and $`\mathrm{\Phi }`$, defined and discussed below in §17. Briefly, the polar angles $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ specify the direction of $`\widehat{𝐚}_1`$ in the Lab Frame. The target is assumed to be rotated around $`\widehat{𝐚}_1`$ by an angle $`\beta `$. The sample ddscat.par file specifies $`\beta =0`$ and $`\mathrm{\Phi }=0`$ (see lines in ddscat.par specifying variables BETA and PHI), and calls for three values of the angle $`\mathrm{\Theta }`$ (see line in ddscat.par specifying variable THETA). DDSCAT.5a chooses $`\mathrm{\Theta }`$ values uniformly spaced in $`\mathrm{cos}\mathrm{\Theta }`$; thus, asking for three values of $`\mathrm{\Theta }`$ between 0 and $`90^{}`$ yields $`\mathrm{\Theta }=0`$, $`60^{}`$, and $`90^{}`$.

Appendix A provides a detailed description of the file ddscat.par.666 Note that the structure of ddscat.par depends on the version of DDSCAT being used. Make sure you update old parameter files before using them with DDSCAT.5a !

## 9 Running DDSCAT.5a Using the Sample ddscat.par File

To execute the program on a UNIX system (running either sh or csh), simply type

ddscat \>& ddscat.out &

which will redirect the “standard output” to the file ddscat.out, and run the calculation in the background. The sample calculation (8x6x4=192 dipole target, 3 orientations, two incident polarizations, with scattering calculated for 14 distinct scattering directions), requires 1.0 cpu seconds on a Sun Ultra-170 workstation.

## 10 Output Files

### 10.1 ASCII files

If you run DDSCAT using the command ddscat \>& ddscat.out & you will have various types of ASCII files when the computation is complete:

* a file ddscat.out;
* a file mtable;
* a file qtable;
* a file qtable2;
* files wxxryyori.avg (one, w00r00ori.avg, for the sample calculation);
* if ddscat.par specified IWRKSC=1, there will also be files wxxryykzzz.sca (3 for the sample calculation: w00r00k000.sca, w00r00k001.sca, w00r00k002.sca).

The file ddscat.out will contain any error messages generated as well as a running report on the progress of the calculation, including creation of the target dipole array. During the iterative calculations, $`Q_{\mathrm{ext}}`$, $`Q_{\mathrm{abs}}`$, and $`Q_{\mathrm{pha}}`$ are printed after each iteration; you will be able to judge the degree to which convergence has been achieved. Unless TIMEIT has been disabled, there will also be timing information.

The file mtable contains a summary of the dielectric constant used in the calculations.

The file qtable contains a summary of the orientationally-averaged values of $`Q_{\mathrm{ext}}`$, $`Q_{\mathrm{abs}}`$, $`Q_{\mathrm{sca}}`$, $`g(1)=\mathrm{cos}(\theta _s)`$, and $`Q_{bk}`$. Here $`Q_{\mathrm{ext}}`$, $`Q_{\mathrm{abs}}`$, and $`Q_{\mathrm{sca}}`$ are the extinction, absorption, and scattering cross sections divided by $`\pi a_{\mathrm{eff}}^2`$. $`Q_{bk}`$ is the differential cross section for backscattering (area per sr) divided by $`\pi a_{\mathrm{eff}}^2`$.

The file qtable2 contains a summary of the orientationally-averaged values of $`Q_{\mathrm{pha}}`$, $`Q_{\mathrm{pol}}`$, and $`Q_{\mathrm{cpol}}`$. Here $`Q_{\mathrm{pha}}`$ is the “phase shift” cross section divided by $`\pi a_{\mathrm{eff}}^2`$ (see definition in Draine 1988). $`Q_{\mathrm{pol}}`$ is the “polarization efficiency factor”, equal to the difference between $`Q_{\mathrm{ext}}`$ for the two orthogonal polarization states. We define a “circular polarization efficiency factor” $`Q_{\mathrm{cpol}}Q_{\mathrm{pol}}Q_{\mathrm{pha}}`$, since an optically-thin medium with a small twist in the alignment direction will produce circular polarization in initially unpolarized light in proportion to $`Q_{\mathrm{cpol}}`$.

For each wavelength and size, DDSCAT.5a produces a file with a name of the formwxxryyori.avg, where index xx (=00, 01, 02….) designates the wavelength and index yy (=00, 01, 02…) designates the “radius”; this file contains $`Q`$ values and scattering information averaged over however many target orientations have been specified (see §17. The file w00r00ori.avg produced by the sample calculation is provided below in Appendix B.

In addition, if ddscat.par has specified IWRKSC=1 (as for the sample calculation), DDSCAT.5a will generate files with names of the form wxxryykzzz.avg, where xx and yy are as before, and index zzz =(000,001,002…) designates the target orientation; these files contain $`Q`$ values and scattering information for each of the target orientations. The structure of each of these files is very similar to that of the wxxryyori.avg files. Because these files may not be of particular interest, and take up disk space, you may choose to set IWRKSC=0 in future work. However, it is suggested that you run the sample calculation with IWRKSC=1.

The sample ddscat.par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 target size, and averaging over 3 target orientations. Running DDSCAT.5a with the sample ddscat.par file will therefore generate files w00r00k000.sca, w00r00k001.sca, and w00r00k002.sca . To understand the information contained in one of these files, please consult Appendix C, which contains an example of the file w00r00k000.sca produced in the sample calculation.

### 10.2 Binary Option

It is possible to output an “unformatted” or “binary” file (dd.bin) with fairly complete information, including header and data sections. This is accomplished by specifying either ALLBIN or ORIBIN in ddscat.par .

Subroutine writebin.f provides an example of how this can be done. The “header” section contains dimensioning and other variables which do not change with wavelength, particle geometry, and target orientation. The header section contains data defining the particle shape, wavelengths, particle sizes, and target orientations. If ALLBIN has been specified, the “data” section contains, for each orientation, Mueller matrix results for each scattering direction. The data output is limited to actual dimensions of arrays; e.g. nscat,4,4 elements of Mueller matrix are written rather than mxscat,4,4. This is an important consideration when writing postprocessing codes.

A skeletal example of a postprocessing code was written by us and is provided in subdirectory DDA/IDL. If you do plan to use the Interactive Data Language (IDL) for postprocessing, you may consider using the netCDF binary file option which offers substantial advantages over the FORTRAN unformatted write. More information about IDL is available at `http://www.rsinc.com/idl`. Unfortunately IDL requires a license and hence is not distributed with DDSCAT.

### 10.3 Machine-Independent Binary File Option: netCDF

The “unformatted” binary file is compact, fairly complete, and (in comparison to ASCII output files) is efficiently written from FORTRAN. However, binary files are not compatible between different computer architectures if written by “unformatted” WRITE from FORTRAN. Thus, you have to process the data on the same computer architecture that was used for the DDSCAT calculations. We provide an option of using netCDF with DDSCAT. The netCDF library defines a machine-independent format for representing scientific data. Together, the interface, library, and format support the creation, access, and sharing of scientific data. For more information see http://www.unidata.ucar.edu/packages/netcdf.

Several major graphics packages (for example IDL) have adopted netCDF as a standard for data transport. For this reason, and because of strong and free support of netCDF over the network by UNIDATA, we have implemented a netCDF interface in DDSCAT. This may become the method of choice for binary file storage of output from DDSCAT.

After the initial “learning curve” netCDF offers many advantages:

* It is easy to examine the contents of the file (using tools such as ncdump).
* Binary files are portable - they can be created on a supercomputer and processed on a workstation.
* Major graphics packages now offer netCDF interfaces.
* Data input and output are an order of magnitude faster than for ASCII files.
* Output data files are compact.

The disadvantages include:

* Need to have netCDF installed on your system.
* Lack of support of complex numbers.
* Nontrivial learning curve for using netCDF.
* Lack of portability of netCDF libraries.

The public-domain netCDF software is available for many operating systems from http://www.unidata.ucar.edu/packages/netcdf . The steps necessary for enabling the netCDF capability in DDSCAT.5a are listed above in §6.4. Once these have been successfully accomplished, and you are ready to run DDSCAT to produce netCDF output, you must choose either the ALLCDF or ORICDF option in ddscat.par; ALLCDF will result in a file which will include the Mueller matrix for every wavelength, particle size, and orientation, whereas ORICDF will result in a file limited to just the orientational averages for each wavelength and target size.

## 11 Dipole Polarizabilities

Option LATTDR specifies that the “Lattice Dispersion Relation” of Draine and Goodman (1993) be employed to determine the dipole polarizabilities.; other possible choices are DRAI88 (prescription used by Draine 1988) and GOBR88 (prescription used by Goedecke & O’Brien 1988 and Hage & Greenberg 1990). In the limit $`|m|kd1`$ (where $`k=2\pi /\lambda `$ is the wave vector and $`d`$ is the lattice spacing) all three options converge to the same limit; for $`|m|kd>0.1`$ there are significant differences among them. Consult the paper by Draine & Goodman (1993) for discussion and comparison of these three prescriptions. Option LATTDR is recommended for general use, based upon the tests presented by Draine & Goodman (1993).

## 12 Dielectric Functions

In order to assign the appropriate dipole polarizabilities, DDSCAT.5a must be given the dielectric constant of the material (or materials) of which the target of interest is composed. Unless the user wishes to use the dielectric function of either solid or liquid H2O (see below), his information is supplied to DDSCAT through a table (or tables), read by subroutine DIELEC in file dielec.f, and providing either the complex refractive index $`m=n+ik`$ or complex dielectric function $`ϵ=ϵ_1+iϵ_2`$ as a function of wavelength $`\lambda `$. Since $`m=ϵ^{1/2}`$, or $`ϵ=m^2`$, the user must supply either $`m`$ or $`ϵ`$, but not both.

The table formatting is intended to be quite flexible. The first line of the table consists of text, up to 80 characters of which will be read and included in the output to identify the choice of dielectric function. (For the sample problem, it consists of simply the statement m = 1.33 + 0.01i.) The second line consists of 5 integers; either the second and third or the fourth and fifth should be zero.

* The first integer specifies which column the wavelength is stored in.
* The second integer specifies which column Re$`(m)`$ is stored in.
* The third integer specifies which column Im$`(m)`$ is stored in.
* The fourth integer specifies which column Re$`(ϵ)`$ is stored in.
* The fifth integer specifies which column Im$`(ϵ)`$ is stored in.

If the second and third integers are zeros, then DIELEC will read Re$`(ϵ)`$ and Im$`(ϵ)`$ from the file; if the fourth and fifth integers are zeros, then Re$`(m)`$ and Im$`(m)`$ will be read from the file.

The third line of the file is used for column headers, and the data begins in line 4. There must be at least 3 lines of data: even if $`ϵ`$ is required at only one wavelength, please supply two additional “dummy” wavelength entries in the table so that the interpolation apparatus will not be confused.

In the event that the user wishes to compute scattering by targets with the refractive index of either solid or liquid H2O, we have included two “built-in” dielectric functions. If H2OICE is specified on line 10 of ddscat.par, DDSCAT will use the dielectric function of H2O ice at T=250K as compiled by Warren (1984). If H2OLIQ is specified on line 10 of ddscat.par, DDSCAT will use the dielectric function for liquid H2O at T=280K using subroutine REFWAT written by Eric A. Smith. For more information see http://atol.ucsd.edu/$``$pflatau/scatlib/refr.html .

## 13 Choice of FFT Algorithm

One major change in going from DDSCAT.4b to 4c and 5a was modification of the code to permit use of the GPFA FFT algorithm developed by Dr. Clive Temperton. DDSCAT continues to offer both the Brenner code as well as the “old” Temperton code as alternative FFT implementations. The “old” Temperton code requires about 11% more memory than either the Brenner or GPFA codes (to use the “old” Temperton algorithm, DDSCAT.f must be compiled with MXMEM=1 rather than MXMEM=0).

To help persuade the user that the GPFA code is an important step forward, we provide a program TSTFFT to compare the performance of different 3-D FFT implementations. To compile, link, and run this program on a Unix system,777 Non-Unix systems: the TSTFFT Fortran source code is in the file MISC.FOR. position yourself in the DDA/src directory and type

make tstfft tstfft

Output results will be written into a file res.dat. Here is a copy of the res.dat file created when run on a Sun Ultrasparc 170:

```
 CPU time (sec) for different 3-D FFT methods
 Machine= Sun Ultrasparc-170
 parameter LVR = 64
 LVR="length of vector registers" in gpfa2f,gpfa3f,gpfa5f
==========================================================
                     Brenner      Temperton      Temperton
   NX  NY  NZ        (FOURX)        (Old)          (GPFA)
    2   2   2       0.000068       0.000065       0.000262
    3   3   3       0.000200       0.000077       0.000096
    4   4   4       0.000064       0.000086       0.000111
    5   5   5       0.000651       0.000292       0.000137
    6   6   6       0.001171       0.000170       0.000223
    8   8   8       0.000350       0.000250       0.000323
    9   9   9       0.004654       0.000300       0.000505
   10  10  10       0.005218       0.000457       0.000575
   12  12  12       0.008664       0.000814       0.000751
   15  15  15       0.023387       0.001583       0.001798
   16  16  16       0.003137       0.001359       0.002386
   18  18  18       0.042196       0.003908       0.003479
   20  20  20       0.043754       0.003881       0.004010
   24  24  24       0.075672       0.010722       0.007198
   25  25  25       0.103492       0.008142       0.011232
   27  27  27       0.160173       0.011126       0.012293
   30  30  30       0.185177       0.017111       0.014631
   32  32  32       0.042211       0.034032       0.023245
   36  36  36       0.322581       0.071013       0.027502
   40  40  40       0.379551       0.133391       0.044182
   45  45  45       0.797107       0.199537       0.069740
   48  48  48       0.668408       0.255849       0.094520
   50  50  50       0.983722       0.287208       0.115833
   54  54  54       1.530319       0.427839       0.142388
   60  60  60       1.804154       0.607189       0.181221
   64  64  64       1.005013       0.806311       0.342762
   72  72  72       3.812469       1.219364       0.362328
   75  75  75       5.364067       1.234528       0.437496
   80  80  80       4.509048       1.576118       0.563839
   81  81  81       8.017456       1.798195       0.560326
   90  90  90      10.201881       2.455399       0.779869

```

It is clear that the GPFA code is generally much faster than the Brenner FFT (by a factor of 13 for the 90$`\times `$90$`\times `$90 case) and significantly faster than the “old” Temperton FFT code (by a factor of 3 for the 90$`\times `$90$`\times `$90 case), although for some cases the differences are not large (e.g., 27$`\times `$27$`\times `$27). Since the GPFA code is memory-efficient as well, it appears to be the method of choice on scalar machines. It appears to also be best on vector machines.

The GPFA code contains a parameter LVR which is set in data statements in the routines gpfa2f, gpfa3f, and gpfa5f. LVR is supposed to be optimized to correspond to the “length of a vector register” on vector machines. As delivered, this parameter is set to 64, which is supposed to be appropriate for Crays other than the C90 (for the C90, 128 is supposed to be preferable).888 In place of “preferable” users are encouraged to read “necessary” – there is some basis for fearing that results computed on a C90 with LVR other than 128 run the risk of being incorrect! Please use LVR=128 if running on a C90! The value of LVR is not critical for scalar machines, as long as it is fairly large. We found little difference between LVR=64 and 128 on a Sparc 10/51 and on an Ultrasparc 170. You may wish to experiment with different LVR values on your computer architecture. To change LVR, you need to edit gpfa.f and change the three data statements where LVR is set.

The choice of FFT implementation is obtained by specifying one of:

* GPFAFT to use the GPFA algorithm (Temperton 1992) – this is recommended!;
* TMPRTN to obtain the “old” Temperton (1983) implementation;
* BRENNR to obtain the FOURX implementation of the Brenner (1969) algorithm;
* CONVEX to use the native FFT routine on a Convex system

## 14 Choice of Iterative Algorithm

As discussed elsewhere (e.g., Draine 1988), the problem of electromagnetic scattering of an incident wave $`𝐄_{inc}`$ by an array of $`N`$ point dipoles can be cast in the form

$$\mathrm{𝐀𝐏}=𝐄$$
(17)
where $`𝐄`$ is a $`3N`$-dimensional (complex) vector of the incident electric field $`𝐄_{inc}`$ at the $`N`$ lattice sites, $`𝐏`$ is a $`3N`$-dimensional (complex) vector of the (unknown) dipole polarizations, and $`𝐀`$ is a $`3N\times 3N`$ complex matrix.

Because $`3N`$ is a large number, direct methods for solving this system of equations for the unknown vector $`𝐏`$ are impractical, but iterative methods are useful: we begin with a guess (typically, $`𝐏=0`$) for the unknown polarization vector, and then iteratively improve the estimate for $`𝐏`$ until equation (17) is solved to some error criterion. The error tolerance may be specified as

$$\frac{|𝐀^{}\mathrm{𝐀𝐏}𝐀^{}𝐄|}{|𝐀^{}𝐄|}<h,$$
(18)
where $`𝐀^{}`$ is the Hermitian conjugate of $`𝐀`$ \[$`(A^{})_{ij}(A_{ji})^{}`$\], and $`h`$ is the error tolerance. We typically use $`h=10^5`$ in order to satisfy eq.(17) to high accuracy. The error tolerance $`h`$ can be specified by the user (see Appendix A).

A major change in going from DDSCAT.4b to 5a is the implementation of several different algorithms for iterative solution of the system of complex linear equations. DDSCAT.5a has been modified to permit solution algorithms to be treated in a fairly “modular” fashion, facilitating the testing of different algorithms. Many algorithms were compared by Flatau (1997)999 A postscript copy of this report – file cg.ps – is distributed with the DDSCAT.5a documentation.; two of them performed well and are made available to the user in DDSCAT.5a. The choice of algorithm is made by specifying one of the options:

* PBCGST – Preconditioned BiConjugate Gradient with STabilitization method from the Parallel Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins.
* PETRKP – the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), as coded in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the algorithm discussed by Draine (1988) and used in previous versions of DDSCAT.

Both methods work well. Our experience suggests that PBCGST is often faster than PETRKP, by perhaps a factor of two. We therefore recommend it, but include the PETRKP method as an alternative.

The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rdd@ukc.ac.uk) and Tim Hopkins (trh@ukc.ac.uk) is a collection of Fortran 77 routines designed to solve systems of linear equations on parallel and scalar computers using a variety of iterative methods (available at http://www.mat.ufrgs.br/pim-e.html). PIM offers a number of iterative methods, including

* Conjugate-Gradients, CG (Hestenes 1952),
* Bi-Conjugate-Gradients, BICG (Fletcher 1976),
* Conjugate-Gradients squared, CGS (Sonneveld 1989),
* the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991),
* the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992)
* the restarted generalized minimal residual, RGMRES (Saad 1986),
* the restarted generalized conjugate residual, RGCR (Eisenstat 1983),
* the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955),
* the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996),
* transpose-free quasi-minimal residual, TFQMR (Freund 1992),
* the Chebyshev acceleration, CHEBYSHEV (Young 1981).

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can be called directly via ddscat.par. It is possible (and was done) to add other options by changing the code in getfml.f . A helpful introduction to conjugate gradient methods is provided by the report “Conjugate Gradient Method Without Agonizing Pain” by Jonathan R. Shewchuk, available as a postscript file: ftp://REPORTS.ADM.CS.CMU.EDU/usr0/anon/1994/CMU-CS-94-125.ps.

## 15 Calculation of Radiative Force and Torque

In addition to solving the scattering problem for a dipole array, DDSCAT.5a can compute the three-dimensional force $`𝐅_{\mathrm{rad}}`$ and torque $`𝚪_{\mathrm{rad}}`$ exerted on this array by the incident and scattered radiation fields. This calculation is carried out, after solving the scattering problem, provided DOTORQ has been specified in ddscat.par. For each incident polarization mode, the results are given in terms of dimensionless efficiency vectors $`𝐐_{pr}`$ and $`𝐐_\mathrm{\Gamma }`$, defined by

$$𝐐_{pr}\frac{𝐅_{\mathrm{rad}}}{\pi a_{\mathrm{eff}}^2u_{\mathrm{rad}}},$$
(19)
$$𝐐_\mathrm{\Gamma }\frac{k𝚪_{\mathrm{rad}}}{\pi a_{\mathrm{eff}}^2u_{\mathrm{rad}}},$$
(20)
where $`𝐅_{\mathrm{rad}}`$ and $`𝚪_{\mathrm{rad}}`$ are the time-averaged force and torque on the dipole array, $`k=2\pi /\lambda `$ is the wavenumber in vacuo, and $`u_{\mathrm{rad}}=E_0^2/8\pi `$ is the time-averaged energy density for an incident plane wave with amplitude $`E_0\mathrm{cos}(\omega t+\varphi )`$. The radiation pressure efficiency vector can be written

$$𝐐_{pr}=Q_{\mathrm{ext}}\widehat{𝐤}Q_{\mathrm{sca}}𝐠,$$
(21)
where $`\widehat{𝐤}`$ is the direction of propagation of the incident radiation, and the vector g is the mean direction of propagation of the scattered radiation:

$$𝐠=\frac{1}{C_{\mathrm{sca}}}𝑑\mathrm{\Omega }\frac{dC_{\mathrm{sca}}(\widehat{𝐧},\widehat{𝐤})}{d\mathrm{\Omega }}\widehat{𝐧},$$
(22)
where $`d\mathrm{\Omega }`$ is the element of solid angle in scattering direction $`\widehat{𝐧}`$, and $`dC_{\mathrm{sca}}/d\mathrm{\Omega }`$ is the differential scattering cross section.

Equations for the evaluation of the radiative force and torque are derived by Draine & Weingartner (1996). It is important to note that evaluation of $`𝐐_{pr}`$ and $`𝐐_\mathrm{\Gamma }`$ involves averaging over scattering directions to evaluate the linear and angular momentum transport by the scattered wave. This evaluation requires appropriate choices of the parameters ICTHM and IPHM – see §22.

## 16 Memory Requirements

Since Fortran-77 does not allow for dynamic memory allocation, the executable has compiled into it the dimensions for a number of arrays; these array dimensions limit the size of the dipole array which the code can handle. The source code supplied to you (which can be used to run the sample calculation with 192 dipoles) is restricted to problems with targets which are maximum extent of 8 lattice spacings along the $`x`$-, $`y`$-, and $`z`$-directions (MXNX=8,MXNY=8,MXNZ=8; i.e, the target must fit within an 8$`\times `$8$`\times `$8=512 cube) and involve at most 9 different dielectric functions (MXCOMP=9). With this dimensioning, the executable requires about 1.3 MB of memory to run on an Ultrasparc system; memory requirements on other hardware and with other compilers should be similar. To run larger problems, you will need to edit the code to change PARAMETER values and recompile.

All of the dimensioning takes place in the main program DDSCAT – this should be the only routine which it is necessary to recompile. All of the dimensional parameters are set in PARAMETER statements appearing before the array declarations. You need simply edit the parameter statements. Remember, of course, that the amount of memory allocated by the code when it runs will depend upon these dimensioning parameters, so do not set them to ridiculously large values! The parameters MXNX, MXNY, MXNZ specify the maximum extent of the target in the $`x`$-, $`y`$-, or $`z`$-directions. Set the parameter MXMEM=0 or 1 depending on whether workspace required by Temperton’s FFT algorithm is to be reserved. The parameter MXCOMP specifies the maximum number of different dielectric functions which the code can handle at one time. The comment statements in the code supply all the information you should need to change these parameters.

The memory requirement for DDSCAT.5a (with the netCDF option disabled) is approximately $`(1059+0.623\mathrm{𝙼𝚇𝙽𝚇}\times \mathrm{𝙼𝚇𝙽𝚈}\times \mathrm{𝙼𝚇𝙽𝚉})`$ kbytes when compiled with MXMEM=0, or $`(1059+0.686\mathrm{𝙼𝚇𝙽𝚇}\times \mathrm{𝙼𝚇𝙽𝚈}\times \mathrm{𝙼𝚇𝙽𝚉})`$ kbytes when compiled with MXMEM=1. Thus with MXMEM=0 a 32$`\times `$32$`\times `$32 calculation requires 21.0 MBytes, while a 48$`\times `$48$`\times `$48 calculation requires 68.4 MBytes. These values are for an Ultrasparc system using the Sun FORTRAN compiler.

## 17 Target Orientation

Recall that we define a “Lab Frame” (LF) in which the incident radiation propagates in the $`+x`$ direction. In ddscat.par one specifies the first polarization state $`\widehat{𝐞}_{01}`$ (which obviously must lie in the $`y,z`$ plane in the LF); DDSCAT automatically constructs a second polarization state $`\widehat{𝐞}_{02}=\widehat{𝐱}\times \widehat{𝐞}_{01}^{}`$ orthogonal to $`\widehat{𝐞}_{01}`$ (here $`\widehat{𝐱}`$ is the unit vector in the $`+x`$ direction of the LF. For purposes of discussion we will always let unit vectors $`\widehat{𝐱}`$, $`\widehat{𝐲}`$, $`\widehat{𝐳}=\widehat{𝐱}\times \widehat{𝐲}`$ be the three coordinate axes of the LF. Users will often find it convenient to let polarization vectors $`\widehat{𝐞}_{01}=\widehat{𝐲}`$, $`\widehat{𝐞}_{02}=\widehat{𝐳}`$ (although this is not mandatory – see §21).

Recall that definition of a target involves specifying two unit vectors, $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$, which are imagined to be “frozen” into the target. We require $`\widehat{𝐚}_2`$ to be orthogonal to $`\widehat{𝐚}_1`$. Therefore we may define a “Target Frame” (TF) defined by the three unit vectors $`\widehat{𝐚}_1`$, $`\widehat{𝐚}_2`$, and $`\widehat{𝐚}_3=\widehat{𝐚}_1\times \widehat{𝐚}_2`$ .

For example, when DDSCAT creates a 8$`\times `$6$`\times `$4 rectangular solid, it fixes $`\widehat{𝐚}_1`$ to be along the longest dimension of the solid, and $`\widehat{𝐚}_2`$ to be along the next-longest dimension.

Orientation of the target relative to the incident radiation can in principle be determined two ways:

1. specifying the direction of $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ in the LF, or
2. specifying the directions of $`\widehat{𝐱}`$ (incidence direction) and $`\widehat{𝐲}`$ in the TF.

DDSCAT.5a uses method 1.: the angles $`\mathrm{\Theta }`$, $`\mathrm{\Phi }`$, and $`\beta `$ are specified in the file ddscat.par. The target is oriented such that the polar angles $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ specify the direction of $`\widehat{𝐚}_1`$ relative to the incident direction $`\widehat{𝐱}`$, where the $`\widehat{𝐱}`$,$`\widehat{𝐲}`$ plane has $`\mathrm{\Phi }=0`$. Once the direction of $`\widehat{𝐚}_1`$ is specified, the angle $`\beta `$ than specifies how the target is to rotated around the axis $`\widehat{𝐚}_1`$ to fully specify its orientation. A more extended and precise explanation follows:

### 17.1 Orientation of the Target in the Lab Frame

DDSCAT uses three angles, $`\mathrm{\Theta }`$, $`\mathrm{\Phi }`$, and $`\beta `$, to specify the directions of unit vectors $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ in the LF (see Fig. 5).

$`\mathrm{\Theta }`$ is the angle between $`\widehat{𝐚}_1`$ and $`\widehat{𝐱}`$.

When $`\mathrm{\Phi }=0`$, $`\widehat{𝐚}_1`$ will lie in the $`\widehat{𝐱},\widehat{𝐲}`$ plane. When $`\mathrm{\Phi }`$ is nonzero, it will refer to the rotation of $`\widehat{𝐚}_1`$ around $`\widehat{𝐱}`$: e.g., $`\mathrm{\Phi }=90^{}`$ puts $`\widehat{𝐚}_1`$ in the $`\widehat{𝐱},\widehat{𝐳}`$ plane.

When $`\beta =0`$, $`\widehat{𝐚}_2`$ will lie in the $`\widehat{𝐱},\widehat{𝐚}_1`$ plane, in such a way that when $`\mathrm{\Theta }=0`$ and $`\mathrm{\Phi }=0`$, $`\widehat{𝐚}_2`$ is in the $`\widehat{𝐲}`$ direction: e.g, $`\mathrm{\Theta }=90^{}`$, $`\mathrm{\Phi }=0`$, $`\beta =0`$ has $`\widehat{𝐚}_1=\widehat{𝐲}`$ and $`\widehat{𝐚}_2=\widehat{𝐱}`$. Nonzero $`\beta `$ introduces an additional rotation of $`\widehat{𝐚}_2`$ around $`\widehat{𝐚}_1`$: e.g., $`\mathrm{\Theta }=90^{}`$, $`\mathrm{\Phi }=0`$, $`\beta =90^{}`$ has $`\widehat{𝐚}_1=\widehat{𝐲}`$ and $`\widehat{𝐚}_2=\widehat{𝐳}`$.

Mathematically:

$`\widehat{𝐚}_1`$ $`=`$ $`\widehat{𝐱}\mathrm{cos}\mathrm{\Theta }+\widehat{𝐲}\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\mathrm{\Phi }+\widehat{𝐳}\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\mathrm{\Phi }`$ (23)
$`\widehat{𝐚}_2`$ $`=`$ $`\widehat{𝐱}\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\beta +\widehat{𝐲}[\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\beta \mathrm{cos}\mathrm{\Phi }\mathrm{sin}\beta \mathrm{sin}\mathrm{\Phi }]`$ (24)
$`+\widehat{𝐳}[\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\beta \mathrm{sin}\mathrm{\Phi }+\mathrm{sin}\beta \mathrm{cos}\mathrm{\Phi }]`$
$`\widehat{𝐚}_3`$ $`=`$ $`\widehat{𝐱}\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\beta \widehat{𝐲}[\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\beta \mathrm{cos}\mathrm{\Phi }+\mathrm{cos}\beta \mathrm{sin}\mathrm{\Phi }]`$ (25)
$`\widehat{𝐳}[\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\beta \mathrm{sin}\mathrm{\Phi }\mathrm{cos}\beta \mathrm{cos}\mathrm{\Phi }]`$
or, equivalently:

$`\widehat{𝐱}`$ $`=`$ $`\widehat{𝐚}_1\mathrm{cos}\mathrm{\Theta }\widehat{𝐚}_2\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\beta +\widehat{𝐚}_3\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\beta `$ (26)
$`\widehat{𝐲}`$ $`=`$ $`\widehat{𝐚}_1\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\mathrm{\Phi }+\widehat{𝐚}_2[\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\beta \mathrm{cos}\mathrm{\Phi }\mathrm{sin}\beta \mathrm{sin}\mathrm{\Phi }]`$ (27)
$`\widehat{𝐚}_3[\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\beta \mathrm{cos}\mathrm{\Phi }+\mathrm{cos}\beta \mathrm{sin}\mathrm{\Phi }]`$
$`\widehat{𝐳}`$ $`=`$ $`\widehat{𝐚}_1\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\mathrm{\Phi }+\widehat{𝐚}_2[\mathrm{cos}\mathrm{\Theta }\mathrm{cos}\beta \mathrm{sin}\mathrm{\Phi }+\mathrm{sin}\beta \mathrm{cos}\mathrm{\Phi }]`$ (28)
$`\widehat{𝐚}_3[\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\beta \mathrm{sin}\mathrm{\Phi }\mathrm{cos}\beta \mathrm{cos}\mathrm{\Phi }]`$
### 17.2 Orientation of the Incident Beam in the Target Frame

Under some circumstances, one may wish to specify the target orientation such that $`\widehat{𝐱}`$ (the direction of propagation of the radiation) and $`\widehat{𝐲}`$ (usually the first polarization direction) and $`\widehat{𝐳}`$ (= $`\widehat{𝐱}\times \widehat{𝐲})`$ refer to certain directions in the TF. Given the definitions of the LF and TF above, this is simply an exercise in coordinate transformation. For example, one might wish to have the incident radiation propagating along the (1,1,1) direction in the TF (example 14 below). Here we provide some selected examples:

1. $`\widehat{𝐱}=\widehat{𝐚}_1`$, $`\widehat{𝐲}=\widehat{𝐚}_2`$, $`\widehat{𝐳}=\widehat{𝐚}_3`$ : $`\mathrm{\Theta }=0`$, $`\mathrm{\Phi }+\beta =0`$
2. $`\widehat{𝐱}=\widehat{𝐚}_1`$, $`\widehat{𝐲}=\widehat{𝐚}_3`$, $`\widehat{𝐳}=\widehat{𝐚}_2`$ : $`\mathrm{\Theta }=0`$, $`\mathrm{\Phi }+\beta =90^{}`$
3. $`\widehat{𝐱}=\widehat{𝐚}_2`$, $`\widehat{𝐲}=\widehat{𝐚}_1`$, $`\widehat{𝐳}=\widehat{𝐚}_3`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =180^{}`$, $`\mathrm{\Phi }=0`$
4. $`\widehat{𝐱}=\widehat{𝐚}_2`$, $`\widehat{𝐲}=\widehat{𝐚}_3`$, $`\widehat{𝐳}=\widehat{𝐚}_1`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =180^{}`$, $`\mathrm{\Phi }=90^{}`$
5. $`\widehat{𝐱}=\widehat{𝐚}_3`$, $`\widehat{𝐲}=\widehat{𝐚}_1`$, $`\widehat{𝐳}=\widehat{𝐚}_2`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =90^{}`$, $`\mathrm{\Phi }=0`$
6. $`\widehat{𝐱}=\widehat{𝐚}_3`$, $`\widehat{𝐲}=\widehat{𝐚}_2`$, $`\widehat{𝐳}=\widehat{𝐚}_1`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =90^{}`$, $`\mathrm{\Phi }=90^{}`$
7. $`\widehat{𝐱}=\widehat{𝐚}_1`$, $`\widehat{𝐲}=\widehat{𝐚}_2`$, $`\widehat{𝐳}=\widehat{𝐚}_3`$ : $`\mathrm{\Theta }=180^{}`$, $`\beta +\mathrm{\Phi }=180^{}`$
8. $`\widehat{𝐱}=\widehat{𝐚}_1`$, $`\widehat{𝐲}=\widehat{𝐚}_3`$, $`\widehat{𝐳}=\widehat{𝐚}_2`$ : $`\mathrm{\Theta }=180^{}`$, $`\beta +\mathrm{\Phi }=90^{}`$
9. $`\widehat{𝐱}=\widehat{𝐚}_2`$, $`\widehat{𝐲}=\widehat{𝐚}_1`$, $`\widehat{𝐳}=\widehat{𝐚}_3`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =0`$, $`\mathrm{\Phi }=0`$
10. $`\widehat{𝐱}=\widehat{𝐚}_2`$, $`\widehat{𝐲}=\widehat{𝐚}_3`$, $`\widehat{𝐳}=\widehat{𝐚}_1`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =0`$, $`\mathrm{\Phi }=90^{}`$
11. $`\widehat{𝐱}=\widehat{𝐚}_3`$, $`\widehat{𝐲}=\widehat{𝐚}_1`$, $`\widehat{𝐳}=\widehat{𝐚}_2`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =90^{}`$, $`\mathrm{\Phi }=0`$
12. $`\widehat{𝐱}=\widehat{𝐚}_3`$, $`\widehat{𝐲}=\widehat{𝐚}_2`$, $`\widehat{𝐳}=\widehat{𝐚}_1`$ : $`\mathrm{\Theta }=90^{}`$, $`\beta =90^{}`$, $`\mathrm{\Phi }=90^{}`$
13. $`\widehat{𝐱}=(\widehat{𝐚}_1+\widehat{𝐚}_2)/\sqrt{2}`$, $`\widehat{𝐲}=\widehat{𝐚}_3`$, $`\widehat{𝐳}=(\widehat{𝐚}_1\widehat{𝐚}_2)/\sqrt{2}`$: $`\mathrm{\Theta }=45^{}`$, $`\beta =180^{}`$, $`\mathrm{\Phi }=90^{}`$
14. $`\widehat{𝐱}=(\widehat{𝐚}_1+\widehat{𝐚}_2+\widehat{𝐚}_3)/\sqrt{3}`$, $`\widehat{𝐲}=(\widehat{𝐚}_1\widehat{𝐚}_2)/\sqrt{2}`$, $`\widehat{𝐳}=(\widehat{𝐚}_1+\widehat{𝐚}_22\widehat{𝐚}_3)/\sqrt{6}`$: $`\mathrm{\Theta }=54.7356^{}`$, $`\beta =135^{}`$, $`\mathrm{\Phi }=30^{}`$.

### 17.3 Sampling in $`\mathrm{\Theta }`$, $`\mathrm{\Phi }`$, and $`\beta `$

The present version, DDSCAT.5a, chooses the angles $`\beta `$, $`\mathrm{\Theta }`$, and $`\mathrm{\Phi }`$ to sample the intervals (BETAMI,BETAMX), (THETMI,THETMX), (PHIMIN,PHIMAX), where BETAMI, BETAMX, THETMI, THETMX, PHIMIN, PHIMAX are specified in ddscat.par . The prescription for choosing the angles is to:

* uniformly sample in $`\beta `$;
* uniformly sample in $`\mathrm{\Phi }`$;
* uniformly sample in $`\mathrm{cos}\mathrm{\Theta }`$.

This prescription is appropriate for random orientation of the target, within the specified limits of $`\beta `$, $`\mathrm{\Phi }`$, and $`\mathrm{\Theta }`$.

Note that when DDSCAT.5a chooses angles it handles $`\beta `$ and $`\mathrm{\Phi }`$ differently from $`\mathrm{\Theta }`$.101010 This is a change from DDSCAT.4a!!. The range for $`\beta `$ is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if you take BETAMI=0, BETAMX=90, NBETA=2 you will get $`\beta =22.5^{}`$ and $`67.5^{}`$. Similarly, if you take PHIMIN=0, PHIMAX=180, NPHI=2 you will get $`\mathrm{\Phi }=45^{}`$ and $`135^{}`$.

Sampling in $`\mathrm{\Theta }`$ is done quite differently from sampling in $`\beta `$ and $`\mathrm{\Phi }`$. First, as already mentioned above, DDSCAT.5a samples uniformly in $`\mathrm{cos}\mathrm{\Theta }`$, not $`\mathrm{\Theta }`$. Secondly, the sampling depends on whether NTHETA is even or odd.

* If NTHETA is odd, then the values of $`\mathrm{\Theta }`$ selected include the extreme values THETMI and THETMX; thus, THETMI=0, THETMX=90, NTHETA=3 will give you $`\mathrm{\Theta }=0,60^{},90^{}`$.
* If NTHETA is even, then the range of $`\mathrm{cos}\mathrm{\Theta }`$ will be divided into NTHETA intervals, and the midpoint of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give you $`\mathrm{\Theta }=41.41^{}`$ and $`75.52^{}`$ \[$`\mathrm{cos}\mathrm{\Theta }=0.25`$ and $`0.75`$\].

The reason for this is that if odd NTHETA is specified, then the “integration” over $`\mathrm{cos}\mathrm{\Theta }`$ is performed using Simpson’s rule for greater accuracy. If even NTHETA is specified, then the integration over $`\mathrm{cos}\mathrm{\Theta }`$ is performed by simply taking the average of the results for the different $`\mathrm{\Theta }`$ values.

If averaging over orientations is desired, it is recommended that the user specify an odd value of NTHETA so that Simpson’s rule will be employed.

## 18 Orientational Averaging

DDSCAT has been constructed to facilitate the computation of orientational averages. How to go about this depends on the distribution of orientations which is applicable.

### 18.1 Randomly-Oriented Targets

For randomly-oriented targets, we wish to compute the orientational average of a quantity $`Q(\beta ,\mathrm{\Theta },\mathrm{\Phi })`$:

$$Q=\frac{1}{8\pi ^2}_0^{2\pi }𝑑\beta _1^1d\mathrm{cos}\mathrm{\Theta }_0^{2\pi }𝑑\mathrm{\Phi }Q(\beta ,\mathrm{\Theta },\mathrm{\Phi }).$$
(29)
To compute such averages, all you need to do is edit the file ddscat.par so that DDSCAT knows what ranges of the angles $`\beta `$, $`\mathrm{\Theta }`$, and $`\mathrm{\Phi }`$ are of interest. For a randomly-oriented target with no symmetry, you would need to let $`\beta `$ run from 0 to $`360^{}`$, $`\mathrm{\Theta }`$ from 0 to $`180^{}`$, and $`\mathrm{\Phi }`$ from 0 to $`360^{}`$.

For targets with symmetry, on the other hand, the ranges of $`\beta `$, $`\mathrm{\Theta }`$, and $`\mathrm{\Phi }`$ may be reduced. First of all, remember that averaging over $`\mathrm{\Phi }`$ is relatively “inexpensive”, so when in doubt average over 0 to $`360^{}`$; most of the computational “cost” is associated with the number of different values of ($`\beta `$,$`\mathrm{\Theta }`$) which are used. Consider a cube, for example, with axis $`\widehat{𝐚}_1`$ normal to one of the cube faces; for this cube $`\beta `$ need run only from 0 to $`90^{}`$, since the cube has fourfold symmetry for rotations around the axis $`\widehat{𝐚}_1`$. Furthermore, the angle $`\mathrm{\Theta }`$ need run only from 0 to $`90^{}`$, since the orientation ($`\beta `$,$`\mathrm{\Theta }`$,$`\mathrm{\Phi }`$) is indistinguishable from ($`\beta `$, $`180^{}\mathrm{\Theta }`$, $`360^{}\mathrm{\Phi }`$).

For targets with symmetry, the user is encouraged to test the significance of $`\beta `$,$`\mathrm{\Theta }`$,$`\mathrm{\Phi }`$ on targets with small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry.

It is important to remember that DDSCAT.4b handled even and odd values of NTHETA differently – see §8 above! For averaging over random orientations odd values of NTHETA are to be preferred, as this will allow use of Simpson’s rule in evaluating the “integral” over $`\mathrm{cos}\mathrm{\Theta }`$.

### 18.2 Nonrandomly-Oriented Targets

Some special cases (where the target orientation distribution is uniform for rotations around the $`x`$ axis = direction of propagation of the incident radiation), one may be able to use DDSCAT.5a with appropriate choices of input parameters. More generally, however, you will need to modify subroutine ORIENT to generate a list of NBETA values of $`\beta `$, NTHETA values of $`\mathrm{\Theta }`$, and NPHI values of $`\mathrm{\Phi }`$, plus two weighting arrays WGTA(1-NTHETA,1-NPHI) and WGTB(1-NBETA). Here WGTA gives the weights which should be attached to each ($`\mathrm{\Theta }`$,$`\mathrm{\Phi }`$) orientation, and WGTB gives the weight to be attached to each $`\beta `$ orientation. Thus each orientation of the target is to be weighted by the factor WGTA$`\times `$WGTB. For the case of random orientations, DDSCAT.5a chooses $`\mathrm{\Theta }`$ values which are uniformly spaced in $`\mathrm{cos}\mathrm{\Theta }`$, and $`\beta `$ and $`\mathrm{\Phi }`$ values which are uniformly spaced, and therefore uses uniform weights WGTB=1./NBETA When NTHETA is even, DDSCAT sets WGTA=1./(NTHETA$`\times `$NPHI) but when NTHETA is odd, DDSCAT uses Simpson’s rule when integrating over $`\mathrm{\Theta }`$ and WGTA= (1/3 or 4/3 or 2/3)/(NTHETA$`\times `$NPHI)

Note that the program structure of DDSCAT may not be ideally suited for certain highly oriented cases. If, for example, the orientation is such that for a given $`\mathrm{\Phi }`$ value only one $`\mathrm{\Theta }`$ value is possible (this situation might describe ice needles oriented with the long axis perpendicular to the vertical in the Earth’s atmosphere, illuminated by the Sun at other than the zenith) then it is foolish to consider all the combinations of $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ which the present version of DDSCAT is set up to do. We hope to improve this in a future version of DDSCAT.

## 19 Target Generation

DDSCAT contains routines to generate dipole arrays representing targets of various geometries, including spheres, ellipsoids, rectangular solids, cylinders, hexagonal prisms, tetrahedra, two touching ellipsoids, and three touching ellipsoids. The target type is determined by the variable CSHAPE (= ELLIPS, CYLNDR, etc., …) on line 7 of ddscat.par, and up to 6 parameters (SHPAR1, SHPAR2, SHPAR3, …) on line 8. The target geometry is most conveniently described in a coordinate system attached to the target which we refer to as the “Target Frame” (TF), so in this section only we will let x,y,z be coordinates in the Target Frame. Once the target is generated, the orientation of the target in the Lab Frame is accomplished as described in §17.

Target geometries currently supported include:

* ELLIPSHomogeneous, isotropic ellipsoid.
“Lengths” SHPAR1, SHPAR2, SHPAR3 in the $`x`$, $`y`$, $`z`$ directions in the TF. $`(x/\mathrm{𝚂𝙷𝙿𝙰𝚁𝟷})^2+(y/\mathrm{𝚂𝙷𝙿𝙰𝚁𝟸})^2+(z/\mathrm{𝚂𝙷𝙿𝙰𝚁𝟹})^2=d^2/4`$, where $`d`$ is the interdipole spacing.
The target axes are set to $`\widehat{𝐚}_1=(1,0,0)`$ and $`\widehat{𝐚}_2=(0,1,0)`$ in the TF.
User must set NCOMP=1 on line 9 of ddscat.par.
A homogeneous, isotropic sphere is obtained by setting SHPAR1 = SHPAR2 = SHPAR3 = diameter/$`d`$.
* ANIELLHomogeneous, anisotropic ellipsoid.
SHPAR1, SHPAR2, SHPAR3 have same meaning as for ELLIPS. Target axes $`\widehat{𝐚}_1=(1,0,0)`$ and $`\widehat{𝐚}_2=(0,1,0)`$ in the TF. It is assumed that the dielectric tensor is diagonal in the TF. User must set NCOMP=3 and provide $`xx`$, $`yy`$, $`zz`$ elements of the dielectric tensor.
* CONELLTwo concentric ellipsoids.
SHPAR1, SHPAR2, SHPAR3 specify the lengths along $`x`$, $`y`$, $`z`$ axes (in the TF) of the outer ellipsoid; SHPAR4, SHPAR5, SHPAR6 are the lengths, along the $`x`$, $`y`$, $`z`$ axes (in the TF) of the inner ellipsoid. The “core” within the inner ellipsoid is composed of isotropic material 1; the “mantle” between inner and outer ellipsoids is composed of isotropic material 2. Target axes $`\widehat{𝐚}_1=(1,0,0)`$, $`\widehat{𝐚}_2=(0,1,0)`$ in TF. User must set NCOMP=2 and provide dielectric functions for “core” and “mantle” materials.
* CYLNDRHomogeneous, isotropic cylinder.
Length/$`d`$=SHPAR1, diameter/$`d`$=SHPAR2, with cylinder axis $`=\widehat{𝐚}_1=(1,0,0)`$ and $`\widehat{𝐚}_2=(0,1,0)`$ in the TF. User must set NCOMP=1.
* UNICYLHomogeneous cylinder with unixial anisotropic dielectric tensor.
SHPAR1, SHPAR2 have same meaning as for CYLNDR. Cylinder axis $`=\widehat{𝐚}_1=(1,0,0)`$, $`\widehat{𝐚}_2=(0,1,0)`$. It is assumed that the dielectric tensor $`ϵ`$ is diagonal in the TF, with $`ϵ_{yy}=ϵ_{zz}`$. User must set NCOMP=2. Dielectric function 1 is for $`𝐄\widehat{𝐚}_1`$ (cylinder axis), dielectric function 2 is for $`𝐄\widehat{𝐚}_1`$.
* RCTNGLHomogeneous, isotropic, rectangular solid.
x, y, z lengths/$`d`$ = SHPAR1, SHPAR2, SHPAR3. Target axes $`\widehat{𝐚}_1=(1,0,0)`$ and $`\widehat{𝐚}_2=(0,1,0)`$ in the TF. User must set NCOMP=1.
* HEXGONHomogeneous, isotropic hexagonal prism.
SHPAR1=Length/$`d`$, SHPAR2 = 2$`\times `$(hexagon side)/$`d`$ = (distance between opposite vertices of hexagon)/$`d`$ = “hexagon diameter”/$`d`$. Target axis $`\widehat{𝐚}_1=(1,0,0)`$ in the TF is along the prism axis, and target axis $`\widehat{𝐚}_2=(0,1,0)`$ in the TF is normal to one of the rectangular faces of the hexagonal prism. User must set NCOMP=1.
* TETRAHHomogeneous, isotropic tetrahedron.
SHPAR1=length/$`d`$ of one edge. Orientation: one face parallel to $`y`$,$`z`$ plane (in the TF), opposite “vertex” is in $`+x`$ direction, and one edge is parallel to $`z`$ axis (in the TF). Target axes $`\widehat{𝐚}_1=(1,0,0)`$ \[emerging from one vertex\] and $`\widehat{𝐚}_2=(0,1,0)`$ \[emerging from an edge\] in the TF. User must set NCOMP=1.
* TWOSPHTwo touching homogeneous, isotropic spheroids, with distinct compositions.
First spheroid has length SHPAR1 along symmetry axis, diameter SHPAR2 perpendicular to symmetry axis. Second spheroid has length SHPAR3 along symmetry axis, diameter SHPAR4 perpendicular to symmetry axis. Contact point is on line connecting centroids. Line connecting centroids is in $`x`$ direction. Symmetry axis of first spheroid is in $`y`$ direction. Symmetry axis of second spheroid is in direction $`\widehat{𝐲}\mathrm{cos}(\mathrm{𝚂𝙷𝙿𝙰𝚁𝟻})+\widehat{𝐳}\mathrm{sin}(\mathrm{𝚂𝙷𝙿𝙰𝚁𝟻})`$, where $`\widehat{𝐲}`$ and $`\widehat{𝐳}`$ are basis vectors in TF, and SHPAR5 is in degrees. If SHPAR6=0., then target axes $`\widehat{𝐚}_1=(1,0,0)`$, $`\widehat{𝐚}_2=(0,1,0)`$. If SHPAR6=1., then axes $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ are set to principal axes with largest and 2nd largest moments of inertia assuming spheroids to be of uniform density. User must set NCOMP=2 and provide dielectric function files for each spheroid.
* TWOELLTwo touching, homogeneous, isotropic ellipsoids, with distinct compositions.
SHPAR1, SHPAR2, SHPAR3=x-length/$`d`$, y-length/$`d`$, $`z`$-length/$`d`$ of one ellipsoid. The two ellipsoids have identical shape, size, and orientation, but distinct dielectric functions. The line connecting ellipsoid centers is along the $`x`$-axis in the TF. Target axes $`\widehat{𝐚}_1=(1,0,0)`$ \[along line connecting ellipsoids\] and $`\widehat{𝐚}_2=(0,1,0)`$. User must set NCOMP=2 and provide dielectric function file names for both ellipsoids. Ellipsoids are in order of increasing x: first dielectric function is for ellipsoid with center at negative $`x`$, second dielectric function for ellipsoid with center at positive $`x`$.
* TWOAELTwo touching, homogeneous, anisotropic ellipsoids, with distinct compositions.
Geometry as for TWOELL; SHPAR1, SHPAR2, SHPAR3 have same meanings as for TWOELL. Target axes $`\widehat{𝐚}_1=(1,0,0)`$ and $`\widehat{𝐚}_2=(0,1,0)`$ in the TF. It is assumed that (for both ellipsoids) the dielectric tensor is diagonal in the TF. User must set NCOMP=6 and provide $`xx`$, $`yy`$, $`zz`$ components of dielectric tensor for first ellipsoid, and $`xx`$, $`yy`$, $`zz`$ components of dielectric tensor for second ellipsoid (ellipsoids are in order of increasing x).
* THRELLThree touching homogeneous, isotropic ellipsoids of equal size and orientation, but distinct compositions.
SHPAR1, SHPAR2, SHPAR3 have same meaning as for TWOELL. Line connecting ellipsoid centers is parallel to $`x`$ axis. Target axis $`\widehat{𝐚}_1=(1,0,0)`$ along line of ellipsoid centers, $`\widehat{𝐚}_2=(0,1,0)`$. User must set NCOMP=3 and provide (isotropic) dielectric functions for first, second, and third ellipsoid.
* THRAELThree touching homogeneous, anisotropic ellipsoids with same size and orientation but distinct dielectric tensors.
SHPAR1, SHPAR2, SHPAR3 have same meanings as for THRELL. Target axis $`\widehat{𝐚}_1=(1,0,0)`$ along line of ellipsoid centers, $`\widehat{𝐚}_2=(0,1,0)`$. It is assumed that dielectric tensors are all diagonal in the TF. User must set NCOMP=9 and provide $`xx`$, $`yy`$, $`zz`$ elements of dielectric tensor for first ellipsoid, $`xx`$, $`yy`$, $`zz`$ elements for second ellipsoid, and $`xx`$, $`yy`$, $`zz`$ elements for third ellipsoid (ellipsoids are in order of increasing $`x`$).
* BLOCKSHomogeneous target constructed from cubic “blocks”.
Number and location of blocks are specified in separate file blocks.par with following structure:
one line of comments (may be blank)
PRIN (= 0 or 1 – see below)
N (= number of blocks)
B (= width/$`d`$ of one block)
$`x`$ $`y`$ $`z`$ (= position of 1st block in units of B$`d`$)
$`x`$ $`y`$ $`z`$ (= position of 2nd block)
…
$`x`$ $`y`$ $`z`$ = position of Nth block
If PRIN=0, then $`\widehat{𝐚}_1=(1,0,0)`$, $`\widehat{𝐚}_2=(0,1,0)`$. If PRIN=1, then $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ are set to principal axes with largest and second largest moments of inertia, assuming target to be of uniform density. User must set NCOMP=1.
* DW199613 block target used by Draine & Weingartner (1996).
Single, isotropic material. Target geometry was used in study by Draine & Weingartner (1996) of radiative torques on irregular grains. $`\widehat{𝐚}_1`$ and $`\widehat{𝐚}_2`$ are principal axes with largest and second-largest moments of inertia. User must set NCOMP=1. Target size is controlled by shape parameter SHPAR(1) = width of one block in lattice units.
* NSPHERMultisphere target consisting of the union of $`N`$ spheres of single isotropic material.
Spheres may overlap if desired. The relative locations and sizes of these spheres are defined in an external file, whose name (enclosed in quotes) is passed through ddscat.par. The length of the file name should not exceed 13 characters. The pertinent line in ddscat.par should read
SHPAR(1) SHPAR(2) ‘filename’ (quotes must be used)
where SHPAR(1) = target diameter in $`x`$ direction (in Target Frame) in units of $`d`$
SHPAR(2)= 0 to have $`a_1=(1,0,0)`$, $`a_2=(0,1,0)`$ in Target Frame.
SHPAR(2)= 1 to use principal axes of moment of inertia tensor for $`a_1`$ (largest $`I`$) and $`a_2`$ (intermediate $`I`$).
filename is the name of the file specifying the locations and relative sizes of the spheres.
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter SHPAR(1), which is the extent of the multisphere target in the $`x`$-direction, in units of the lattice spacing $`d`$. The file ‘filename’ should have the following structure:
$`N`$ (= number of spheres)
one line of comments (may be blank)
$`x_1`$ $`y_1`$ $`z_1`$ $`a_1`$ (arb. units)
$`x_2`$ $`y_2`$ $`z_2`$ $`a_2`$ (arb. units)
…
$`x_N`$ $`y_N`$ $`z_N`$ $`a_N`$ (arb. units)
where $`x_j`$, $`y_j`$, $`z_j`$ are the coordinates of the center, and $`a_j`$ is the radius of sphere $`j`$.
Note that $`x_j`$, $`y_j`$, $`z_j`$, $`a_j`$ ($`j=1,\mathrm{},N`$) establish only the shape of the $`N`$sphere target. For instance, a target consisting of two touching spheres with the line between centers parallel to the $`x`$ axis could equally well be described by lines 3 and 4 being
0 0 0 0.5
1 0 0 0.5
or
0 0 0 1
2 0 0 1
The actual size (in physical units) is set by the value of $`a_{\mathrm{eff}}`$ specified in ddscat.par, where, as always, $`a_{\mathrm{eff}}(3V/4\pi )^{1/3}`$, where $`V`$ is the volume of material in the target.
User must set NCOMP=1.
* PRISM3Triangular prism of homogeneous, isotropic material.
SHPAR1, SHPAR2, SHPAR3, SHPAR4 $`=a/d`$, $`b/a`$, $`c/a`$, $`L/a`$
The triangular cross section has sides of width $`a`$, $`b`$, $`c`$. $`L`$ is the length of the prism. $`d`$ is the lattice spacing. The triangular cross-section has interior angles $`\alpha `$, $`\beta `$, $`\gamma `$ (opposite sides $`a`$, $`b`$, $`c`$) given by $`\mathrm{cos}\alpha =(b^2+c^2a^2)/2bc`$, $`\mathrm{cos}\beta =(a^2+c^2b^2)/2ac`$, $`\mathrm{cos}\gamma =(a^2+b^2c^2)/2ab`$. In the Target Frame, the prism axis is in the $`\widehat{𝐱}`$ direction, the normal to the rectangular face of width $`a`$ is (0,1,0), the normal to the rectangular face of width $`b`$ is $`(0,\mathrm{cos}\gamma ,\mathrm{sin}\gamma )`$, and the normal to the rectangular face of width $`c`$ is $`(0,\mathrm{cos}\gamma ,\mathrm{sin}\gamma )`$.
* FRMFILList of dipole locations and “compositions” obtained from a file.
This option causes DDSCAT to read the target geometry information from a file shape.dat instead of automatically generating one of the geometries listed above. The shape.dat file is read by routine REASHP (file reashp.f). The user can customize REASHP as needed to conform to the manner in which the target geometry is stored in file shape.dat. However, as supplied, REASHP expects the file shape.dat to have the following structure:

+ one line containing a description; the first 67 characters will be read and printed in various output statements.
    + N = number of dipoles in target
    + $`a_{1x}`$ $`a_{1y}`$ $`a_{1z}`$ = x,y,z components of $`𝐚_\mathrm{𝟏}`$
    + $`a_{2x}`$ $`a_{2y}`$ $`a_{2z}`$ = x,y,z components of $`𝐚_\mathrm{𝟐}`$
    + (line containing comments)
    + $`dummy`$ IXYZ(1,1) IXYZ(1,2) IXYZ(1,3) ICOMP(1,1) ICOMP(1,2) ICOMP(1,3)
    + $`dummy`$ IXYZ(2,1) IXYZ(2,2) IXYZ(2,3) ICOMP(2,1) ICOMP(2,2) ICOMP(2,3)
    + $`dummy`$ IXYZ(3,1) IXYZ(3,2) IXYZ(3,3) ICOMP(3,1) ICOMP(3,2) ICOMP(3,3)
    +
    + $`dummy`$ IXYZ(J,1) IXYZ(J,2) IXYZ(J,3) ICOMP(J,1) ICOMP(J,2) ICOMP(J,3)
    +
    + $`dummy`$ IXYZ(N,1) IXYZ(N,2) IXYZ(N,3) ICOMP(N,1) ICOMP(N,2) ICOMP(N,3)

where the number $`dummy`$ is ignored, IXYZ(J,1-3) = $`x/d`$, $`y/d`$, $`z/d`$ (TF) for $`J`$-th dipole, and ICOMP(J,1-3) = composition identifier to specify dielectric function appropriate for $`x,y,z`$ directions at location of $`J`$-th dipole. For the common case of a single isotropic material, ICOMP(J,1-3) = 1 1 1 .

The user should be able to easily modify these routines, or write new routines, to generate targets with other geometries. The user should first examine the routine target.f and modify it to call any new target generation routines desired. Alternatively, targets may be generated separately, and the target description (locations of dipoles and “composition” corresponding to x,y,z dielectric properties at each dipole site) read in from a file by invoking the option FRMFIL in ddscat.f.

It is often desirable to be able to run the target generation routines without running the entire DDSCAT code. We have therefore provided a program CALLTARGET which allows the user to generate targets interactively; to create this executable just type111111 Non-Unix sites: The source code for CALLTARGET is in the file MISC.FOR. You must compile and link this to ERRMSG, GETSET, LAPACKBLAS, LAPACKSUBS, PRINAXIS, REASHP, TAR2EL, TAR2SP, TAR3EL, TARBLOCKS, TARCEL, TARCYL, TARELL, TARGET, TARHEX, TARNSP, TARREC, TARTET, and WRIMSG. The source code for ERRMSG is in SRC2.FOR, GETSET is in SRC5.FOR, and the rest are in SRC8.FOR and SRC9.FOR .
make calltarget . The program calltarget is to be run interactively; the prompts are self-explanatory. You may need to edit the code to change the device number IDVOUT as for DDSCAT (see §6.1 above).

After running, calltarget will leave behind an ASCII file target.out which is a list of the occupied lattice sites in the last target generated. The format of target.out is the same as the format of the shape.dat files read if option FRMFIL is used (see above). Therefore you can simply
mv target.out shape.dat
and then use DDSCAT with the option FRMFIL in order to input a target shape generated by CALLTARGET.

## 20 Scattering Directions

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix are reported in the output files wxxryykzzz.sca . The scattering direction is specified through angles $`\theta _s`$ and $`\varphi _s`$ (not to be confused with the angles $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ which specify the orientation of the target relative to the incident radiation!).

The scattering angle $`\theta _s`$ is simply the angle between the incident beam (along direction $`\widehat{𝐱}`$) and the scattered beam ($`\theta _s=0`$ for forward scattering, $`\theta _s=180^{}`$ for backscattering).

The scattering angle $`\varphi _s`$ specifies the orientation of the “scattering plane” relative to the $`\widehat{𝐱},\widehat{𝐲}`$ plane in the Lab Frame. When $`\varphi _s=0`$ the scattering plane is assumed to coincide with the $`\widehat{𝐱},\widehat{𝐲}`$ plane. When $`\varphi _s=90^{}`$ the scattering plane is assumed to coincide with the $`\widehat{𝐱},\widehat{𝐳}`$ plane. Within the scattering plane the scattering directions are specified by $`0\theta _s180^{}`$.

Scattering directions for which the scattering properties are to be calculated are set in the file ddscat.par by specifying one or more scattering planes (determined by the value of $`\varphi _s`$) and for each scattering plane, the number and range of $`\theta _s`$ values. The only limitation is that the number of scattering directions not exceed the parameter MXSCA in DDSCAT.f (in the code as distributed it is set to MXSCA=1000).

## 21 Incident Polarization State

Recall that the “Lab Frame” is defined such that the incident radiation is propagating along the $`\widehat{𝐱}`$ axis. DDSCAT.5a allows the user to specify a general elliptical polarization for the incident radiation, by specifying the (complex) polarization vector $`\widehat{𝐞}_{01}`$. The orthonormal polarization state $`\widehat{𝐞}_{02}=\widehat{𝐱}\times \widehat{𝐞}_{01}^{}`$ is generated automatically if ddscat.par specifies IORTH=2.

For incident linear polarization, one can simply set $`\widehat{𝐞}_{01}=\widehat{𝐲}`$ by specifying (0,0) (1,0) (0,0) in ddscat.par; then $`\widehat{𝐞}_{02}=\widehat{𝐳}`$ For polarization mode $`\widehat{𝐞}_{01}`$ to correspond to right-handed circular polarization, set $`\widehat{𝐞}_{01}=(\widehat{𝐲}+i\widehat{𝐳})/\sqrt{2}`$ by specifying (0,0) (1,0) (0,1) in ddscat.par (DDSCAT.5a automatically takes care of the normalization of $`\widehat{𝐞}_{01}`$); then $`\widehat{𝐞}_{02}=(i\widehat{𝐲}+\widehat{𝐳})/\sqrt{2}`$, corresponding to left-handed circular polarization.

## 22 Averaging over Scattering: $`g(1)=\mathrm{cos}\theta _s`$, etc.

DDSCAT automatically carries out numerical integration of various scattering properties. In particular, it calculates the mean value $`g(1)=\mathrm{cos}\theta _s`$ for the scattered intensity for each incident polarization state. This is accomplished by evaluating the scattered intensity for ICTHM different values of $`\theta _s`$ (including $`\theta _s=0`$ and $`\theta _s=\pi `$), and taking a weighted sum. For each value of $`\theta _s`$ except $`0`$ and $`\pi `$, the scattering intensity is evaluated for IPHM different values of the scattering angle $`\varphi _s`$. The angular integration over $`\theta _s`$ is accomplished using Simpson’s rule (with uniform intervals in $`\mathrm{cos}\theta _s`$), so ICTHM should be an odd number. The angular integration over $`\varphi _s`$ is accomplished by sampling uniformly in $`\varphi _s`$ with uniform weighting; IPHM can be either even or odd.

The following quantities are evaluated by this angular integration:

* $`𝐠=\mathrm{cos}\theta _s\widehat{𝐱}+\mathrm{sin}\theta _s\mathrm{cos}\varphi _s\widehat{𝐲}+\mathrm{sin}\theta _s\mathrm{sin}\varphi _s\widehat{𝐳}`$ (see §15);
* $`𝐐_\mathrm{\Gamma }`$ (see §15).

It is important that the user recognize that accurate evaluation of these angular averages requires adequate sampling over scattering angles. For small values of the size parameter $`x=2\pi a_{\mathrm{eff}}/\lambda `$, the angular distribution of scattered radiation has a dipolar character and the sampling in $`\theta _s`$ and $`\varphi _s`$ does not need to be very fine, so ICTHM and IPHM need not be large. For larger values of the size parameter $`x`$, however, higher multipoles in the scattered radiation field become important, and finer sampling in $`\theta _s`$ and $`\varphi _s`$ is required. We do not have any foolproof prescription to offer, since the scattering pattern will depend upon the target geometry and dielectric constant in addition to overall size parameter. However, as a very rough guide, we suggest that the user specify values of ICTHM and IPHM satisfying

$`\mathrm{𝙸𝙲𝚃𝙷𝙼}`$ $`>`$ $`5(1+x),`$ (30)
$`\mathrm{𝙸𝙿𝙷𝙼}`$ $`>`$ $`2(1+x).`$ (31)
The sample ddscat.par file supplied has $`\mathrm{𝙸𝙲𝚃𝙷𝙼}=\mathrm{𝟹𝟹}`$ and $`\mathrm{𝙸𝙿𝙷𝙼}=12`$; the above criteria would suggest that this would be suitable for $`x<5`$.

The cpu time required for evaluation of these angular averages is proportional to $`[2+\mathrm{𝙸𝙿𝙷𝙼}(\mathrm{𝙸𝙲𝚃𝙷𝙼}2)]`$. Since the computational time spent in evaluating these angular integrals can be a significant part of the total, it is important to choose values of ICTHM and IPHM which will provide a suitable balance between accuracy (in this part of the overall calculation) and cpu time.

Within one scattering plane, the scattered intensity tends to have approximately $`(1+x)`$ peaks for $`0\theta _s\pi `$, so that the above prescription for ICTHM would have at least 5 sampling points per maximum. The angular distribution over $`\varphi _s`$ is usually not as structured as that over $`\theta _s`$ so we suggest that IPHM need not be as large as ICTHM. We have refrained from “hard-wiring” the values of ICTHM and IPHM because we are not confident of the reliability of the recommended criteria (30,31) – it is up to the user to specify appropriate values of ICTHM and IPHM according to the requirements of the problem being addressed.

## 23 Mueller Matrix for Scattering in Selected Directions

### 23.1 Two Orthogonal Incident Polarizations (IORTH=2)

DDSCAT.5a internally computes the scattering properties of the dipole array in terms of a complex scattering matrix $`f_{ml}(\theta _s,\varphi _s)`$ (Draine 1988), where index $`l=1,2`$ denotes the incident polarization states, $`m=1,2`$ denotes the scattered polarization state, and $`\theta _s`$,$`\varphi _s`$ specify the scattering direction. Normally DDSCAT is used with IORTH=2 in ddscat.par, so that the scattering problem will be solved for both incident polarization states ($`l=1`$ and 2); in this subsection it will be assumed that this is the case.

Incident polarization states $`l=1,2`$ correspond to polarization states $`\widehat{𝐞}_{01}`$, $`\widehat{𝐞}_{02}`$; recall that polarization state $`\widehat{𝐞}_{01}`$ is user-specified, and $`\widehat{𝐞}_{02}=\widehat{𝐱}\times \widehat{𝐞}_{01}^{}`$. Scattered polarization state $`m=1`$ corresponds to linear polarization of the scattered wave parallel to the scattering plane ($`\widehat{𝐞}_1=\widehat{𝐞}_s=\widehat{\theta }_s`$) and $`m=2`$ corresponds to linear polarization perpendicular to the scattering plane (in the $`+\widehat{\varphi }_s`$ direction). The scattering matrix $`f_{ml}`$ was defined (Draine 1988) so that the scattered electric field $`𝐄_s`$ is related to the incident electric field $`𝐄_i(0)`$ at the origin (where the target is assumed to be located) by

$$\left(\begin{array}{c}𝐄_s\widehat{\theta }_s\\ 𝐄_s\widehat{\varphi }_s\end{array}\right)=\frac{\mathrm{exp}(i𝐤𝐫)}{kr}\left(\begin{array}{cc}f_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right)\left(\begin{array}{c}𝐄_i(0)\widehat{𝐞}_{01}^{}\\ 𝐄_i(0)\widehat{𝐞}_{02}^{}\end{array}\right).$$
(32)
The 2$`\times `$2 complex amplitude scattering matrix (with elements $`S_1`$, $`S_2`$, $`S_3`$, and $`S_4`$) is defined so that (see Bohren & Huffman 1983)

$$\left(\begin{array}{c}𝐄_s\widehat{\theta }_s\\ 𝐄_s\widehat{\varphi }_s\end{array}\right)=\frac{\mathrm{exp}(i𝐤𝐫)}{ikr}\left(\begin{array}{cc}S_2& S_3\\ S_4& S_1\end{array}\right)\left(\begin{array}{c}𝐄_i(0)\widehat{𝐞}_i\\ 𝐄_i(0)\widehat{𝐞}_i\end{array}\right),$$
(33)
where $`\widehat{𝐞}_i`$, $`\widehat{𝐞}_i`$ are (real) unit vectors for incident polarization parallel and perpendicular to the scattering plane (with the customary definition of $`\widehat{𝐞}_i=\widehat{𝐞}_i\times \widehat{𝐱}`$).

From (32,33) we may write

$$\left(\begin{array}{cc}S_2& S_3\\ S_4& S_1\end{array}\right)\left(\begin{array}{c}𝐄_i(0)\widehat{𝐞}_i\\ 𝐄_i(0)\widehat{𝐞}_i\end{array}\right)=i\left(\begin{array}{cc}f_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right)\left(\begin{array}{c}𝐄_i(0)\widehat{𝐞}_{01}^{}\\ 𝐄_i(0)\widehat{𝐞}_{02}^{}\end{array}\right)$$
(34)
Let

$`a`$ $``$ $`\widehat{𝐞}_{01}^{}\widehat{𝐲},`$ (35)
$`b`$ $``$ $`\widehat{𝐞}_{01}^{}\widehat{𝐳},`$ (36)
$`c`$ $``$ $`\widehat{𝐞}_{02}^{}\widehat{𝐲},`$ (37)
$`d`$ $``$ $`\widehat{𝐞}_{02}^{}\widehat{𝐳}.`$ (38)
Note that since $`\widehat{𝐞}_{01},\widehat{𝐞}_{02}`$ could be complex (i.e., elliptical polarization), the quantities $`a,b,c,d`$ are complex. Then

$$\left(\begin{array}{c}\widehat{𝐞}_{01}^{}\\ \widehat{𝐞}_{02}^{}\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}\widehat{𝐲}\\ \widehat{𝐳}\end{array}\right)$$
(39)
and eq. (34) can be written

$$\left(\begin{array}{cc}S_2& S_3\\ S_4& S_1\end{array}\right)\left(\begin{array}{c}𝐄_i(0)\widehat{𝐞}_i\\ 𝐄_i(0)\widehat{𝐞}_i\end{array}\right)=i\left(\begin{array}{cc}f_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right)\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}𝐄_i(0)\widehat{𝐲}\\ 𝐄_i(0)\widehat{𝐳}\end{array}\right)$$
(40)
The incident polarization states $`\widehat{𝐞}_i`$ and $`\widehat{𝐞}_i`$ are related to $`\widehat{𝐲}`$, $`\widehat{𝐳}`$ by

$$\left(\begin{array}{c}\widehat{𝐲}\\ \widehat{𝐳}\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\varphi _s& \mathrm{sin}\varphi _s\\ \mathrm{sin}\varphi _s& \mathrm{cos}\varphi _s\end{array}\right)\left(\begin{array}{c}\widehat{𝐞}_i\\ \widehat{𝐞}_i\end{array}\right);$$
(41)
substituting (41) into (40) we obtain

$$(\begin{array}{cc}S_2& S_3\\ S_4& S_1\end{array})(\begin{array}{c}𝐄_i(0)\widehat{𝐞}_i\\ 𝐄_i(0)\widehat{𝐞}_i\end{array})=i(\begin{array}{cc}f_{11}& f_{12}\\ f_{21}& f_{22}\end{array})(\begin{array}{cc}a& b\\ c& d\end{array})(\begin{array}{cc}\mathrm{cos}\varphi _s& \mathrm{sin}\varphi _s\\ \mathrm{sin}\varphi _s& \mathrm{cos}\varphi _s\end{array})(\begin{array}{c}𝐄_i(0)\widehat{𝐲}\\ 𝐄_i(0)\widehat{𝐳}\end{array})$$
(42)
Eq. (42) must be true for all $`𝐄_i(0)`$, so we obtain an expression for the complex scattering amplitude matrix in terms of the $`f_{ml}`$:

$$\left(\begin{array}{cc}S_2& S_3\\ S_4& S_1\end{array}\right)=i\left(\begin{array}{cc}f_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right)\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{cc}\mathrm{cos}\varphi _s& \mathrm{sin}\varphi _s\\ \mathrm{sin}\varphi _s& \mathrm{cos}\varphi _s\end{array}\right).$$
(43)
This provides the 4 equations used in subroutine GETMUELLER to compute the amplitude scattering matrix elements:

$`S_1`$ $`=`$ $`i\left[f_{21}(b\mathrm{cos}\varphi _sa\mathrm{sin}\varphi _s)+f_{22}(d\mathrm{cos}\varphi _sc\mathrm{cos}\varphi _s)\right],`$ (44)
$`S_2`$ $`=`$ $`i\left[f_{11}(a\mathrm{cos}\varphi _s+b\mathrm{sin}\varphi _s)+f_{12}(c\mathrm{cos}\varphi _s+d\mathrm{sin}\varphi _s)\right],`$ (45)
$`S_3`$ $`=`$ $`i\left[f_{11}(b\mathrm{cos}\varphi _sa\mathrm{sin}\varphi _s)+f_{12}(d\mathrm{cos}\varphi _sc\mathrm{sin}\varphi _s)\right],`$ (46)
$`S_4`$ $`=`$ $`i\left[f_{21}(a\mathrm{cos}\varphi _s+b\mathrm{sin}\varphi _s)+f_{22}(c\mathrm{cos}\varphi _s+d\mathrm{sin}\varphi _s)\right].`$ (47)
It is both convenient and customary to characterize both incident and scattered radiation by 4 “Stokes parameters” – the elements of the “Stokes vector”. There are different conventions in the literature; we adhere to the definitions of the Stokes vector ($`I`$,$`Q`$,$`U`$,$`V`$) adopted in the excellent treatise by Bohren & Huffman (1983), to which the reader is referred for further detail. Here are some examples of Stokes vectors $`(I,Q,U,V)`$:

* $`(1,0,0,0)I`$ : unpolarized light.
* $`(1,1,0,0)I`$ : 100% linearly polarized with $`𝐄`$ parallel to the scattering plane;
* $`(1,1,0,0)I`$ : 100% linearly polarized with $`𝐄`$ perpendicular to the scattering plane;
* $`(1,0,1,0)I`$ : 100% linearly polarized with $`𝐄`$ at +45 relative to the scattering plane;
* $`(1,0,1,0)I`$ : 100% linearly polarized with $`𝐄`$ at -45 relative to the scattering plane;
* $`(1,0,0,1)I`$ : 100% right circular polarization (i.e., negative helicity);
* $`(1,0,0,1)I`$ : 100% left circular polarization (i.e., positive helicity).

It is convenient to describe the scattering properties in terms of the $`4\times 4`$ Mueller matrix relating the Stokes parameters $`(I_i,Q_i,U_i,V_i)`$ and $`(I_s,Q_s,U_s,V_s)`$ of the incident and scattered radiation:

$$\left(\begin{array}{c}I_s\\ Q_s\\ U_s\\ V_s\end{array}\right)=\frac{1}{k^2r^2}\left(\begin{array}{cccc}S_{11}& S_{12}& S_{13}& S_{14}\\ S_{21}& S_{22}& S_{23}& S_{24}\\ S_{31}& S_{32}& S_{33}& S_{34}\\ S_{41}& S_{42}& S_{43}& S_{44}\end{array}\right)\left(\begin{array}{c}I_i\\ Q_i\\ U_i\\ V_i\end{array}\right).$$
(48)
This $`4\times 4`$ Mueller matrix for scattering by a single particle is also referred to as the “scattering matrix”, and sometimes as the “phase matrix”.

Once the amplitude scattering matrix elements $`S_1`$, $`S_2`$, $`S_3`$, and $`S_4`$ have been obtained, the Mueller matrix elements can be computed (Bohren & Huffman 1983):

$`S_{11}`$ $`=`$ $`\left(|S_1|^2+|S_2|^2+|S_3|^2+|S_4|^2\right)/2,`$
$`S_{12}`$ $`=`$ $`\left(|S_2|^2|S_1|^2+|S_4|^2|S_3|^2\right)/2,`$
$`S_{13}`$ $`=`$ $`\mathrm{Re}\left(S_2S_3^{}+S_1S_4^{}\right),`$
$`S_{14}`$ $`=`$ $`\mathrm{Im}\left(S_2S_3^{}S_1S_4^{}\right),`$
$`S_{21}`$ $`=`$ $`\left(|S_2|^2|S_1|^2+|S_3|^2|S_4|^2\right)/2,`$
$`S_{22}`$ $`=`$ $`\left(|S_1|^2+|S_2|^2|S_3|^2|S_4|^2\right)/2,`$
$`S_{23}`$ $`=`$ $`\mathrm{Re}\left(S_2S_3^{}S_1S_4^{}\right),`$
$`S_{24}`$ $`=`$ $`\mathrm{Im}\left(S_2S_3^{}+S_1S_4^{}\right),`$
$`S_{31}`$ $`=`$ $`\mathrm{Re}\left(S_2S_4^{}+S_1S_3^{}\right),`$
$`S_{32}`$ $`=`$ $`\mathrm{Re}\left(S_2S_4^{}S_1S_3^{}\right),`$
$`S_{33}`$ $`=`$ $`\mathrm{Re}\left(S_1S_2^{}+S_3S_4^{}\right),`$
$`S_{34}`$ $`=`$ $`\mathrm{Im}\left(S_2S_1^{}+S_4S_3^{}\right),`$
$`S_{41}`$ $`=`$ $`\mathrm{Im}\left(S_4S_2^{}+S_1S_3^{}\right),`$
$`S_{42}`$ $`=`$ $`\mathrm{Im}\left(S_4S_2^{}S_1S_3^{}\right),`$
$`S_{43}`$ $`=`$ $`\mathrm{Im}\left(S_1S_2^{}S_3S_4^{}\right),`$
$`S_{44}`$ $`=`$ $`\mathrm{Re}\left(S_1S_2^{}S_3S_4^{}\right).`$ (49)
These matrix elements are computed in DDSCAT and passed to subroutine WRITESCA which handles output of scattering properties. As delivered, WRITESCA writes out 6 selected elements: $`S_{11}`$, $`S_{21}`$, $`S_{31}`$, $`S_{41}`$ (these 4 elements describe the intensity and polarization state for scattering of unpolarized incident radiation), $`S_{12}`$, and $`S_{13}`$. In addition, WRITESCA writes out the linear polarization $`P`$ of the scattered light for incident unpolarized light:

$$P=\frac{(S_{21}^2+S_{31}^2)^{1/2}}{S_{11}}.$$
(50)
Of course, other elements $`S_{ij}`$ may be of interest. It is relatively straightforward for the user to modify subroutine WRITESCA to write out whatever elements of the Mueller matrix (or the scattering amplitude matrix) are desired.

### 23.2 One Incident Polarization State Only (IORTH=1)

In some cases it may be desirable to limit the calculations to a single incident polarization state – for example, when each solution is very time-consuming, and the target is known to have some symmetry so that solving for a single incident polarization state may be sufficient for the required purpose. In this case, set IORTH=1 in ddscat.par.

When IORTH=1, only $`f_{11}`$ and $`f_{21}`$ are available; hence, DDSCAT cannot automatically generate the Mueller matrix elements. In this case, the output routine WRITESCA writes out the quantities $`|f_{11}|^2`$, $`|f_{21}|^2`$, $`\mathrm{Re}(f_{11}f_{21}^{})`$, and $`\mathrm{Im}(f_{11}f_{21}^{})`$ for each of the scattering directions.

Note, however, that if IPHI is greater than 1, DDSCAT will automatically set IORTH=2 even if ddscat.par specified IORTH=1: this is because when more than one value of the target orientation angle $`\mathrm{\Phi }`$ is required, there is no additional “cost” to solve the scattering problem for the second incident polarization state, since when solutions are available for two orthogonal states for some particular target orientation, the solution may be obtained for another target orientation differing only in the value of $`\mathrm{\Phi }`$ by appropriate linear combinations of these solutions. Hence we may as well solve the “complete” scattering problem so that we can compute the complete Mueller matrix.

## 24 Graphics and Postprocessing

### 24.1 IDL

At present, we do not offer a comprehensive package for DDSCAT data postprocessing and graphical display in IDL. However, there are several developments worth mentioning: First, we offer several output capabilities from within DDSCAT: ASCII (see §10.1), FORTRAN unformatted binary (see §10.2), and netCDF portable binary (see §10.3). Second, we offer several skeleton IDL utilities:

* bhmie.pro is our translation to IDL of the popular Bohren-Huffman code which calculates efficiencies for spherical particles using Mie theory.
* readbin.pro reads the FORTRAN unformatted binary file written by routine writebin.f. The variables are stored in a common block.
* readnet.pro reads NetCDF portable binary file and should be the method of choice for IDL users. It offers random data access.
* mie.pro is an example of an interface to binary files and the Bohren-Huffman code. It plots a comparison of DDSCAT results with scattering by equivalent radius spheres.

At present the IDL code is experimental.

## 25 Miscellanea

Additional source code, refractive index files, etc., contributed by users will be located in the directory DDA/misc. These routines and files should be considered to be not supported by Draine and Flatau – caveat receptor! These routines and files should be accompanied by enough information (e.g., comments in source code) to explain their use.

## 26 Finale

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure of the ddscat.par file is intended to be simple and suggestive so that, after reading the above notes once, the user may not have to refer to them again.

The file rel\_notes in DDA/doc lists known bugs in DDSCAT.5a . An up-to-date version will be maintained in the directory where the publicly available source code is located (see the README file). Users are encouraged to provide B. T. Draine (draine@astro.princeton.edu) with their email address; bug reports, and any new releases of DDSCAT, will be made known to those who do!

P. J. Flatau maintains the WWW page “SCATTERLIB - Light Scattering Codes Library” with URL http://atol.ucsd.edu/$``$pflatau. The SCATTERLIB Internet site is a library of light scattering codes. Emphasis is on providing source codes (mostly FORTRAN). However, other information related to scattering on spherical and non-spherical particles is collected: an extensive list of references to light scattering methods, refractive index, etc. This URL page contains a section on the discrete dipole approximation.

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. If you wish to cite this User Guide, we suggest the following citation:

Draine, B.T., & Flatau, P.J. 2000, “User Guide for the Discrete Dipole Approximation Code DDSCAT (Version 5a10)”, http://xxx.lanl.gov/abs/astro-ph/0008151v2

Finally, the authors have one special request: We would very much appreciate preprints and (especially!) reprints of any papers which make use of DDSCAT!

## 27 Acknowledgments

The routine ESELF making use of the FFT was originally written by Jeremy Goodman, Princeton University Observatory. The FFT routine FOURX is based on a FFT routine written by Norman Brenner (Brenner 1969). The routine REFICE was written by Steven B. Warren, based on Warren (1984). The routine REFWAT was written by Eric A. Smith. The GPFAPACK package was written by Clive Temperton, and generously made available by him for use with DDSCAT. Subroutine CXFFT3 is based on a FFT routine written by Clive Temperton (Temperton 1983). We make use of routines from the LAPACK package (Anderson et al. 1995), the result of work by Jack Dongarra and others at the Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University. We are indebted to all of these authors for making their code available.

We wish also to acknowledge bug reports and suggestions from DDSCAT users, including Henrietta Lemke, Timo Nousianen, and Mike Wolff.

Development of DDSCAT was supported in part by National Science Foundation grants AST-8341412, AST-8612013, AST-9017082, AST-9319283, AST-9619429 to BTD, in part by support from the Office of Naval Research Young Investigator Program to PJF, and in part by DuPont Corporate Educational Assistance to PJF.

## Appendix A Understanding and Modifying ddscat.par

In order to use DDSCAT to perform the specific calculations of interest to you, it will be necessary to modify the ddscat.par file. Here we list the sample ddscat.par file, followed by a discussion of how to modify this file as needed. Note that all numerical input data in DDSCAT is read with free-format READ(IDEV,\*)... statements. Therefore you do not need to worry about the precise format in which integer or floating point numbers are entered on a line. The crucial thing is that lines in ddscat.par containing numerical data have the correct number of data entries, with any informational comments appearing after the numerical data on a given line.

```
’ =================== Parameter file ===================’
’**** PRELIMINARIES ****’
’NOTORQ’= CMTORQ*6 (DOTORQ, NOTORQ) -- either do or skip torque calculations
’PBCGST’= CMDSOL*6 (PBCGST, PETRKP) -- select solution method
’GPFAFT’= CMETHD*6 (GPFAFT, BRENNR, TMPRTN, CONVEX)
’LATTDR’= CALPHA*6 (LATTDR, LDRISO, GOBR88, DRAI88)
’NOTBIN’= CBINFLAG (ALLBIN, ORIBIN, NOTBIN)
’NOTCDF’= CNETFLAG (ALLCDF, ORICDF, NOTCDF)
’RCTNGL’= CSHAPE*6 (FRMFIL,ELLIPS,CYLNDR,RCTNGL,HEXGON,TETRAH,UNICYL,UNIELL)
8 6 4 = shape parameters PAR1, PAR2, PAR3
1         = NCOMP = number of dielectric materials
’TABLES’= CDIEL*6 (TABLES,H2OICE,H2OLIQ; if TABLES, then filenames follow...)
’diel.tab’
’**** CONJUGATE GRADIENT DEFINITIONS ****’
0       = INIT (TO BEGIN WITH |X0> = 0)
1.00e-5 = ERR = MAX ALLOWED (NORM OF |G>=AC|E>-ACA|X>)/(NORM OF AC|E>)
’**** ANGLES FOR CALCULATION OF CSCA, G’
33      = ICTHM (number of theta values for evaluation of Csca and g)
12      = IPHM (number of phi values for evaluation of Csca and g)
’**** Wavelengths (micron) ****’
6.283185 6.283185 1 ’INV’ = wavelengths (first,last,how many,how=LIN,INV,LOG)
’**** Effective Radii (micron) **** ’
1. 1. 1  ’LIN’ = eff. radii (first, last, how many, how=LIN,INV,LOG)
’**** Define Incident Polarizations ****’
(0,0) (1.,0.) (0.,0.) = Polarization state e01 (k along x axis)
2 = IORTH  (=1 to do only pol. state e01; =2 to also do orth. pol. state)
1 = IWRKSC (=0 to suppress, =1 to write ".sca" file for each target orient.
’**** Prescribe Target Rotations ****’
 0.   0.  1  = BETAMI, BETAMX, NBETA (beta=rotation around a1)
 0.  90.  3  = THETMI, THETMX, NTHETA (theta=angle between a1 and k)
 0.   0.  1  = PHIMIN, PHIMAX, NPHI (phi=rotation angle of a1 around k)
’**** Specify Scattered Directions ****’
0.  0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A
90. 0. 180. 30 = phi, ... for plane B

```

| Lines | Comments |
| --- | --- |
| 1-2 | comment lines – no need to change. |
| 3 | NOTORQ if torque calculation is not required; |
|  | DOTORQ if torque calculation is required. |
| 4 | PBCGST is recommended; other option is PETRKP (see §14). |
| 5 | GPFAFT recommended; other options are BRENNR, TMPRTN, CONVEX13). |
| 6 | LATTDR (LATTice Dispersion Relation) is recommended (§11). |
| 7 | ALLBIN for unformatted binary dump (§10.2; |
|  | ORIBIN for unformatted binary dump of orientational averages only; |
|  | NOTBIN for no unformatted binary output. |
| 8 | ALLCDF for output in netCDF format (must have netCDF option enabled; cf. §10.3); |
|  | ORICDF for orientational averages in netCDF format (must have netCDF option enabled). |
|  | NOTCDF for no output in netCDF format; |
| 9 | specify choice of target shape (see §19 for description of options RCTNGL, ELLIPS, TETRAH, …) |
| 10 | shape parameters SHPAR1, SHPAR2, SHPAR3, … (see §19). |
| 11 | number of different dielectric constant tables (§12). |
| 12 | name(s) of dielectric constant table(s) (one per line). |
| 13 | comment line – no need to change. |
| 14 | 0 is recommended value of parameter INIT. |
| 15 | ERR = error tolerance $`h`$: maximum allowed value of $`|A^{}EA^{}AP|/|A^{}E|`$ \[see eq.(18)\]. |
| 16 | comment line – no need to change. |
| 17 | ICTHM – number of $`\theta _s`$ values for angular averages (§22). |
| 18 | IPHM – number of $`\varphi _s`$ values for angular averages. |
| 19 | comment line – no need to change. |
| 20 | $`\lambda `$ – first, last, how many, how chosen. |
| 21 | comment line – no need to change. |
| 22 | $`a_{\mathrm{eff}}`$ – first, last, how many, how chosen. |
| 23 | comment line – no need to change. |
| 24 | specify x,y,z components of (complex) incident polarization $`\widehat{𝐞}_{01}`$21) |
| 25 | IORTH = 1 to do one polarization state only; |
|  | 2 to do second (orthogonal) incident polarization as well. |
| 26 | IWRKSC = 0 to suppress writing of “.sca” files; |
|  | 2 to enable writing of “.sca” files. |
| 27 | comment line – no need to change. |
| 28 | $`\beta `$ (see §17) – first, last, how many . |
| 29 | $`\mathrm{\Theta }`$ – first, last, how many. |
| 30 | $`\mathrm{\Phi }`$ – first, last, how many. |
| 31 | comment line – no need to change. |
| 32 | $`\varphi _s`$ for first scattering plane, $`\theta _{s,min}`$, $`\theta _{s,max}`$, how many $`\theta _s`$ values; |
| 33,… | $`\varphi _s`$ for 2nd,… scattering plane, … |

## Appendix B wxxryyori.avg Files

The file w00r00ori.avg contains the results for the first wavelength (w00) and first target radius (r00) averaged over orientations (ori.avg). The w00r00ori.avg file generated by the sample calculation should look like the following:

```
 DDSCAT --- DDSCAT.5a.10 [00.11.03]
 TARGET --- Rectangular prism; NX,NY,NZ=   8   6   4
 GPFAFT --- method of solution
 LATTDR --- prescription for polarizabilies
 RCTNGL --- shape
 NAT0 =     192 = number of dipoles
  AEFF=   1.00000 = effective radius (physical units)
  WAVE=   6.28319 = wavelength (physical units)
   K*A=   1.00000 = 2*pi*aeff/lambda
 n= ( 1.3300 ,  0.0100),  epsilon= (  1.7688 ,  0.0266)  for material 1
   TOL= 1.000E-05 = error tolerance for CCG method
 ICTHM= 33 = theta values used in comp. of Qsca,g
 IPHIM= 12 = phi values used in comp. of Qsca,g
 ( 1.00000  0.00000  0.00000) = target axis A1 in Target Frame
 ( 0.00000  1.00000  0.00000) = target axis A2 in Target Frame
 ( 0.27942  0.00000  0.00000) = k vector (latt. units) in Lab Frame
 ( 0.00000, 0.00000)( 1.00000, 0.00000)( 0.00000, 0.00000)=inc.pol.vec. 1 in LF
 ( 0.00000, 0.00000)( 0.00000, 0.00000)( 1.00000, 0.00000)=inc.pol.vec. 2 in LF
   0.000   0.000 = beta_min, beta_max ;  NBETA = 1
   0.000  90.000 = theta_min, theta_max; NTHETA= 3
   0.000   0.000 = phi_min, phi_max   ;   NPHI = 1
 Results averaged over   3 target orientations
                   and   2 incident polarizations
          Qext      Qabs       Qsca     g=<cos>    Qbk       Qpha
 JO=1:  1.330E-01 3.269E-02 1.003E-01  2.345E-01 5.552E-03 4.836E-01
 JO=2:  9.063E-02 2.397E-02 6.666E-02  2.673E-01 3.430E-03 4.159E-01
 mean:  1.118E-01 2.833E-02 8.347E-02  2.476E-01 4.491E-03 4.497E-01
 Qpol=  4.234E-02                                   dQpha= 6.768E-02
        Qsca*g(1)  Qsca*g(2)  Qsca*g(3)
 JO=1:  2.352E-02  1.126E-03 -7.190E-10
 JO=2:  1.781E-02  2.768E-03  3.244E-11
 mean:  2.066E-02  1.947E-03 -3.433E-10
            ** Mueller matrix elements for selected scattering directions **
theta   phi   S_11       S_21       S_31       S_41       S_12       S_13      Pol.
  0.0   0.0 5.167E-02  7.916E-03  9.582E-11  1.267E-11  7.916E-03  9.239E-11  0.15320
 30.0   0.0 4.046E-02  4.518E-04  2.540E-11  4.619E-11  4.518E-04  2.339E-11  0.01117
 60.0   0.0 2.169E-02 -1.023E-02  8.752E-11  5.628E-11 -1.023E-02  6.967E-11  0.47184
 90.0   0.0 1.193E-02 -1.192E-02  2.068E-11  6.778E-11 -1.192E-02 -1.481E-11  0.99888
120.0   0.0 1.170E-02 -5.702E-03 -4.564E-11  3.005E-11 -5.702E-03  5.392E-11  0.48743
150.0   0.0 1.403E-02  9.758E-04 -2.008E-11 -1.155E-11  9.758E-04  1.665E-11  0.06956
180.0   0.0 1.411E-02  3.333E-03 -1.390E-11  1.454E-11  3.333E-03  1.643E-11  0.23621
  0.0  90.0 5.167E-02 -7.916E-03 -7.574E-10  4.831E-11 -7.916E-03 -7.428E-10  0.15320
 30.0  90.0 4.487E-02 -1.280E-02 -5.036E-04  8.979E-05 -1.281E-02 -3.843E-04  0.28550
 60.0  90.0 3.040E-02 -2.062E-02 -7.990E-04  1.643E-04 -2.063E-02 -3.773E-04  0.67857
 90.0  90.0 1.991E-02 -1.989E-02 -7.361E-04  1.827E-04 -1.991E-02 -3.807E-05  0.99926
120.0  90.0 1.626E-02 -1.175E-02 -4.415E-04  1.359E-04 -1.176E-02  1.634E-04  0.72285
150.0  90.0 1.478E-02 -5.293E-03 -1.716E-04  6.554E-05 -5.295E-03  1.131E-04  0.35833
180.0  90.0 1.411E-02 -3.333E-03  3.067E-10 -2.982E-11 -3.333E-03 -3.009E-10  0.23621

```

## Appendix C wxxryykzzz.sca Files

The w00r00k000.sca file contains the results for the first wavelength (w00), first target radius (r00), and first orientation (k000). The w00r00k000.sca file created by the sample calculation should look like the following:

```
 DDSCAT --- DDSCAT.5a.10 [00.11.03]
 TARGET --- Rectangular prism; NX,NY,NZ=   8   6   4
 GPFAFT --- method of solution
 LATTDR --- prescription for polarizabilies
 RCTNGL --- shape
 NAT0 =     192 = number of dipoles
  AEFF=   1.00000 = effective radius (physical units)
  WAVE=   6.28319 = wavelength (physical units)
   K*A=   1.00000 = 2*pi*aeff/lambda
 n= ( 1.3300 ,  0.0100),  epsilon= (  1.7688 ,  0.0266)  for material 1
   TOL= 1.000E-05 = error tolerance for CCG method
 ICTHM= 33 = theta values used in comp. of Qsca,g
 IPHIM= 12 = phi values used in comp. of Qsca,g
 ( 1.00000  0.00000  0.00000) = target axis A1 in Target Frame
 ( 0.00000  1.00000  0.00000) = target axis A2 in Target Frame
 ( 0.27942  0.00000  0.00000) = k vector (latt. units) in TF
 ( 0.00000, 0.00000)( 1.00000, 0.00000)( 0.00000, 0.00000)=inc.pol.vec. 1 in TF
 ( 0.00000, 0.00000)( 0.00000, 0.00000)( 1.00000, 0.00000)=inc.pol.vec. 2 in TF
 BETA =  0.000 = rotation of target around A1
 THETA=  0.000 = angle between A1 and k
  PHI =  0.000 = rotation of A1 around k
          Qext      Qabs       Qsca     g=<cos>    Qbk       Qpha
 JO=1:  1.110E-01 3.028E-02 8.077E-02  3.504E-01 1.858E-03 4.668E-01
 JO=2:  8.651E-02 2.441E-02 6.209E-02  3.650E-01 1.362E-03 4.197E-01
 mean:  9.878E-02 2.735E-02 7.143E-02  3.567E-01 1.610E-03 4.432E-01
 Qpol=  2.454E-02                                   dQpha= 4.711E-02
        Qsca*g(1)  Qsca*g(2)  Qsca*g(3)
 JO=1:  2.830E-02  9.829E-10 -2.070E-09
 JO=2:  2.266E-02  1.827E-09  9.086E-11
 mean:  2.548E-02  1.405E-09 -9.896E-10
            ** Mueller matrix elements for selected scattering directions **
theta   phi   S_11       S_21       S_31       S_41       S_12       S_13      Pol.
  0.0   0.0 4.987E-02  5.372E-03  8.485E-11 -1.807E-11  5.372E-03  8.482E-11  0.10772
 30.0   0.0 4.040E-02 -9.275E-04 -4.151E-11  1.973E-11 -9.275E-04 -4.057E-11  0.02296
 60.0   0.0 2.191E-02 -1.060E-02  1.292E-10  1.107E-11 -1.060E-02  1.205E-10  0.48383
 90.0   0.0 1.058E-02 -1.056E-02 -2.798E-11  7.613E-11 -1.056E-02 -1.444E-10  0.99794
120.0   0.0 7.506E-03 -3.995E-03  3.887E-11 -3.829E-11 -3.995E-03 -4.776E-12  0.53222
150.0   0.0 5.921E-03 -9.862E-06  1.709E-11 -8.195E-13 -9.862E-06 -1.717E-11  0.00167
180.0   0.0 5.058E-03  7.791E-04  6.177E-12 -8.607E-12  7.791E-04 -4.347E-12  0.15403
  0.0  90.0 4.987E-02 -5.372E-03 -4.752E-10  4.820E-11 -5.372E-03 -5.005E-10  0.10772
 30.0  90.0 4.285E-02 -1.032E-02 -5.283E-10 -1.079E-11 -1.032E-02 -5.098E-10  0.24082
 60.0  90.0 2.736E-02 -1.773E-02 -2.169E-10  7.855E-11 -1.773E-02 -2.703E-10  0.64798
 90.0  90.0 1.540E-02 -1.539E-02 -2.506E-11  2.845E-11 -1.539E-02 -1.588E-10  0.99942
120.0  90.0 9.786E-03 -6.745E-03  1.827E-11  3.636E-12 -6.745E-03 -7.262E-11  0.68918
150.0  90.0 6.389E-03 -1.820E-03  4.853E-11  2.225E-11 -1.820E-03 -5.195E-11  0.28491
180.0  90.0 5.058E-03 -7.791E-04  5.552E-11  2.060E-12 -7.791E-04 -5.362E-11  0.15403

```