View Factor Calculation with Monte Carlo Method

Hi! I am trying to find view factors between a disk and a surface shown in the fig. that i added. So i ned 2 different random numbers for that one for teta angle and the other for azimut angle. First i am checking the teta angle is it is in the range if it is in the range and next i 'm gonna check azimut angle if it is in the range or not. If azimut in the range also it is a hit so it will be registried as a hit. Im trying to do this without for loop cause i'm gonna go up to 10^9 for mt random numbers. See NN=10,10^2,10^3,...10^9. My problem is after i checked teta angle i want my azimut calculated with my new teta(azimut depend on teta) .That means if c==d i want to throw these elements from the teta and my teta_new will builted up without them. but i dont know how can i throw these elements and built my new teta which is teta_new. If i can throw them i will built up teta_new and azimut_new and so i can find how many hits will be registered. Thanks!!

clear;

a=2;

r=4;

h=1;

TH_max =atan((a+r)/h);

TH_min=atan((a-r)/h);

NN1=10000;

RN1=rand(1,NN1)./max(rand(1,NN1));

TH=asin(sqrt(RN1));

x1=(h.^2.*(tan(TH)).^2+a.^2-r.^2)./(2.*a);

azimut_max=atan(sqrt(h.^2.*(tan(TH).^2-x1.^2)./x1));

c=TH(1,:)<TH_max ; d=TH_min<TH(1,:);

TH_new=???

NN2=sum(c==d);

RN2=rand(1,NN2)./max(rand(1,NN2));

azimut=RN2*2*pi;

azimut_max_new=atan(sqrt(h.^2.*(tan(TH_new).^2-x1.^2)./x1));

REG=sum(azimut<=azimut_max_new);

VF = REG/NN2;

H=h/a;

R=r/a;

Z=1+H.^2+R.^2;

VF_exact=0.5.*(1-((Z-2.*R.^2)./sqrt(Z.^2-4.*R.^2)));