Using 'solve' on a vector of line equations?

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!

Answer 1

[David Goodmanson]

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)

[Jonathan]

Hello David,

Many thanks for that. Yes, you are correct that "there is a separate t value for each row of the hits{1} and hits{2}".

I've sanity-checked it to confirm it gives the expected results and same results as before, which it does. It is much faster than my previous approach!

Best wishes,

Jonathan