How to simulate RMS delay distribution?

Hello everyone, my name is Duy and i'm doing a project about VLC.

I tried the code in page 465 from the book: "Optical Wireless Communications: System and Channel Modelling with MATLAB" but it didn't work. I dont know if you can help me or not, but I would be very greatful if somebody could. I really really appreciate that. Thank you so much.

p/s: Here are the code :D

%%
C=3e8*1e-9;
%time will be measured in ns in the program
theta=70;
%semi-angle at half power
m=-log10(2)/log(cosd(theta));
%Lambertian order of emission
P_LED=20;
%transmitted optical power by individual LED
nLED=60;
%number of LED array nLED*nLED
P_total=1;
%Total transmitted power
Adet=1e-4;
%detector physical area of a PD
rho=0.8;
%reflection coefficient
Ts=1;
%gain of an optical filter;
index=1.5;
%refractive index of a lens at a PD
FOV=60;
%Fov of a receiver
G_Con=(index^2)/(sind(FOV).^2);
%Gain of an optical concentrator

%%
lx=5; ly=5; lz=3-0.85;
%Room dimension in meter
[XT,YT,ZT]=meshgrid([-lx/4 lx/4], [-ly/4 ly/4],lz/2);
%position of transmitter
Nx=lx*10; Ny=ly*10; Nz=round(lz*10);
%number of grid in each surface
dA=(lz*ly)/(Ny*Nz);
%calculation grid area
x=linspace(-lx/2,lx/2,Nx);
y=linspace(-ly/2,ly/2,Ny);
z=linspace(-lz/2,lz/2,Nz);

%%
TP1=[-lx/4 -ly/4 lz/2];
TP2=[lx/4 ly/4 lz/2];
TP3=[lx/4 -ly/4 lz/2];
TP4=[-lx/4 ly/4 lz/2];
%Position of transmitter 1, 2, 3, 4 respectively
TPV=[0 0 -1];
%Transmitter position vector
RPV=[0 0 1];
%Receiver position vector

%%
WPV1=[1 0 0];
WPV2=[0 1 0];
WPV3=[-1 0 0];
WPV4=[0 -1 0];
%Position of wall 1, 2, 3, 4 respectively
delta_t=1/2;
%time resolution in ns, use in the form of 1/2^m

for ii=1:Nx
    for jj=1:Ny
        RP=[x(ii) y(jj) -lz/2];
        t_vector=0:30/delta_t; %time vector in ns
        h_vector=zeros(1,length(t_vector));
        %receiver position vector
        %LOS channel gain
        D1=sqrt(dot(TP1-RP, TP1-RP));
        cosphi=lz/D1;
        tau0=D1/C;
        index=find(round(tau0/delta_t)==t_vector);
        if abs(acosd(cosphi))<=FOV
            h_vector(index)=h_vector(index)+(m+1)*Adet.*cosphi.^(m+1)./(2*pi.*D1.^2);
            h_vector(index)= fliplr(h_vector(index)) + fliplr(flipud(h_vector(index))) + flipud(h_vector(index));
        end
        %%reflection from wall 1
        count=1;
        for kk=1:Ny
            for ll=1:Nz
                WP1=[-lx/2 y(kk) z(ll)];
                D1=sqrt(dot(TP1-WP1, TP1-WP1));
                cos_phi=abs(WP1(3)-TP1(3))/D1;
                cos_alpha=abs(TP1(1)-WP1(1))/D1;
                
                D2=sqrt(dot(WP1-RP, WP1-RP));
                cos_beta=abs(WP1(1)-RP(1))/D2;
                cos_psi=abs(WP1(3)-RP(3))/D2;
                tau1=(D1+D2)/C;
                index=find(round(tau1/delta_t)==t_vector);
                if abs(acosd(cos_psi))<=FOV
                    h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*cos_phi^m*cos_alpha*cos_beta*cos_psi/(2*pi^2*D1^2*D2^2);
                end
                WP2=[x(ii) -ly/2 z(ll)];
                RP2=[x(ii) y(kk) lz/2];
                D21=sqrt(dot(TP2-WP2, TP2-WP2));
                cos_phi2=abs(WP2(3)-TP2(3))/D21;
                cos_alpha2=abs(WP2(2)-TP2(2))/D21;
                
                D22=sqrt(dot(WP2-RP, WP2-RP));
                cos_beta2=abs(WP2(2)-RP2(2))/D22;
                cos_psi2=abs(WP2(3)-RP2(3))/D22;
                tau2=(D21+D22)/C;
                index=find(round(tau2/delta_t)==t_vector);
                if abs(acosd(cos_psi))<=FOV
                    h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*cos_phi2^m*cos_alpha2*cos_beta2*cos_psi2/(2*pi^2*D21^2*D22^2);
                end
                WP3=[lx/2 y(kk) z(ll)];
                RP3=[x(ii) y(kk) lz/2];
                D31=sqrt(dot(TP3-WP3, TP3-WP3));
                cos_phi3=abs(WP3(3)-TP3(3))/D31;
                cos_alpha3=abs(WP3(1)-TP3(1))/D31;
                
                D32=sqrt(dot(WP3-RP, WP3-RP));
                cos_beta3=abs(WP3(1)-RP3(1))/D32;
                cos_psi3=abs(WP3(3)-RP3(3))/D32;
                tau3=(D31+D32)/C;
                index=find(round(tau3/delta_t)==t_vector);
                if abs(acosd(cos_psi))<=FOV
                    h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*cos_phi3^m*cos_alpha3*cos_beta3*cos_psi3/(2*pi^2*D31^2*D32^2);
                end
                WP4=[x(ii) ly/3 z(ll)];
                RP4=[x(ii) y(kk) lz/2];
                D41=sqrt(dot(TP4-WP4, TP4-WP4));
                cos_phi4=abs(WP4(3)-TP4(3))/D41;
                cos_alpha4=abs(WP4(2)-TP4(2))/D41;
                
                D42=sqrt(dot(WP4-RP, WP4-RP));
                cos_beta4=abs(WP4(2)-RP4(2))/D42;
                cos_psi4=abs(WP4(3)-RP4(3))/D42;
                tau4=(D41+D42)/C;
                index=find(round(tau4/delta_t)==t_vector);
                if abs(acosd(cos_psi))<=FOV
                    h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*cos_phi4^m*cos_alpha4*cos_beta4*cos_psi4/(2*pi^2*D41^2*D42^2);
                end
                count=count+1;
            end
        end
        
      t_vector=t_vector*delta_t;
      mean_delay(ii,jj)=sum((h_vector).^2.*t_vector)/sum(h_vector.^2);
      Drms(ii,jj)=sqrt(sum((t_vector-mean_delay(ii,jj)).^2.*h_vector.^2)/sum(h_vector.^2));
    end
end

surfc(x,y,Drms);
%surf(x,y,mean_delay);
axis([-lx/2 lx/2 -ly/2 ly/2 min(min(Drms)) max(max(Drms))]);

[Busra AYDIN]

hello, it may be late for your question but this is fix of the code for the other researchers:

clc

clear all

%%

C=3e8*1e-9;

% time will be measured in ns in the program

theta=70;

% semi-angle at half power

m=-log10(2)/log10(cosd(theta));

% Lambertian order of emission

% P_LED=20;

% %transmitted optical power by individual LED

% nLED=60;

% number of LED array nLED*nLED

P_total=1;

% Total transmitted power

Adet=1e-4;

% detector physical area of a PD

rho=0.8;

% reflection coefficient

Ts=1;

% gain of an optical filter; ignore if no filter is used

index=1.5;

% refractive index of a lens at a PD; ignore if no lens is used

FOV=60;

% FOV of a receiver

G_Con=(index^2)/(sind(FOV).^2);

% gain of an optical concentrator; ignore if no lens is used

%%

%%%%%%

lx=5; ly=5; lz=3-0.85;

% room dimension in metre

%[XT,YT,ZT]=meshgrid([−lx/4 lx/4],[−ly/4 ly/4],lz/2);

% position of Transmitter (LED);

Nx=lx*3; Ny=ly*3; Nz=round(lz*3);

% number of grid in each surface

dA=lz*ly/(Ny*Nz);

% calculation grid area

x=-lx/2:lx/Nx:lx/2;

y=-ly/2:ly/Ny:ly/2;

z=-lz/2:lz/Nz:lz/2;

%%

% first transmitter calculation

TP1=[0 0 lz/2];

% transmitter position

TPV=[0 0 -1];

% transmitter position vector

RPV=[0 0 1];

% receiver position vector

%%

WPV1=[1 0 0];

WPV2=[0 1 0];

WPV3=[-1 0 0];

WPV4=[0 -1 0];

% wall vectors

delta_t=1/2;

% time resolution in ns, use in the form of 1/2ˆm

for ii=1:Nx+1

for jj=1:Ny+1

RP=[x(ii) y(jj) -lz/2];

t_vector=0:25/delta_t; % time vector in ns

h_vector=zeros(1,length(t_vector));

% receiver position vector

% LOS channel gain

D1=sqrt(dot(TP1-RP,TP1-RP));

cosphi=lz/D1;

tau0=D1/C;

index=find(round(tau0/delta_t)==t_vector);

if abs(acosd(cosphi))<=FOV

h_vector(index)=h_vector(index)+(m+1)*Adet.*cosphi.^(m+1)./(2*pi.*D1.^2);

end

%% reflection from first wall

count=1;

for kk=1:Ny+1

for ll=1:Nz+1

WP1=[-lx/2 y(kk) z(ll)];

D1=sqrt(dot(TP1-WP1,TP1-WP1));

cos_phi=abs(WP1(3)-TP1(3))/D1;

cos_alpha=abs(TP1(1)-WP1(1))/D1;

D2=sqrt(dot(WP1-RP,WP1-RP));

cos_beta=abs(WP1(1)-RP(1))/D2;

cos_psi=abs(WP1(3)-RP(3))/D2;

tau1=(D1+D2)/C;

index=find(round(tau1/delta_t)==t_vector);

if abs(acosd(cos_psi))<=FOV

h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*...

cos_phi^m*cos_alpha*cos_beta*cos_psi/(2*pi^2*D1^2*D2^2);

end

count=count+1;

end

end

%% Reflection from second wall

count=1;

for kk=1:Nx+1

for ll=1:Nz+1

WP2=[x(kk) -ly/2 z(ll)];

D1=sqrt(dot(TP1-WP2,TP1-WP2));

cos_phi= abs(WP2(3)-TP1(3))/D1;

cos_alpha=abs(TP1(2)-WP2(2))/D1;

D2=sqrt(dot(WP2-RP,WP2-RP));

cos_beta=abs(WP2(2)-RP(2))/D2;

cos_psi=abs(WP2(3)-RP(3))/D2;

tau2=(D1+D2)/C;

index=find(round(tau2/delta_t)==t_vector);

if abs(acosd(cos_psi))<=FOV

h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*...

cos_phi^m*cos_alpha*cos_beta*cos_psi/(2*pi^2*D1^2*D2^2);

end

count=count+1;

end

end

%% Reflection from third wall

count=1;

for kk=1:Ny+1

for ll=1:Nz+1

WP3=[lx/2 y(kk) z(ll)];

D1=sqrt(dot(TP1-WP3,TP1-WP3));

cos_phi= abs(WP3(3)-TP1(3))/D1;

cos_alpha=abs(TP1(1)-WP3(1))/D1;

D2=sqrt(dot(WP3-RP,WP3-RP));

cos_beta=abs(WP3(1)-RP(1))/D2;

cos_psi=abs(WP3(3)-RP(3))/D2;

tau3=(D1+D2)/C;

index=find(round(tau3/delta_t)==t_vector);

if abs(acosd(cos_psi))<=FOV

h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*...

cos_phi^m*cos_alpha*cos_beta*cos_psi/(2*pi^2*D1^2*D2^2);

end

count=count+1;

end

end

%% Reflection from fourth wall

count=1;

for kk=1:Nx+1

for ll=1:Nz+1

WP4=[x(kk) ly/2 z(ll)];

D1=sqrt(dot(TP1-WP4,TP1-WP4));

cos_phi= abs(WP4(3)-TP1(3))/D1;

cos_alpha=abs(TP1(2)-WP4(2))/D1;

D2=sqrt(dot(WP4-RP,WP4-RP));

cos_beta=abs(WP4(2)-RP(2))/D2;

cos_psi=abs(WP4(3)-RP(3))/D2;

tau4=(D1+D2)/C;

index=find(round(tau4/delta_t)==t_vector);

if abs(acosd(cos_psi))<=FOV

h_vector(index)=h_vector(index)+(m+1)*Adet*rho*dA*cos_phi^m*cos_alpha*cos_beta*cos_psi/(2*pi^2*D1^2*D2^2);

end

count=count+1;

end

end

t_vector=t_vector*delta_t;

mean_delay(ii,jj)=sum((h_vector).^2.*t_vector)/sum(h_vector.^2);

Drms1(ii,jj)=sqrt(sum((t_vector-mean_delay(ii,jj)).^2.*h_vector.^2)/sum(h_vector.^2));

end

end

surf(x,y, Drms1);

% surf(x,y,mean_delay);

% plot(t_vector*delta_t, h_vector)

% save Drms_1LED mean_delay Drms x y z theta m P_total rho FOV

TP1;

axis([-lx/2 lx/2 -ly/2 ly/2 min(min(Drms1)) max(max(Drms1))]);