Anyone knows what's wrong with this piece of Matlab code? It's outputting nonsense as I can manually find a vector that actually intersects with the rectangle but the program isn't adding it up. If it did, then it should show up in the prob even if it's 1/300000.
function [prob, P] = MCProbSeq(d, alpha, dir)
%MCPROBSEQ Generates the probability of hitting a rectangular target
% at a given distance d with apex angle alpha, the axis of the cone
% pointing at direction of vector dir, with a Monte Carlo method.
% The random direction is assumed to be uniformly distributed within the
% cone.
dist = [0 0 d]';
a = dir/sqrt(sum(dir.^2));%Unit vector
phi = acos(dot(a, [0 0 1]'));
if alpha >= pi
error('Apex angle too big');
elseif d < 0
error('d must be positive scalar')
end
if (pi/2 < phi < 3*pi/2 && d == 0) %Point blank range
prob = 1;
return
end
%Rectangular approximation based on a person with 1.8 m height.
len = 1.8;
wid = 0.88/1.8;
apex = alpha/2;
n = [0 0 1]';%Normal vector of the plane
%Using p_0 = [0 0 0]'
A = null(a');
u = A(:,1);%Orthogonal unit vectors
v = A(:,2);
samples = 300000;
t = 2*pi*rand(samples,1);%Uniform distribution
z = cos(apex) + (1-cos(apex))*rand(samples,1);%Uniform distribution
P = [sqrt(1-z.^2).*cos(t) sqrt(1-z.^2).*sin(t) z];%Uniform distribution
count = 0;
for i=1:samples
P(i,:) = (P(i,1)*u + P(i,2)*v + P(i,3)*a)';
if dot(P(i,:)', n) ~= 0 %Is the line parallel?
k = dot(-dist, n)/dot(P(1,:),n);
if k >= 0 %Is it pointing at it?
point = k*P(1,:)' + dist;
if abs(point(1,1)) <= len/2 && abs(point(2,1)) <= wid/2
%Is the point on the rectangle?
count = count + 1;
end
end
end
end
prob = count/length(P(:,1));
end