How do I check a each element of a 2 x n matrix for values to remove?

I am trying to simulate a fluorescing vertical 1-D sample, emitting light from -pi/2 to pi/2, passing through a single thin lens, and then striking a detector. I have gotten this to work when the angle is zero (light is traveling horizontally and parallel to other light rays), but when I try to add the angle, it gets bad. To statistically model the process, I have made random number generation loop that creates a random position on the sample (r) and a random angle (r'), and vector v=(r,r'), with r = zero at the center and the sample going up h/2 and down -h/2.

I tried setting up a loop that would check where the ray is located (r) and direction it is aiming (r'), and then check if it is ouside of the lens. For example, diameter lens = Diam= 10, object distance = s_o=100, v=(2,-0.8), the condition would be: v((s_o*tan(v(2,j)) > (Diam/2)+(v(1,j))))=[ ]. I am doing this because it is necesary to omit values that are outside of the lens focal length ratio and diameter, f/.

This is not working. I have tried multiple ways including changing the for loop and if loops to multiple individual loops and even made them switch-case for a while.

function [I,S] = RayTraceAngle(s_i,s_o,f,Diam,L_sample,L_detector,n_pix,n)
%This function receives the optical setup information and determins
%statistical probability of the distribution of photons on a sensor.
%   s_i is the distance from lens to image plane
%   s_o is the distance from lens to object plane
%   f is the focal length of the lens
%   L_sample is the sample height
%   L_detector is the detector height
%   n_pix is the number of vertical pixels.
%   Diam is the lens diameter
%   n is the number of random values or iterations. Reccomend at least
%   125,000
format long
L=[1-s_i/f s_o*(1-s_i/f)+s_i;-1/f 1-s_o/f];
r1 = zeros(2,n); %allocate space for initial vector
r2 = zeros(2,n);%allocate space for final vector
h_pix=L_detector/n_pix;%height of pixel
for i = 1:n
    %%%%%%%%%%%%%begin random position generation%%%%%%%%%%%%%
    a=randi([-1,1],1,1); %generating a integer between -1 and 1
    while a ==0         %excluding zero
    b=randi([-1,1],1,1);%generating another number if a is zero
        a=a+b;          %adding new random number to a
    end
    c=rand(1)*a;    %multiplying random number between 0 and 1 by + or -
    rad_samp = L_sample/2;  %distance from optical axis to edge of sample
    r=rad_samp*c;   %gives position from optical axis to point on sample
    %%%%%%%%%%%%%end random position generation%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%begin random angle generation%%%%%%%%%%%%%
    d=randi([-1,1],1,1); %generating a integer between -1 and 1
    while d ==0         %excluding zero
    e=randi([-1,1],1,1);    %generating another number if a is zero
        d=d+e;          %adding new random number to d
    end
    r_prime=rand(1)*d*pi/2;    %multiplying random number between 0 and 1 by + or -
    v=[r;r_prime];  %the initial vector
   
    
    
    [row,col]= size(v);
    for j=1:1:length(col)
        if (v(1,j) > 0) && (v(2,j)>0)
            v((s_o*tan(v(2,j)) > (Diam/2)-(v(1,j))))=[];
        end
    end
    [row,col]= size(v);
    for j=1:1:length(col)
        if (v(1,j) < 0) && (v(2,j)>0)
            v((s_o*tan(v(2,j)) > (Diam/2)-(v(1,j))))=[];
        end
    end
    [row,col]= size(v);
    for j=1:1:length(col)
        if (v(1,j) > 0) && (v(2,j)<0)
            v((s_o*tan(v(2,j)) > (Diam/2)+(v(1,j))))=[];
        end
    end
    [row,col]= size(v);
    for j=1:1:length(col)
        if (v(1,j) < 0) && (v(2,j)<0)
            v((s_o*tan(v(2,j)) > (Diam/2)+(v(1,j))))=[];
        end
    end

end
    %%%%%%%%%%%%%end random angle generation%%%%%%%%%%%%
    [row,col]= size(v);
    for i = 1:1:length(col)
        r1(:,i)=v;  %initial vector
    r2(:,i)=L*v;    %the found vector
    
I = r2;
S=r1;
edges = -(L_detector/2):h_pix:(L_detector/2);   %creates the edges for the histogram bins with pixel height
histogram(I(1,:),edges)     %creates histogram which includes Left bin edge but not right
end

Answer 1

[Adam Wyatt]

Is your lens a spherical lens, in which case there is a very simple solution (Google intersection of line with sphere: Wiki).

Don't use for-loops - create your matrix of rays using vectorised commands.

Then solve for the intersection of the rays with the lens - this is a quadratic equation with 3 possible sets of solutions:

1. Result is complex - ray misses lens
2. Only a single root - ray is tangent to surface
3. Two roots - ray travels through sphere

You can simply remove any rays for 1 & 2. Then for those with 3 solutions, you can select those that fall within the size of the lens, and select the correct solution (shortest or longest path length depending on whether the lens in concave or convex).

If you want a quick and dirty approach - calculate the intersection of the ray with a plane through the centre of the lens and then select points that fall within the silhouette of the lens on that plane. This will throw out some of the rays that would hit the very edges of the lens, but that's probably fine.

indx = ((x.^2 + y.^2)<R^2)

where indx is a logical array of which rays to accept, x & y are the co-ordinates of the intersection of the ray with the plane and R is the radius of the lens.

[PSilver]

Thank you for the info. I changed my for loops to vectors and it is a lot cleaner. I also used your indx method. Worked great. Thanks again!