10 views (last 30 days)

Show older comments

Hello,

I'm working on some some Monte Carlo code that models photons moving and scattering in e.g. air or water. I have code that works, but I'd like to speed it up if possible.

I have a pair of coordinates that correspond to two locations where each photon scattered. A vector drawn between these two points traces the path the photon takes, and I already know that this path intercepts my "target", which is parallel to the z-plane and at a user-specified distance. The locations are called hits{1} and hits{2}, the photon moves from {1} to {2}. The columns 1 to 3 in these cells represent the x, y, and z coordinates, respectively.

To find the coordinates of where the photon paths intercept the target, I use the symbolic tool box function syms and solve. First I define the parametric line equations for all the photon paths:

syms t

line = hits{1}(:,1:3) +t*(hits{2}(:,1:3) - hits{1}(:,1:3)); % Line equations

Then I use solve to find the value of t when the z-coordinate is equal to the specified distance at which the target is located ("t_target"). This value of t is then substituted into the line equation using the subs function to get the (x,y) coordinates for each intercept point.

I have got code that does this successfully, however, so far I can only get it to solve each line equation one at a time in a for loop:

coordinates = sym(zeros(size(line,1),3)); % Empty symbolic array for intercept coordinates

weights = NaN(size(line,1),1); % Empty (NaN) array for the corresponding hit weights

t_target = NaN(size(line,1),1); % Empty (NaN) array for values of parameter t

for counter = 1:size(t_target,1)

t_target(counter,:) = solve(line(counter,3) == trgt_dist, t);

coordinates(counter,:) = subs(line(counter,:),t,t_target(counter,:));

end

Would it be faster to vectorise the process and do away with the for loop, effectively running "solve" on all lines at once?

I tried:

t_target(1:end,:) = solve(line(1:end,3) == trgt_dist, t);

but get an error because the r.h.s. yields an empty syms array.

Any help would be appreciated!

David Goodmanson
on 7 Jul 2020

Edited: David Goodmanson
on 7 Jul 2020

Hi Jonathan,

It's not totally clear if the geometry is such that there is a separate t value for each row of the hits{1} and hits{2} matrices. Assuming that's the case, then the line equations are (with some abbreviations for hits{1} and hits{2})

h1 + t.*(h2-h1)

where h1 and h2 are nx3 and t is nx1. The z coordinate is supposed to be t_target which is also nx1 so

h1(:,3) + t.*(h2(:,3)-h1(:,3)) = t_target

with the solution

t = (t_target -h1(:,3))./(h2(:,3)-h1(:,3))

xyzresult = h1 + t.*(h2-h1)

xyresult = xyzresult(:,1:2)

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!