室内可见光通信一次反射功率分布matlab代码

theta=70; %半功率角

m=-log10（2）/log10（cosd（theta））; %兰博辐射

P_LED=20;

nLED=60;

P_total=nLED*nLED*P_LED;

Adet=1e-4;
%PD面积
rho=0.8; %反射系数

Ts=1;

index=1.5;
%透镜反射
FOV=60;

G_Con=（index^2）/（sind（FOV）.^2）;

lx=5; ly=5; lz=2.15;

Nx=lx*12; Ny=ly*12; Nz=round（lz*12）;

dA=lz*ly/（Ny*Nz）;

x=linspace（0lxNx）;
y=linspace（0lyNy）;
z=linspace（0lzNz）;
[XRYRZR]=meshgrid（xy0）;
%%
%计算反射
TP1=[1.251.252.15];

WPV1=[1 0 0]; %第一面墙

for ii=1:Nx
    for jj=1:Ny
        H_LOS（iijj）=0;
        RP_LOS=[x（ii） y（jj） 0];
        % receiver position vector
        D_LOS=sqrt（dot（RP_LOS-TP1RP_LOS-TP1））;
        cos_LOS_phi = abs（TP1（3）- RP_LOS（3））/D_LOS;
        receiver_angle= acosd（cos_LOS_phi）;
        if receiver_angle<=FOV
         H_LOS（iijj） = H_LOS（iijj）+（m+1）*Adet*cos_LOS_phi^（m+1）/（2*pi*D_LOS^2）;
        end
    end
end
for ii=1:Nx
for jj=1:Ny
RP=[x（ii） y（jj） 0];

h1（iijj）=0; %信道增益

for kk=1:Ny
for ll=1:Nz
WP1=[0 y（kk） z（ll）];

D1=sqrt（dot（TP1-WP1TP1-WP1））;
cos_phi=abs（WP1（3）- TP1（3））/D1;
cos_alpha=abs（TP1（1）- WP1（1））/D1;
D2=sqrt（dot（WP1-RPWP1-RP））;
cos_beta=abs（WP1（1）- RP（1））/D2;
cos_psi=abs（WP1（3）- RP（3））/D2;
if abs（acosd（cos_psi））<=FOV
h1（iijj）=h1（iijj）+（m+1）*Adet*rho*dA*cos_phi^m*cos_alpha*cos_beta*cos_psi/（2*pi^2*D1^2*D2^2）;
end
end
end
end
end
%%

WPV2=[0 1 0];
% 墙2
for ii=1:Nx
for jj=1:Ny
RP1=[x（ii） y（jj） 0];

h2（iijj）=0;

for kk=1:Ny
for ll=1:Nz
WP1=[lx y（kk） z（ll）];
D1=sqrt（dot（TP1-WP1TP1-WP1））;
cos_phi=abs（WP1（3）- TP1（3））/D1;
cos_alpha=abs（TP1（1）- WP1（1））/D1;
D2=sqrt（dot（WP1-RP1WP1-RP1））;
cos_beta=abs（WP1（1）- RP1（1））/D2;
cos_psi=abs（WP1（3）- RP1（3））/D2;
if abs（acosd（cos_psi））<=FOV
h2（iijj）=h2（i