====================================================
Release 0.5 of Electric Field Monte Carlo code (EMC)
Dr. Min Xu
see COPYRIGHT
====================================================
1. Pre-requirements
EMC depends on the installation of certain softwares.
(a) a C++ compiler
GNU C++ should be fine. The compiler should support templates.
(b) MIEV package by Dr. Warren J. Wiscombe
(wiscombe@climate.gsfc.nasa.gov).
This is an excellent code for Mie calculations.
(c) mtwist-0.5 package by Prof. Geoff Kuenning at Harvey Mudd College
(geoff@cs.hmc.edu).
This is an Mersenne Twist Pseudorandom Number Generator
Package. The Mersenne Twist method for generating pseudorandom
numbers is an extremely fast, flexible, and desirable approach to
random-number generation. It has superb statistical properties and
a ridiculously long period (2^19937-1).
(d) Netcdf package (optional)
Netcdf API is used in some examples coming with EMC to store
simulation outputs. The benefit of using netcdf files is the
platform-independence of these files.
I have included parts of (a-c) in the distribution only for the purpose
your convenience. Full credit should be given to respective authors.
2. EMC
EMC package is implemented in C++ and contained in three headers
files. We use x as the size parameter of the particle, m the complex
relative refractive index, S1 and S2 are the diagonal elements of the
amplitude scattering matrix (whose off-diagonal elements are zero) in
convention of Hulst.
dmiev.h: a C interface for MIEV package, used by scatterer.h
scatterer.h: a C++ header implementing a class "scatterer". The key
methods include:
// N and nslot can be increased for a higher accuracy
// of the forward and inverse table.
scatterer(double x, dcmplx m, int N=10001, int nslot=10000);
// obtain S1 and S2 at a list of cosines of angles of
// length n
phasef(int n, double* mu, dcmplx* s1, dcmplx* s2);
// obtain S1 and S2 at the cosine of one angle mu
phasef(double mu, dcmplx* s1, dcmplx* s2);
// a quicker version using the lookup table pre-computed
phasef_lu(double mu, dcmplx* s1, dcmplx* s2);
// yield the scattering angle given the probability
// within (0,1) using the inverse table pre-computed
draw_mu(double p);
pol_montecarlo.h: the core header file of EMC implementing a
class "photonPacket". The key methods include:
// initialization
photonPacket(const scatterer* sct, char* fname=NULL,
unsigned long seed=0);
// launch the photons. The incident electric field is
// given by (E1, E2) where light is propagating in the
// direction (u, v, w), E1 is in the direction (l, m,
// n), and E2 is in the direction specified by the cross
// product of the above two directions. Light is
// incident upon the position (x,y,z) at time t. The
// incident light intensity must be unity (|E1|^2 +
// |E2|^2=1) and the directions must be unit vectors.
void launch(dcmplx E1=1, dcmplx E2=0, double l=1,
double m=0, double n=0, double u=0, double v=0,
double w=1, double x=0, double y=0, double z=0,
double t=0);
// move to next scattering or absorption event
void move();
// absorption event
void absorb();
// scattering event
void scatter();
// next event estimator of (Ed1, Ed2) with the
// direction vector Q at time td crossing the boundary
// z=zd in the direction (ud, vd, wd)
void pointestimator(double *td, double *deposit,
double Q[3][3], dcmplx *Ed1, dcmplx *Ed2, double zd,
double ud=0, double vd=0, double wd=1);
The photon status is stored in the public data members
of the class, they are:
xold, yold, zold, told:
the position and time of the previous scattering
event
x, y, z, t, weight:
the current position and time of photon and its weight
E1, E2:
electric field components
P[3][3]:
(l, m, n) the unit vector for E1 direction
(p, q, r) the unit vector for E2 direction
(u, v, w) the propagation direction
nsct:
number of scattering events