Rotate Existing Vector

Question

Hello and apologies if the title of the question is not very precise.

Question: I am reading the document talking about the simulation of photons in tissues using a Monte Carlo simulation. The exact title is "MCNP - A general Monte Carlo N-Partcle Transport Code".

Link to MCNP - A general Monte Carlo N-Partcle Transport Code

When photons are scattered, a new direction for the photon is sampled. When the medium is isotropic the direction is random so the math for this case are not hard, but when the medium is not isotropic, the angle of deflection needs to be computed from a function such as the Heney-Greenstein function.

Anyway to be short, this function returns a $\cos(\theta)$ which is the angle of deflection between the existing photon direction and the new desired direction. In the paper, this term is called $\mu_{lab}$. In the document I am referring to, they also compute a $\phi$ angle by sampling two random uniformly distributed variables, which are inscribed in the unit disk. These are called $\epsilon_1$ and $\epsilon_2$ in the paper.

So the formula they use in this paper (see reference below page 2-38) to compute the new photon direction (using the three aforementioned variables) is:

And where $u_0 v_0 w_0$ are the coordinates of the "old photon direction" and $u v w$ the new photon direction after scattering. My problem, is that I have no idea how they derived this formula.

So the way I understand how this can be done is by "reconstructing" a coordinate system in which the old direction is the z unit vector, and express the coordinates of the new direction within this frame? Is that correct. But could someone put me on the right track so that I can derive this formula (understand how they get there)? I don't need a full answer, I am happy to make an effort, I just need someone to put me not the right track (and I will publish the answer when I have one).

It would be fantastic. Your help is greatly appreciated.

EDIT:

Following John's comments, if we consider that these formulas are an application of the Rodrigues' rotation formula where the last term of the equation is $0$ because $k.v = 0$ (the two vectors are orthogonal), we have for this formula:

$R = \cos(\theta) v + \sin(\theta) (k \times v)$

So using the equations provided in the paper, for u, v, and w (the coordinates of the rotate vector), I found for k:

$k_x = { - { \epsilon_1 v_0 } \over { \sqrt { (\epsilon_1^2 + \epsilon_2^2)(1 - w_0^2) } } } = -{{v_0 \cos(\phi) } \over {\sqrt { (1 - w_0^2) } }}$,

$k_y = { { \epsilon_1 u_0 } \over { \sqrt { (\epsilon_1^2 + \epsilon_2^2)(1 - w_0^2) } } } = {{u_0 \cos(\phi) } \over {\sqrt { (1 - w_0^2) } }}$,

$k_z = { { \epsilon_2 } \over { \sqrt { (\epsilon_1^2 + \epsilon_2^2)(1 - w_0^2) } } } = {{\sin(\phi) } \over {\sqrt { (1 - w_0^2) } }}$,

The problem is that I can't find "how" they found this k? I looked into the Gram-Schmidt method but when I develop the terms using this method, I don't end up with the same result. So I have no idea which method they use to construct k. I searched, I really did, nothing comes to my mind, so if someone could give it a little push, it would be greatly appreciated.

Answer 1

First, the dangling $o$ in the numerator for the expression for $u$ should be $v_o$. (Otherwise, the angle between $(u_o,v_o,w_o)$ and $(u,v,w)$ would not be $\theta$.) There is a mistake in the link you provide.

Start with unit vector $\mathbf{\hat{x}}\equiv(u_o, v_o, w_o)$. We wish to compute all unit vectors $\mathbf{\hat{y}}$ such that the angle between $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ is $\theta$, i.e., $\mathbf{\hat{x}}\cdot\mathbf{\hat{y}}=\cos\theta$. It is easy to see that $\mathbf{\hat{y}}$ is a unit vector such that $\mathbf{\hat{x}}\cdot\mathbf{\hat{y}}=\cos\theta$ if and only if $$ \mathbf{\hat{y}} = \cos\theta \mathbf{\hat{x}} + \sin\theta \mathbf{\hat{z}} $$ for some unit vector $\mathbf{\hat{z}}$ orthogonal to $\mathbf{\hat{x}}$. (Proof: Suppose $\mathbf{\hat{y}} = \cos\theta \mathbf{\hat{x}} + \sin\theta \mathbf{\hat{z}}$ where $\mathbf{\hat{x}}\cdot\mathbf{\hat{z}}=0$. Then $\mathbf{\hat{x}}\cdot\mathbf{\hat{y}}=\cos\theta$. Conversely, suppose $\mathbf{\hat{y}}$ is a unit vector with $\mathbf{\hat{x}}\cdot\mathbf{\hat{y}}=\cos\theta$. Then $$ \mathbf{\hat{y}} = \cos\theta\mathbf{\hat{x}} + \frac{\mathbf{\hat{y}}-\cos\theta\mathbf{\hat{x}}}{\| \mathbf{\hat{y}}-\cos\theta\mathbf{\hat{x}}\|}\| \mathbf{\hat{y}}-\cos\theta\mathbf{\hat{x}}\|. $$ But $\|\mathbf{\hat{y}} - \cos\theta\mathbf{\hat{x}}\| = \sqrt{(\mathbf{\hat{y}} - \cos\theta\mathbf{\hat{x}})\cdot(\mathbf{\hat{y}} - \cos\theta\mathbf{\hat{x}})} = \sqrt{1-2\cos\theta + \cos^2\theta} = \sin\theta$.)

Next, we characterize all unit vectors orthogonal to $\mathbf{\hat{x}}$. Since we are in $\mathbb{R}^3$, the set of vectors orthogonal to $\mathbf{\hat{x}}$ is a 2-dimensional dimensional vector space. If we find two unit vectors $\mathbf{\hat{h}_1}$ and $\mathbf{\hat{h}_2}$ that are orthogonal to one another and to $\mathbf{\hat{x}}$, then every unit vector orthogonal to $\mathbf{\hat{x}}$ can be written as $$ \alpha \mathbf{\hat{h}_1} + \beta \mathbf{\hat{h}_2} $$ for some $\alpha$, $\beta$ where $\alpha^2+\beta^2=1$. (The proof is easy.)

Here are two vectors orthogonal to $\mathbf{\hat{x}}\equiv(u_o,v_o,w_o)$ and to each other: \begin{eqnarray} (-v_o,u_o,0) \\ (u_o w_o, v_o w_o, -u_o^2 - v_o^2). \end{eqnarray} Lots of other vectors could have been chosen. I found the first of these two vectors by examining the formulas you provided. I computed the second as the vector cross product of $\mathbf{\hat{x}}$ and the first vector. (This ensures mutual orthogonality.)

Each of these two vectors has length $\sqrt{u_o^2 + v_o^2} = \sqrt{1-w_o^2}$. We can normalize each to obtain: \begin{eqnarray} \mathbf{\hat{h}_1} &=& \frac{1}{\sqrt{1-w_o^2}}(-v_o,u_o,0) \\ \mathbf{\hat{h}_2} &=& \frac{1}{\sqrt{1-w_o^2}}(u_o w_o, v_o w_o, -(1- w_o^2)). \end{eqnarray}

Putting all this together, we can compute $\mathbf{\hat{y}}$ as \begin{eqnarray} \mathbf{\hat{y}}\equiv(u,v,w) &=& \cos\theta\ \ (u_o,v_o,w_o) + \alpha\sin\theta\frac{1}{\sqrt{1-w_o^2}}(-v_o,u_o,0) \\ & & + \beta\sin\theta\ \ \frac{1}{\sqrt{1-w_o^2}}(u_o w_o, v_o w_o, -(1-w_o^2)). \end{eqnarray} Now just set \begin{eqnarray} \alpha &=& \frac{\xi_2}{\sqrt{\xi_1^2 + \xi_2^2}} \\ \beta &=& \frac{\xi_1}{\sqrt{\xi_1^2 + \xi_2^2}} \\ \cos\theta &=& \mu_{lab} \\ \sin\theta &=& \sqrt{1-\mu_{lab}^2} \end{eqnarray} and we have the formulae you gave.

Answer 2

Well, there are two parts.

Let's write this in vectors. Your expression says

$$ \mathbf v = \mu_\text{lab} \mathbf v_0 + \mathbf h $$ where $\mathbf h$ is a complicated expression. But let's look at just this much. Imagine that $\mathbf v_0$ points to the north pole; this expression says that you're going to build $\mathbf v$ to have some northerly component, namely $\cos \theta$ times the old northerly component. So perhaps a better way to write that is

$$ \mathbf v = (\cos \theta)\mathbf v_0 + (\sin \theta) \mathbf k $$ where $\mathbf k$ is going to be a unit vector perpendicular to $\mathbf v_0$.

How do they determine $\mathbf k$? It looks to me as if it's basically Rodrigues' formula, disguised with lots of greek letters. The first square root (which appears in all three formulas) is just $\sin \theta$, for instance. The terms like $u_0 w_0$ come from the cross product. Try looking up R's formula, and you should see where all this comes from.

If you can find a copy of the 3rd edition of Computer Graphics: Principles and Practice, there's a pretty extensive derivation of R's formula in the chapter on transformations in 3 dimensions. (Disclaimer: I'm an author. There are surely other books with good derivations as well; this just just the one I happen to know.)