How to solve this integral equation with integral2 or Monte Carlo methods?

I am trying to plot this equation as a function of lambda. For that I generating PPP(Poisson Point Process) with density lambda. r0 is radial distance of the nearest point and ri's are those point that lie in ring with inner radius r0 and outer radius rfov. Thus I want to generate realization of PPP with density lambda, compute this integral for each lambda. beta, psifov, h is known.

Answer 1

[HIMANSHU GAUTAM]

clc;
clear all;
% close all;
%% Channel parameters
% psi_half is semi angle of LED.
psi_half = pi/180*60; 
% m is Lambert’s mode number expressing directivity of the source beam.
m = -log10(2)/log10(cos(psi_half));
% psi_FOV is field of view angle of Receiver(Photo detector).
% * From the constant radiance theorem (also known as Etendue limit), the FOV of the receiver
%   is related to the collection area of the lens Acoll and the photodetector area.
%   A_coll*sin(psi_FOV/2)<=A_PD
%
%  * Gain of optical (G_C) concentrator is related to FOV angle as:
%                
%       G_C = mu_C^2/(sin(psi_FOV))^2 ; 0<=theta_rx<=psi_FOV
%           =0........................; theta_rx> psi_FOV
%
% * It is apparent that the concentrator gain increases when the FOV is reduced.
% * Constraint on FOV angle: psi_FOV<pi/2 .
psi_FOV = pi/3;
% mu_C is refractive index of optical concentrator lens. It's value lie in
% range [1 2] for VLC.
mu_C = 1.5;
% G_F is the gain of blue filter.
G_F = 1;
% For 20 MHz Noise Power is:
noisepower = 1e-20;
P_tx = 1;
% Efeective area of photodiode is represented by A_PD
A_PD = .01;
height = 2.5; % in units of meter.
a=5; % considering square room of area 20*20 meter^2.
area=4*a^2;
%% Implementation of Analytical results.
alpha = (m+1)*height^(m+1)/(2*pi);
K = alpha^2;
%s = sym('s','real');
r_fov_anlytical = height*tan(psi_FOV);
d_fov_anlytical = height/cos(psi_FOV);
beta = 2*(m+3);
lambdaaxis_analytical= .05:.08:1.1;%10.^(-2:.3:1); % in units of per meter^2.
avgSINR=zeros(1,length(lambdaaxis_analytical));
for lambdaIdx=1:length(lambdaaxis_analytical)
 %% generation of optical attocels
    NoOfTxs=poissrnd(lambdaaxis_analytical(lambdaIdx)*area);
    locationattocels=unifrnd(-a,a,2,NoOfTxs);
    heightOfTxs=ones(1,NoOfTxs)*height;
    locationattocels=[locationattocels;heightOfTxs];
    locationMs=[0,0]; % Location of Mobile Station.
    %% Distance Calculation
    distanceBS2MS_analytical=sqrt((locationattocels(1,:)-locationMs(1)).^2+(locationattocels(2,:)-locationMs(2)).^2+(locationattocels(3,:)-0).^2);
    radius_vec = sqrt((locationattocels(1,:)).^2+(locationattocels(2,:)).^2);
    %% Throwing out optical transmitter that doen't lie in FOV cone.
     indices_analytical = find(distanceBS2MS_analytical>d_fov_anlytical);
     r_i_vec = radius_vec;
     indices_for_rlessthan_r_0 =  find(r_i_vec>r_fov_anlytical);
     r_i_vec(indices_for_rlessthan_r_0) = []; %removing those BS that have distance less than d_0.
     distanceBS2MS_analytical(indices_analytical) = []; % Taking only those base stations that fall in FOV radius.
    %% Finding nearest BS.
    [~,iselecteddBS]=min(distanceBS2MS_analytical); % Selecting nearest basestations.
    d_0=(distanceBS2MS_analytical(iselecteddBS)).^(-beta/2);
    r_0 = min(r_i_vec);% radius of point where nearest bs is located.
    %% Using Monte Carlo integration method.
    
    r_vec = r_0:.01:r_fov_anlytical;
    inner_integrand = @(x)(1-exp((-x.*(r_vec.^2+height.^2).^(-beta/2))./(d_0))).*r_vec;
    integral1 = trapz(r_vec,(1-exp((-x.*(r_vec.^2+height.^2).^(-beta/2))./(d_0))).*r_vec);%sum((1-exp((-s.*(r_vec.^2+height.^2).^(-beta/2))./(d_0))).*r_vec);
   
end
plot(lambdaaxis_analytical,avgSINR,'b');    
    
    
    
    
    
    
    
    

This is the code I tried but it didn't gave the correct result. I wanna do it by integral2. Additionally if r_0 is empty then I wnt to make avg sinr zero.