Free F77 computer codes together with
an overview of soon available codes * Coreshell (nano)particles *
Nanoshells *
Multilayered sphere *
Spheres, spheroids, rods, pillar, cylinders, Chebyshev particles * Tmatrix 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 notforprofit 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.
Important!!!
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.
Singlescatterer 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
Tmatrix elements for the scattering off a single multilayered sphere.
The number LCS of concentric sphere layers can be arbitrary.
The program has builtin 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 multiplescattering programs since the very beginning.
See, for instance,
and more recently a nanophotonic application in
The program has also been used by

Christina Graf and Alfons van Blaaderen,
Metallodielectric Colloidal CoreShell Particles for Photonic Applications,
Langmuir 18, 524534 (2002),
and in
 J. J. Penninkhof, A. Moroz, A. Polman, and A. van Blaaderen,
Optical properties of spherical and oblate spheroidal gold shell colloids,
J. Phys. Chem. C 112(11), 41464150 (2008)
[story behind this article],
 A. Moroz, Localized resonances of composite particles,
J. Phys. Chem. C 113(52), 2160421610 (2009),
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, 702708 (2000) by extending it to the case of
(multi)coated spheres. For the results obtained by this program, see:

A. Moroz, Metallodielectric diamond and zincblende photonic
structures,
Phys. Rev. B 66, 115109 (2002)
[condmat/0209188]
[pdf].
(See discussion subsection "Fabrication" and Fig. 23 therein for optical trapping
of metal coredielectric shell particles.)
 A. van der Horst, P. D. J. van Oostrum, A. Moroz, A. van Blaaderen, and M. Dogterom,
Highrefractive index particles
trapped in dynamic arrays of dualbeam optical tweezers,
Appl. Opt. 47(17),
31963202 (2008).
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 threedimensional optical trapping of partially silvered
silica microparticles was then investigated two years later by
P. Jordan et al,
Opt. Lett. 29, 24882490 (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 nonradiative 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, 440444 (1976) and
provides a new transfer matrix alternative.
Important for the description of inelastic lightscattering (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 fluorescenceactivated 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 nearfield
optical microscopy, in the study of the effects of light absorption
and amplification on the stimulated transition rates of the electricdipole
emission of atoms or molecules embedded in micro or nanostructured 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
online Appendix A of my
Chem. Phys. 317(1), 115 (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. BakerJarvis, and P. Kabos,
A double negative (DNG) composite medium composed of magnetodielectric spherical
particles embedded in a matrix,
IEEE Trans. Antennas Propag. 51, 25962603 (2003) suggested that a
negativerefraction 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, 37173734 (2005)]
we have recently expanded on the ideas of
Holloway et al. and showed that even a composite of inherently
nonmagnetic 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
MaxwellGarnett 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 http://www.giss.nasa.gov/~crmim
you will often find that its nonNAG alternative for a nonabsorbing dielectric particle
and both NAG and nonNAG 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 Tmatrix code for absorbing
particles,
Appl. Opt. 44,
36043609 (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)]:
 spheroids/ellipsoids
 Chebyshev particles (you can roughly imagine a Chebyshev particle
as a corrugated sphere with the Chebyshev particle order yielding the number
of corrugations)
 finite cylinders
 cylinders with a hemispheric ends
 truncated spheres
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 online 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 longwavelength approximations against the exact
Tmatrix results has been provided in
 A. Moroz, Depolarization field of spheroidal particles,
J. Opt. Soc. Am. B 26(3), 517527 (2009).
(Accompanying F77 code sphrd.f to determine extinction efficiency
in various longwavelength approximations discussed in the article.
See here for online supplementary material.
)
Electric intensity around an axially symmetric particle
Indispensible tool for the SERS, nearfield imaging, nearField Raman Microscopy,
nearfield scanning optical microscopy (NSOM),
design of nearfield optical probes, etc.
The particle options are identical to that in
scattering from an axially symmetric particle.
Freespace periodic Greens functions calculations
The set of programs calculates a freespace 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 ddimensional space,
including the ddimensional
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, 1124711282 (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. 8889. 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. 102103.
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. 118119.
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 multiplescattering 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 freespace Green's function given the total lattice sum D_L
generated by dlsumf2in3.f.
The code has been employed in my publication
Multiplescattering programs for scattering off
a periodic arrangement of scatterers
An optimization program for a finite periodic stacks of parallel layers
Implemenets the wellknown 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 metallodielectric
stacks in my publication
Onedimensional local density of states
The source code ldosl1d.f calculates the local density of
states for a toy model
of a onedimensional photonic crystal in one space dimension, which
has been discussed in my article
Twodimensional 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 twodimensional
periodic arrangement of infinitelength 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.
Twodimensional bulk KKR method
Solves the problem of scattering in an infinite periodic
arrangement of discs filling in entire twodimensional space.
Calculates the twodimensional 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
Threedimensional layered KKR method
Solves the problem of scattering in three dimensions off a finite stack of
axially symmetric scatterers arranged
periodically along parallel (twodimensional) planes.
Similar to the twodimensional 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:
 is applicable to lattices with more than one scatterer in the
lattice primitive cell (the socalled complex lattice case)
 can be applied not only to the lattices of spheres but also to a lattice
of arbitrary axiallysymmetric particles
(The sphere options are identical
to that in scattering from a multilayered sphere;
a general axiallysymmetric scatterer options are identical to those
in scattering from an axially symmetric particle.)
 in addition to the reflection, transmission, absorption it also outputs the
change in the phase of the specular beam upon reflection
 includes the option of a multilayered sphere. The number LCS
of concentric sphere layers can be arbitrary
 by controlling the parameter value of NSTACK, you can switch between
12 different stacking sequences
 by controling the parameter value of 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.
 checks the stability of the Ewald summation
 has an option to output reflection, transmission, absorption of
individual diffraction orders separately
You will find many theoretical results generated by this program in
For a comparison of theory and experiment, see

K. P. Velikov, A. Moroz, and A. van Blaaderen,
Photonic crystals of coreshell
colloidal particles,
Appl. Phys. Lett. 80, 4951 (2002)
[pdf].
 D. A. Mazurenko, A. Moroz, C. M. Graf, A. van Blaaderen, and
J. I. Dijkhuis, Threedimensional silicagold coreshell photonic crystal:
linear reflection and ultrafast nonlinear optical properties,
Photonic Crystal Materials and Nanostructures, Photonics Europe 2004,
Proceedings of SPIE Vol. 5450 (2004), pp. 569577.
[pdf]
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.
Threedimensional bulk KKR method
Solves the problem of scattering in three dimensions in an infinite
threedimensional lattice of axially symmetric scatterers.
Threedimensional variant of the twodimensional bulk KKR method with identical
material and multilayered sphere options.
Regarding the bandstructure calculations, the bulk KKR method is much
more robust and stable than the threedimensional 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:

A. Moroz and C. Sommers,
Photonic band gaps of threedimensional facecentered cubic lattices,
J.
Phys.: Condens. Matter 11, 9971008 (1999)
[physics/9807057]
[story behind this article],

A. Moroz, Threedimensional complete photonic bandgap structures
in the visible,
Phys. Rev. Lett. 83, 52745277 (1999),
 A. Moroz, Metallodielectric diamond and zincblende photonic
structures,
Phys. Rev. B 66, 115109 (2002)
[condmat/0209188]
[pdf].
The code will be released upon publication of
 E. C. M. Vermolen, J. H. J. Thijssen, A. Moroz, M. Megens, and A. van Blaaderen,
Comparing photonic band structure calculation methods for diamond and pyrochlore crystals,
Optics Express 17(9), 69526961 (2009).
PhD theses where my codes have been used
Description of some common parameters in
alphabetical order
ILCS
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
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
LMAX is a convergence parameter. It defines the cutoff on the value of the
orbital angular momentum of spherical waves taken into
account in the calculation.
NFIN
The parameter defines the dimension of a corresponding material data file.
NMAT
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, 10991119 (1983).
Regarding other materials, see
Luxpop Index of refraction values.
NOUT
Defines the output unit.
NP
The parameter specifies the shape of particle from a predefined list.
NSTACK
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.
YNPERFCON
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.
Requirements
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.)
Important!!!
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 freespace onedimensional periodic Green function
in two space dimensions
[see A. Moroz, Exponentially convergent lattice sums,
Opt. Lett. 26,
11191121 (2001)
[pdf]], some
additional instructions can be found here.
Disclaimer
WHILE THE COMPUTER PROGRAMS HAVE BEEN TESTED FOR A VARIETY OF CASES,
IT IS NOT INCONCEIVABLE THAT THEY MAY CONTAIN UNDETECTED ERRORS. ALSO,
INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
THE AUTHORS DISCLAIM ALL LIABILITY FOR
ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAMS.
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 yahoo.com
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

Top10 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)