740 likes | 930 Views
Monte Carlo, Euler and Quaternions. Jessica C. Ramella-Roman. Stokes vector reference. Yesterday, Meridian plane Three steps “Rotate” to scattering plane Rotational matrix Scatter Scattering matrix “ Rotate” to a new plane Rotational matrix. Today Monte Carlo.
E N D
Monte Carlo, Euler and Quaternions Jessica C. Ramella-Roman
Stokes vector reference • Yesterday, Meridian plane • Three steps • “Rotate” to scattering plane • Rotational matrix • Scatter • Scattering matrix • “Rotate” to a new plane • Rotational matrix
Today Monte Carlo • Reference frame is a triplet of unit vectors • Rotation are about an axis and follow a specific order http://www.grc.nasa.gov/WWW/K-12/airplane
Today Monte Carlo • Reference frame is a triplet of unit vectors • Rotation of the frame is done in two steps.
Euler and Quaternion • The only difference between today programs is the way we handle these rotations • Euler angles rotation Quaternionsalgebra
Euler Monte Carlo • The photon polarization reference frame is tracked at any time via a triplet of unit vectors that are rotated by an azimuth and scattering angle according to a predefined order • Introduced by Bartelet al. (*) * S. Bartel, A. Hielsher, “Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media,” Applied Optics, Vol. 39, No. 10; (2000).
Euler matrices • Euler's rotation theorem: any rotation may be described using three angles. • And three rotation matrices • About z axis • About x’ axis • About z’ axis Drawings, Wolfram, Mathworld
Euler cnt. • We only need TWO vectors and TWO rotations • Two unit vectors vand u • The third unit vector is defined by the cross product of vand u and is calculated only when a photon reaches a boundary.
Euler cnt. • Advantage • Stokes vector is rotated only once for each scattering event instead of twice as in the meridian plane method. • Simple to implement and visualize • Drawback – • Gimbal lock - makes the rotation fail for angles exactly equal to 90˚ (rare event).
Launch • The two unit vector v and u • v and u define the starting Stokes vector reference plane. • The unit vectorurepresents the direction of photon propagation.
Move • The photon is moved a distance Ds, the new coordinates of the photon are
Drop • Reduction of photon weight (W) • According to material albedo • albedo = ms/(ms+ ma) • absorbed = W*(1-albedo)
Scatter • The rejection method establishes the scattering angle q and azimuth anglef. • The Stokes vector must be rotated by an angle f into the scattering plane before scattering can occur.
Referencing to scattering plane • Done as in meridian plane Monte Carlo • Multiply the Stokes vector by matrix R
Scatter cnt • The interaction with a spherical particle is achieved by multiplying the Stokes vector with a scattering matrix M(q)
M(q) elements • s11, s12 , s33 , s34 calculated with Mie theory • Expressed as • S1, S2* – Mie scattering • *http://omlc.ogi.edu/calc/mie_calc.html • C. Bohren and D. R. Huffman, • Absorption and scattering of light by small • particles, (Wiley Science Paperback Series,1998).
Reference frame • After scattering the Stokes vector the reference coordinate system vu must be updated. • The two rotations, for the angles f and q can be obtained using Euler’s rotational matrices.
Rotational matrix ROT • Rodrigues’s rotational matrix ROT • Accomplishes the general case of rotating any vector by an angle y about an axis K • Where K is the rotational axis, • c=cos(y), s=sin(y) and v=1-cos(y).
Vector rotation • First the unit vector v is rotated about the vector u by an angle f • Multiply v by the rotational matrix ROT(u,f); u remains unchanged
Vector rotation • Second u is rotated about the newly generated v by an angle q. • This is done multiplying the unit vector u by the rotational matrix ROT(v,q).
Update of direction cosines • This is done as in standard Monte Carlo • q and f If |uz|≈ 1
Update of direction cosines If |uz|≠ 1
Photon life • The life of a photon ends when the photon passes through a boundary or when its weight W value falls below a threshold. • Roulette is used to terminate the photon packet when W £Wth. • Gives the photon packet one chance of surviving • If the photon packet does not survive the roulette, the photon weight is reduced to zero and the photon is terminated.
Boundaries • Two final rotations of the Stokes vector are necessary to put the photon status of polarization in the detector reference frame. • This will be achieved with matrix multiplication of two rotational matrix by stokes vector
Boundaries –R(e) • w is reconstructed as the cross product of v and u. • An angle e is needed to rotate the Stokes vector into a scattering plane
Boundaries –R(e) • This is the ONLY reason we need the vector u,v, andw
Boundaries –R(e) • This rotation is about the direction of propagation of the photon, i.e. the axis u • Resulting position of v
Boundaries –R(y) • A second rotation of an angle y about the Z axis • Put the photon reference frame in the detector reference frame
Boundaries –R(y) • Transmission • Reflection
Boundaries End • Stokes vector in the detector frame of reference • We need vector u,v, andwto obtain these angles
Quaternion Monte Carlo • Quaternion Monte Carlo simply uses Quaternion Algebra to handle vector rotation. • Advantage • Stokes vector is rotated only once for each scattering event instead of twice as in the meridian plane method. • No issue with Gimbal lock • Optimized for computer simulations
Quaternions • A quaternion is a 4-tuple of real numbers; it is an element of R4. • Quaternion can also be defined as the sum of a scalar part q0 and a vector part Q in R3of the form
Quaternionscnt. • The vector part Q is the rotational axis and the scalar part is the angle of rotation. • Multiplication of a vector t by the quaternion is equivalent to rotating the vector t around the vector Q of an angle q0.
Vector rotation • First the unit vector v is rotated about the vector u by an angle f
Rotation - f • The first rotation is about the vector u by an angle f. • This is done generating the quaternion qf • And then using the quaternion operator • q∗vqon the vector v • q∗ is the complex conjugate of q
Vector rotation • Second u is rotated about the newly generated v by an angle q.
Rotation - q • The second rotation is about the vector v by an angle q. • This is done generating the quaternion qq • And then using the quaternion operator • qq*uqqon the vectoru
Rotations • These steps are repeated for every scattering event. • At the boundaries the last aligning rotations are the same as in Euler Monte Carlo.
Testing • Comparison with Evans Code • In 1991 Evans designed an adding doubling code that can calculate both the radiance and flux of a polarized light beam exiting the atmosphere • plane parallel slab of thickness L = 4/µs, • absorption coefficient µa=0, • unpolarized incident beam • wavelength l = 0.632 nm.
Evans - Reflectance mode Incident I=[1 0 0 0] Reflectance mode, comparison between Evans adding-doubling code and the meridian plane Monte Carlo program. The results do not include the final rotation for a single detector.
Evans – Transmission mode Incident I=[1 0 0 0] Transmission mode, comparison between Evans adding doubling code and the meridian plane Monte Carlo program. The results are not corrected for a single detector
How to run the code • Download the code • http://faculty.cua.edu/ramella/MonteCarlo/index.html
How to run the code • Download the code • http://faculty.cua.edu/ramella/MonteCarlo/index.html
Download the code Programs
Make Command line Compilation/linking step
Run Command line This program will run a full Mueller matrix Monte Carlo launching four vectors [1 1 0 0] H 0 degree polarized [1 -1 0 0] V 90 degrees polarized [1 0 1 0] P 45 degrees polarized [1 0 0 1] R –Right circular
Memory usage • Intel (64-bit) – Macbook pro 2.3Ghz I7 8GB RAM • Real mem 860 KB • Virtual mem 17 MB
Output is a set of Stokes vector images for each launched Stokes vector outHI -> Horizontal incident polarization and I portion of the Stokes vector outHQ -> Horizontal incident polarization and Q portion of the Stokes vector out
Inside the code GNU Licence