Free F77 computer codes together with an overview of soon available codes * Core-shell (nano)particles * Nanoshells * Multilayered sphere * Spheres, spheroids, rods, pillar, cylinders, Chebyshev particles * T-matrix code * Photonic crystals * Optical trapping * Lattice sums * Ewald summation * Free software

Free F77 computer codes together with an overview of soon available codes

The codes may be used without limitations in any not-for-profit scientific research. The only request is that in any publication using the codes, the source of the code be acknowledged and relevant references be made.


Download new material data file Audat.dat.txt to be used with NFIN=77. (For security reasons, *.dat files cannot be posted; therefore the additional *.txt extension. You have to rename the data file back to Audat.dat after the download.)
The earlier data file followed a typing error in Palik I, p. 294, where the wavelengths of 1.384 and 1.378 micrometers should be interchanged.

Single-scatterer programs

Scattering from a multilayered sphere

The routine returns the total extinction, dipole extinction, quadrupole extinction, (elastic) scattering, and absorption cross sections, albedo, and, optionally, the T-matrix elements for the scattering off a single multilayered sphere. The number LCS of concentric sphere layers can be arbitrary. The program has built-in checks for convergence of the output values. By changing the parameter NMAT value, you can switch between plurality of different materials. Currently, the following options are available: dispersionless dielectric, Drude metal, Ag, Au, ZnS, Cu, Al, Pt, and Si. The option of a perfect conductor (the Dirichlet boundary condition) is activated by the .TRUE. value of the logical parameter YNPERFCON. The algorithm has been described in details in Sec. 3 of my article The routine has been an integral part of my multiple-scattering programs since the very beginning. See, for instance, and more recently a nanophotonic application in The program has also been used by and in where you will find a comparison between theory and experiment.

You can download the source code "sphere.f" here. The instructions how to run and compile it, together with the definition of its parameters, are written within the comment lines at the beginning of the source file. You will also need to download material data file Audat.dat. (For security reasons, *.dat files cannot be posted; therefore the additional *.txt extension. You have to rename the data file back to Audat.dat after the download.)

You can also download the code limited Windows executable sphere, which calculates the scattering off a single SiO2@Au coated sphere in air. Refractive index of SiO2 in the sphere core is taken to be 1.45. Refractive indices of gold are taken to be gold bulk data without any adjustable parameter as read from the data file Audat.dat, which is also to be downloaded.

Optical trapping of a multilayered sphere

The program opttrap.f calculates the axial trapping efficiency Q for a multilayered sphere. The axial trapping efficiency Q, or the dimensionless normalized axial force, when multiplied by the focused laser beam power P, determines the resulting axial force F on a sphere immersed in a homogeneous medium with refractive index n according to the formula F=(n/c) PQ. The program was made according to the theory of P. A. Maia Neto and H. M. Nussenzveig, Theory of optical tweezers, Europhys. Lett. 50, 702-708 (2000) by extending it to the case of (multi)coated spheres. For the results obtained by this program, see: In the first of the two articles, the possibility of optical trapping of coated (Ag core and silica shell) nanospheres was demonstrated. The possibility of optical trapping of metallic particles came as a surpsise, since, at that time, (i) metallic particles were known to be usually repelled from the beam focus by the scattering force and (ii) one usually contemplated only optical trapping of dielectric particles in the most frequently used size regime between 0.2 micrometers and 1.0 micrometers and with refractive indices 1.43 (silica) and 1.57 (polystyrene). For completeness, a complementary case of a three-dimensional optical trapping of partially silvered silica microparticles was then investigated two years later by P. Jordan et al, Opt. Lett. 29, 2488-2490 (2004).

Radiation power of a dipole radiating inside or outside a multilayered sphere

Dipole can be located either outside or embedded in any layer of the multilayered sphere. Provides both radiative and non-radiative decay rates. Avoids the cumbersome algorithm of H. Chew, P. J. McNulty, and M. Kerker, Raman and fluorescent scattering by molecules embedded in concentric spheres, J. Opt. Soc. Am. 66, 440-444 (1976) and provides a new transfer matrix alternative. Important for the description of inelastic light-scattering (fluorescence or Raman) spectroscopy for characterizing single micrometer sized particles (chemical speciation) of both inorganic and organic compounds, aerosols and particulates, in LIDAR applications for remote sensing of both molecular and particulate constituents of atmosphere, for identification of biological particles in fluorescence-activated flow of cytomeres, to monitor specific cell functions, or in the cell identification and sorting systems, in the investigation thermal radiation from spherical microparticles, engineering of the radiative decay for biophysical and biomedical applications, imaging of buried saturated fluorescent molecules and imaging of surfaces in near-field optical microscopy, in the study of the effects of light absorption and amplification on the stimulated transition rates of the electric-dipole emission of atoms or molecules embedded in micro- or nano-structured spheres, stimulated Raman scattering, the interplay between lasing and stimulated Raman scattering, etc.
Theoretical details can be found in my article Some additional results have been published in follow ups

Here you can download F77 codes CHEW and CHEWFS. The first one, CHEW, calculates the radiative and nonradiative (due to Ohmic losses) decay rates (the latter are sometime also called energy transfer rates), whereas CHEWFS calculates induced frequency shifts and total decay rates directly from Green´s function at coinciding spatial arguments. The sphere options are identical to that in scattering from a multilayered sphere. Thus, you can run it for both homogeneous and (multi)coated spheres.
You can also download a corresponding limited Windows executables. Chew.exe (download also material data file Audat.dat, which contains gold bulk data without any adjustable parameter) calculates the electric dipole radiated power and the dipole power loss due to Ohmic losses for a coated SiO2@X@SiO2 sphere in water. Refractive index of SiO2 is taken to be 1.45, that of water 1.33, and that of X you can supply yourself. For a Windows executable generated from CHEWFS, see an on-line Appendix A of my Chem. Phys. 317(1), 1-15 (2005).
In case you are using a Linux or Unix machine, here are the respective makefiles mkchew and mkchewfs. Some of the referred routines are attached at the end of the source files.

Effective medium parameters for a composite of (coated) spheres in the long wavelength limit

Even though all the constituents media are nonmagnetic, the effective medium can nevertheless exhibit magnetic properties. In fact, a recent article by C. L. Holloway, E. F. Kuester, J. Baker-Jarvis, and P. Kabos, A double negative (DNG) composite medium composed of magnetodielectric spherical particles embedded in a matrix, IEEE Trans. Antennas Propag. 51, 2596-2603 (2003) suggested that a negative-refraction material can, in principle, be made as a composite of (not necessarily magnetic) spheres. Together with Vassilis Yannopapas [see our manuscript J. Phys.: Condens. Matter. 17, 3717-3734 (2005)] we have recently expanded on the ideas of Holloway et al. and showed that even a composite of inherently non-magnetic homogeneous spheres can provide a negative refractive index metamaterial. The absence of inherently magnetic materials in our porposal is a key difference which makes it possible to achieve a negative refractive index band even in the deep infrared region. Our results are explained in the context of the extended Maxwell-Garnett theory (see my F77 code EFFE2P) and reproduced by the ab initio calculations based on multiple scattering theory.
Here you can also download a limited Windows executable effe (download also material data file Audat.dat, which contains gold bulk data without any adjustable parameter), which calculates an effective medium properties for a composite of coated Au@SiO2 spheres in air with volume fraction of 34%. The code and executable are to be used for wavelengths typically at least one order of magnitude longer than the sphere radius.

Scattering from an axially symmetric particle

If you use Mishchenko´s code at you will often find that its non-NAG alternative for a nonabsorbing dielectric particle and both NAG and non-NAG alternatives in a strongly absorbing case fail to converge, complain about singular matrices, and do not provide any output. An improvement of the Mishchenko code is provided which expands its range especially in the direction of strongly absorbing (e.g., metallic) particles (see A. Moroz, Improvement of Mishchenko's T-matrix code for absorbing particles, Appl. Opt. 44, 3604-3609 (2005) (accompanying F77 code is available here). Rather simple amendment described in the article allows one to extend Mishchenko´s code, which has been used to generate beautiful benchmarking results for dielectric particles of various axisymmetric shapes such as [see free online book by M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University Press, Cambridge, 2002)]: to the very interesting field of metallic nanoparticles. The latter are a hot subject with many potential applications for (bio)sensing, nanolitography, multicolor labelling, nanorulers, SERS, etc.

Additional improvements not available on-line yet:
Further particle shapes (a cylindrical cone and a sphere cut by a plane) have been included. By changing the value of parameter NMAT you can switch between plurality of different materials: dispersionless dielectric, Drude metal, Ag, Au, ZnS, Cu, Al, Pt, and Si. The option of a perfect conductor (the Dirichlet boundary condition) is activated by the .TRUE. value of the logical parameter YNPERFCON.
A comparison of theory and experiment for the case of metal coated spheroids have been dealt with in and a comparison between various long-wavelength approximations against the exact T-matrix results has been provided in

Electric intensity around an axially symmetric particle

Indispensible tool for the SERS, near-field imaging, near-Field Raman Microscopy, near-field scanning optical microscopy (NSOM), design of near-field optical probes, etc. The particle options are identical to that in scattering from an axially symmetric particle.

Free-space periodic Greens functions calculations

The set of programs calculates a free-space Green´s function of the Helmholtz equation in a given space dimension d which satisfies a Bloch periodic boundary condition in a subspace of the d-dimensional space, including the d-dimensional embedding space itself. The programs employ exponentially convergent representations of lattice sums as derived for different cases by Ewald, Ham and Segall, Kambe, Ozorio de Almeida, and, in the case of d=2 and Bloch periodicity in one dimension, by

and, in the case of d=3 and Bloch periodicity in one dimension, by

The total lattice sum D_L is conventionally written as a sum of three partial contributions

D_L = D_L^1 + D_L^2 + D_L^3

If you look at my J. Phys. A: Math. Gen. 39, 11247-11282 (2006), the D_L^1 term for different subperiodic cases is given by Eqs. 83,85,86. The only special function is the incomplete gamma function which is easily calculated numerically on using recurrences given by Eqs. 88-89. There is nothing mysterious or difficult in calculating incomplete gamma function compared to e.g. calculation of the Bessel functions.

The D_L^2 term for different subperiodic cases is given by Eqs. 102-103. There is no special function involved. An integral therein is calculated easily numerically again on using recurrence given by Eq. 106.

The D_L^3 term for different subperiodic cases is given by Eqs. 118-119. There is no special function involved and the sum is calculated numerically in a straightforward way.

A historical overview, relevant definitions and references can be found in the latter article and also here.
The set of programs forms an integral part of my multiple-scattering programs for scattering off a periodic arrangement of scatterers. The numerical code OLA for a 1D periodicity in 2D is available here.

For a 2D periodicity in 3D, see an Appendix of the monograph by

which provides a full listing of the Pendry's LEED code. The routines of Pendry's LEED code can be used even by a pragmatic engineer as black box routines. The Pendry's LEED code comprises among others also, up to minor amendmends, the program lines that I have arranged in the code dlsumf2in3.f, which performs the Kambe's lattice sums.

My code gff2in3.f can be then used as a main code. It generates the lattice points, calls the routines that calculate spherical harmonics and spherical Bessel functions, and helps you to obtain the free-space Green's function given the total lattice sum D_L generated by dlsumf2in3.f. The code has been employed in my publication

Multiple-scattering programs for scattering off a periodic arrangement of scatterers

An optimization program for a finite periodic stacks of parallel layers

Implemenets the well-known Abeles transfer matrix method (reproduced in Sec. 1.6 "Wave propagation in a stratified medium" and Chapter 14 of Born and Wolf). Can be used to optimize a finite stack parameters with respect to reflectivity, transmittivity, and absorption. I have used it to optimize the properties of metallo-dielectric stacks in my publication

One-dimensional local density of states

The source code ldosl1d.f calculates the local density of states for a toy model of a one-dimensional photonic crystal in one space dimension, which has been discussed in my article

Two-dimensional layered KKR method

Solves the problem of scattering in two dimensions off a finite stack of discs arranged periodically along parallel lines. Used to model the properties of a two-dimensional periodic arrangement of infinite-length circular cylinders in the plane normal to the cylinder axis. Calculates reflection, transmission, and absorption of both classical (water, acoustic, and electromagnetic) and quantum (electron) waves. It is unique in a sense that it allows for an optimization of the stack parameters with respect to a maximal band gap, minimal absorbtion, etc. On a PC, up to cca 15.000 different configurations can be scanned overnight with program selecting the top configurations, whether you want those with minimal absorption, the largest bandgap, or otherwise. For a comparison of theory and experiment, see some preliminary results in Following the link rta1in2k, you can download a limited Windows executable (download also material data file Audat.dat) which calculates reflection, transmission, and absorption (RTA) off a finite 2D photonic crystal stack of silver rods. The executable has been created using Compaq Visual F77 ==> an additional gain in speed could have been obtained by compiling the source on a UNIX machine with a suitable optimization.

Two-dimensional bulk KKR method

Solves the problem of scattering in an infinite periodic arrangement of discs filling in entire two-dimensional space. Calculates the two-dimensional band structure of both classical (water, acoustic, and electromagnetic) and quantum (electron) waves. It has been developed in collaboration with Han van der Lem. For the results in the electromagnetic case, see

Three-dimensional layered KKR method

Solves the problem of scattering in three dimensions off a finite stack of axially symmetric scatterers arranged periodically along parallel (two-dimensional) planes. Similar to the two-dimensional layered KKR method, but in one dimension higher. Improves its predecessor, the computer code of V. Yannopapas, N. Stefanou, and A. Modinos, Heterostructures of photonic crystals: frequency bands and transmission coefficients, Comput. Phys. Commun. 113, 49 (1998) [The scalar predecessor of the latter program has been developed in the book by J. B. Pendry, Low Energy Electron Diffraction (Academic Press, London, 1974)], in the following aspects: You will find many theoretical results generated by this program in For a comparison of theory and experiment, see Results for nonspherical particles can be found in Sec. 5 of A limited Windows executable which incorporates the lattice sums within a photonic LKKR code for the calculation of reflection, transmission and absorption of an electromagnetic plane wave incident on a square array of finite length cylinders arranged on a homogeneous slab of finite thickness is available following the link . In order for the executable to run, it is required that you download material data file Audat.dat.txt (for security reasons, *.dat files cannot be posted; therefore the additional *.txt extension), which you have to save as Audat.dat in the same directory as the executable.

Three-dimensional bulk KKR method

Solves the problem of scattering in three dimensions in an infinite three-dimensional lattice of axially symmetric scatterers. Three-dimensional variant of the two-dimensional bulk KKR method with identical material and multilayered sphere options. Regarding the band-structure calculations, the bulk KKR method is much more robust and stable than the three-dimensional layered KKR method. The reason is that, unlike in the latter method, no matrix inversion is involved in calculating the relevant secular matrix. The program was heavily employed to generate results in several my articles, which representative selection is as follows: The code will be released upon publication of

PhD theses where my codes have been used

Description of some common parameters in alphabetical order

Usually, only a single sphere layer comprises dispersive material. ILCS defines the number of the layer (when numbered from the core outwards). Hence ILCS=1 if the core material is dispersive, ILCS=LCS if the outermost shell material is dispersive.

LCS defines the total number of concentric sphere layers. Sphere core is also considered to be a layer. Therefore, LCS=1 for a homogeneous spheres, LCS=2 for a spheres with a single shell, etc. The number LCS can be arbitrary.

LMAX is a convergence parameter. It defines the cut-off on the value of the orbital angular momentum of spherical waves taken into account in the calculation.

The parameter defines the dimension of a corresponding material data file.

By changing the parameter NMAT value, you can switch between plurality of different materials. Currently, the following options are available: dispersionless dielectric, Drude metal, Ag, Au, ZnS, Cu, Al, and Si. The option of a perfect conductor (the Dirichlet boundary condition) is activated by the .TRUE. value of the logical parameter YNPERFCON. The dielectric constant of metals for wavelengths up to cca 2 micrometers is determined by a linear interpolation between experimental data according to Palik's Handbook of Optical Constants of Solids. For longer wavelengths a Drude fit is used according to parameters by M. A. Ordal et al, Optical properties of the metals Al, Co, Cu, Au, Fe, Pb, Ni, Pd, Pt, Ag, Ti, V, and W in the infrared and far infrared, Appl. Opt. 22, 1099-1119 (1983). Regarding other materials, see Luxpop Index of refraction values.

Defines the output unit.

The parameter specifies the shape of particle from a predefined list.

By changing the value of NSTACK you can switch between different stacking sequences. Currently 12 in three space dimensions and 2 in two space dimensions.

The option of a perfect conductor (the Dirichlet boundary condition) is activated by the .TRUE. value of the logical parameter YNPERFCON. In all other cases, the value of YNPERFCON should be .FALSE.


Most of my programs use the AMOS package of complex Bessel functions routines which you can, together with their dependencies, freely download at Netlib. Before using the AMOS package, machine constants in the real function D1MACH(I) and in the integer function I1MACH(I) have to be adjusted to your computer. When using Microsoft Windows, you can activate the same machine constants as listed therein under the option of Silicon Graphics IRIS.

In writing my programs, I have occasionally made use of routines from Numerical Recipes. For copyright reasons, these routines are not distributed here. (I have not received any reply yet from Numerical Recipes Software on my request for a permission to distribute the routines.)


My programs are in double precision. Therefore, before you make use of the relevant Numerical Recipes routines, which are by default in single precision, the latter have to be transformed into double precision, i.e., you have to change change declaration real to real*8 and complex to complex*16.

In the case of the free-space one-dimensional periodic Green function in two space dimensions [see A. Moroz, Exponentially convergent lattice sums, Opt. Lett. 26, 1119-1121 (2001) [pdf]], some additional instructions can be found here.



I would highly appreciate any comments and/or bug reports and informing me of any problems encountered with any of the above codes. If you choose to use them, please send email to

wavescattering  at

registering as a user; registered users will be notified when updates to the code(s) are made. If you publish results obtained using the above codes, please acknowledge the source of the code.

[ Main page | My personal webpage | Full list of my scientific publications | List of my publications in Scitation | Top-10 of my most frequently cited articles ]
[ Selected Links on Photonics, Photonic Crystals, Numerical Codes, Free Software ]

© Alexander Moroz, September 23, 2003 (last updated on September 25, 2010)

eXTReMe Tracker