Radiation view factor using Monte Carlo

I have the following code for the calcualtion of view factor from a planar patch to a sphere in a frontal configuration using Monte Carlo method for random point generation. When comparing my approximation with the result from analytical expression of view factor my result is way off. Could someone tell me any possible fix, maybe something to do with point generation or with cosines calculation? Thanks in advance

%patch to sphere frontal%
clear all
clc

%geometrical parameters%
R=100;%sphere radius%
H=200;%distance patch-spehere%
h=H/R;
W=5;
L=10;
A1=W*L;%patch area%

%iterations%
k=10;
N=2^k;

%gpoint generation 1-patch,2-sphere%
x1=zeros(N,1);
y1=-W/2+W*rand(N,1);
z1=-L/2+L*rand(N,1);
p1=[x1, y1, z1];

theta=linspace(0,pi(),N)';
phi=linspace(0,2*pi(),N)';
x2 = H+R .* sin(phi).* cos(theta);
y2 = R .* sin(phi) .* sin(theta);
z2 = R .* cos(phi);
p2 = [x2, y2, z2];

%cosines and distance, cosines for angle between normal vector to each surface anad ray between p1-p2%
d=p2-p1;%distance%
A2=4*pi()*R^2;%sphere area%
centro=[H,0,0]';%spehere center%
dir_esfera=p2-centro';
norm_esfera=dir_esfera;
vx=[1,0,0];
norm_patch=repelem(vx,N,1);

for i = 1:N
    r12(i)=norm(d(i));
    cosb1(i)=dot(d(i),norm_patch(i),2)/(norm(norm_patch(i))*r12(i));
    cosb2(i)=dot(d(i),norm_esfera(i),2)/(norm(norm_esfera(i))*r12(i));
    cosb2_valid=fillmissing(cosb2,'constant',rand(1));
    f(i)=(cosb1(i)*cosb2_valid(i))/(pi()*r12(i)^2);
end

fi=sum(f);
F12=(A2/N)*fi %view factor%

F12 = 0.7057

%analytical%
F12a=1/h^2

F12a = 0.2500

Answer 1

[Alan Stevens]

I think you probably need cos(theta), rather than theta, to be linearly spaced.

[Hugo Blanco Quintana]

Thanks for your reply Alan

I tried that but now I keep getting negative values which are also way off

This is the current state of the points generation, any furteher ideas?

thanks