Monte Carlo, Euler and Quaternions

1. Monte Carlo, Euler and Quaternions Jessica C. Ramella-Roman

2. Stokes vector reference • Yesterday, Meridian plane • Three steps • “Rotate” to scattering plane • Rotational matrix • Scatter • Scattering matrix • “Rotate” to a new plane • Rotational matrix

3. 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

4. Today Monte Carlo • Reference frame is a triplet of unit vectors • Rotation of the frame is done in two steps.

5. Euler and Quaternion • The only difference between today programs is the way we handle these rotations • Euler angles rotation Quaternionsalgebra

6. 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).

7. 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

8. 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.

9. 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).

10. Monte Carlo flow chart

11. 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.

12. Move • The photon is moved a distance Ds, the new coordinates of the photon are

13. Drop • Reduction of photon weight (W) • According to material albedo • albedo = ms/(ms+ ma) • absorbed = W*(1-albedo)

14. 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.

15. Referencing to scattering plane • Done as in meridian plane Monte Carlo • Multiply the Stokes vector by matrix R

16. Scatter cnt • The interaction with a spherical particle is achieved by multiplying the Stokes vector with a scattering matrix M(q)

17. 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).

18. 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.

19. 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).

20. 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

21. 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).

22. Update of direction cosines • This is done as in standard Monte Carlo • q and f If |uz|≈ 1

23. Update of direction cosines If |uz|≠ 1

24. 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.

25. 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

26. 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

27. Boundaries –R(e) • This is the ONLY reason we need the vector u,v, andw

28. Boundaries –R(e) • This rotation is about the direction of propagation of the photon, i.e. the axis u • Resulting position of v

29. Boundaries –R(y) • A second rotation of an angle y about the Z axis • Put the photon reference frame in the detector reference frame

30. Boundaries –R(y) • Transmission • Reflection

31. Boundaries End • Stokes vector in the detector frame of reference • We need vector u,v, andwto obtain these angles

32. 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

33. 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

34. 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.

35. Vector rotation • First the unit vector v is rotated about the vector u by an angle f

36. 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

37. Vector rotation • Second u is rotated about the newly generated v by an angle q.

38. 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

39. Rotations • These steps are repeated for every scattering event. • At the boundaries the last aligning rotations are the same as in Euler Monte Carlo.

40. 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.

41. 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.

42. 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

43. How to run the code • Download the code • http://faculty.cua.edu/ramella/MonteCarlo/index.html

44. How to run the code • Download the code • http://faculty.cua.edu/ramella/MonteCarlo/index.html

45. Download the code Programs

46. Make Command line Compilation/linking step

47. 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

48. Memory usage • Intel (64-bit) – Macbook pro 2.3Ghz I7 8GB RAM • Real mem 860 KB • Virtual mem 17 MB

49. 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

50. Inside the code GNU Licence