How to simulate Illuminance distribution?

Hi everyone! My name is Alberto and I am going to start a VLC-LED based project in MATLAB.

As you know, LEDs are Lambertian sources, and the project consist in simulate 1,2...n LED sources in the ceiling of a room, to see how much illuminance they can provide.

I want to get something like Figures 3, 4 and 5 in this document: http://nrl.northumbria.ac.uk/712/1/paper_csndsp2010_0506.pdf

I have been trying this morning but results are not as expected.

These are the codes for 1 LED and for 4 LEDs, respectively:

% semi-angle at half power
theta=30;

%Lambertian order of emission
m=-log10(2)/log10(cosd(theta));

%Center luminous intensity y total según número de LEDs
I0=0.73;
I0_total=60*60*I0;

% room dimension in metre
lx=5; ly=5; lz=3;

%the distance between source and receiver plane
h=3;

% position of LED1
XTrans1=2.5; YTrans1=2.5;

% number of grid in the receiver plane (cada 2 centímetros)
Nx=lx*50; Ny=ly*50;
x=0:lx/Nx:lx; 
y=0:ly/Ny:ly;
[XRec,YRec]=meshgrid(x,y);

% distance vector from source 1
Vector_Distancia=sqrt((XRec-XTrans1(1,1)).^2+(YRec-YTrans1(1,1)).^2+h^2);

% vector que contiene el angulo de irradiancia (del LED) para cada posición
coseno_phi=h./Vector_Distancia;

% FLUJO LUMINOSO o ILUMINANCIA (en lux)
E_lux=(I0_total*coseno_phi.^m)./(Vector_Distancia.^2);

meshc(x,y,E_lux);
xlabel('X (m)');
ylabel('Y (m)');
zlabel('Illuminance(lx)');

And here for 4 LEDs:

% semi-angle at half power
  theta=15;

    %Lambertian order of emission
    m=-log10(2)/log10(cosd(theta));

    %Center luminous intensity y total según número de LEDs
    I0=0.73;
    I0_total=60*60*I0;

    % room dimension in metre
    lx=5; ly=5; lz=3;

    %the distance between source and receiver plane
    h=3;

    % position of LED1
    XTrans1=1.25; YTrans1=1.25;
    % position of LED2
    XTrans2=1.25; YTrans2=3.75;
    % position of LED3
    XTrans3=3.75; YTrans3=1.25;
    % position of LED4
    XTrans4=3.75; YTrans4=3.75;

    % number of grid in the receiver plane (cada 5 centímetros)
    Nx=lx*20; Ny=ly*20;
    x=0:lx/Nx:lx; 
    y=0:ly/Ny:ly;
    [XRec,YRec]=meshgrid(x,y);

    % distance vector from source 1
    D1=sqrt((XRec-XTrans1(1,1)).^2+(YRec-YTrans1(1,1)).^2+h^2);
    % distance vector from source 2
    D2=sqrt((XRec-XTrans2(1,1)).^2+(YRec-YTrans2(1,1)).^2+h^2);
    % distance vector from source 3
    D3=sqrt((XRec-XTrans3(1,1)).^2+(YRec-YTrans3(1,1)).^2+h^2);
    % distance vector from source 4
    D4=sqrt((XRec-XTrans4(1,1)).^2+(YRec-YTrans4(1,1)).^2+h^2);

    % vector que contiene el angulo de irradiancia (del LED) para cada posición
    coseno_phi1=h./D1;
    coseno_phi2=h./D2;
    coseno_phi3=h./D3;
    coseno_phi4=h./D4;

    % FLUJO LUMINOSO o ILUMINANCIA (en lux)
    E_lux1=(I0_total*coseno_phi1.^m)./(D1.^2);
    E_lux2=(I0_total*coseno_phi2.^m)./(D2.^2);
    E_lux3=(I0_total*coseno_phi3.^m)./(D3.^2);
    E_lux4=(I0_total*coseno_phi4.^m)./(D4.^2);

    E_lux=E_lux1+E_lux2+E_lux3+E_lux4;

    meshc(x,y,E_lux);
    xlabel('X (m)');
    ylabel('Y (m)');
    zlabel('Illuminance(lx)');

I dont know if you can help me or not, but I would be very greatful if somebody could.

Thank you very much in advance.

[Ivo Resende]

Hi Alberto,

Why did you do: E_lux=(I0_total*coseno_phi.^m)./(Vector_Distancia.^2) instead of E_lux=(I0_total/((Vector_Distancia.^2).*coseno_phi) as it is in the document (equation 3)?

[Bhagyesh Shiyani]

hi evo,

i think both will give you the same answer but multiplying first with I0_total will generate array and dividing array(multiplied Ans) with array(vector_distance) will be more faster for simulation.

[AHAD]

Can i have for 6 led

Answer 1

[Pham Duc Thinh]

Dear Abelto,

Here is my answer, i hope that it is useful for you "

% semi-angle at half power theta=30; %Lambertian order of emission m=-log10(2)/log10(cosd(theta)); %Center luminous intensity y total según número de LEDs I0=0.73; I0_total=60*60*I0; % room dimension in metre lx=5; ly=5; lz=3; %the distance between source and receiver plane h=2.15; % position of LED1 [XT,YT]=meshgrid(2.5,2.5); % So luong mat luoi Nx=lx*50; Ny=Nx; x=linspace(0,lx,50); y=linspace(0,ly,50); [XR,YR]=meshgrid(x,y); % distance vector from source 1 Vector_Distance=sqrt((XR-XT(1,1)).^2+(YR-YT(1,1)).^2+h^2); % vector que contiene el angulo de irradiancia (del LED) para cada posición cos_phi=h./Vector_Distance; % FLUJO LUMINOSO o ILUMINANCIA (en lux) E_lux=(I0_total*cos_phi.^m)./(Vector_Distance.^2); meshc(x,y,E_lux); xlabel('X (m)'); ylabel('Y (m)'); zlabel('Illuminance(lx)');

Hint : Find this Book : Optical Wireless Communications: System and Channel Modelling with MATLAB.pdf

All you need is inside it. Good luck to you ! :)

[Alberto López]

Hi Pham! First of all, thank you very much for your repply. I read that book a month ago and I finally got it, but I dont know why I0 is 0.73 and not other number.

Another thing I want to do and I cant, is to represent the illuminance in the same graphic as X,Y,Z. In fact, I want to represent a 3-axis "cube" and inside it, represent the illuminance in lux in every point in the room, not only in a predefined plane.

How can I do it? Thank you again :D

[Pham Duc Thinh]

Give me something that i can communicate with you :) (For example : Skype Chat)

[Alberto López]

How can I give it that by private?

[Pham Duc Thinh]

I want to talk, not by chat, because i feel tired when i type too much :)

You know, voice is always better than typing.

[Alberto López]

Yes I know :) but now I am not at home, so I cant talk

Answer 2

[Pham Duc Thinh]

:) sorry for late reply :D Do you have simulation code for BER and Channel Delay Spread ?

[Alberto López]

Dont worry, mate! No, I dont have it but actually, I dont need it. I need how to represent three axis (x,y,z) forming a cube, and inside it represent the illuminance. How can i do it?

Answer 3

[Vladimir Fedotov]

If you managed to write this code, can you share it please?

[Toby Bennett]

Does this help?

clear all
clc 
close all
%%%%%%%%%QAM%%%%%%%%%%
M=32;               %   size of the QAM lavel
k=log2(M);          % number of bits per QAM symbol or signal
n=30000,  % number of bits or transmission bits
numSamplespersymbol=1;
datain=randi([0 1], n,1);  % random data generation
 
figure % This will plot first 20 samples of the signal
stem(datain(1:20),'filled');
title ('Random Bits');
xlabel('Bit Index');
ylabel('Binary value');
%%Reshape binary to integer value signal (bi2de)
datainMatrix=reshape(datain, length(datain)/k,k);
datasymbolsIn=bi2de(datainMatrix);
 
%%%plot first 10 symbols %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
figure
stem(datasymbolsIn(1:20),'filled');
title ('RandomSymbols');
xlabel('Symbols Index');
ylabel('Symbol value');
 
%%%QAM Modulation of the signal , these next few lines of code will be on how to modulate
%%%signal using QAMmod command in Matlab
 
dataMod=qammod(datasymbolsIn,M); % Ninary coding, phase offset =0;
 
%%%%%%%%%Wireless channel, add noise to simulate wireless channel, Add channel, AWGN
EbN0=15;   %% energy per bit 
 
snr= EbN0+10*log10(k)-10*log10(numSamplespersymbol);  %% translate into SNR from EbNo
% snr=15;
 
ReceivedSig=awgn(dataMod,snr,'measured'); % For binary
% figure
%%%Plot Transmit symbols and received symbols
 
scatterplot(dataMod,k,0)
title ('Transmit symbols');
hold on
%text(real(dataMod), imag(dataMod), dec2bin(datasymbolsIn))
scatterplot(ReceivedSig,1,0,'g.')  % this command will print for gray code data
% text(real(ReceivedSig), imag(ReceivedSig), dec2bin(datasymbolsIn))
title ('Received symbols after transmission');
 
%%%Demodulation or Data recovery at the Receiver after the channel%%%
dataSymOut=qamdemod(ReceivedSig,M);   % use of demod command to demodulate the QAM data
%convert decimel to binary
dataOutMatrix=de2bi(dataSymOut,k);   %Decimel to bianary
dataOut=dataOutMatrix(:);           %Recoved data after demodulation
 
%%%Compute the BER using biterr command
[numErrors,ber]=biterr(datain,dataOut);
fprintf('\nThe binary coding bit error rate at Physical layer = %5.2e, based on %d errors\n',...
 ber,numErrors')
%%%Compute the BER using biterr command for Gray code 
% [numErrorsG,berG]=biterr(datain,dataOutG);
% fprintf('\nThe Gray coding bit error rate = %5.2e, based on %d errors\n',...
% berG,numErrorsG')
 
%%Eye diagram, vitualisation of filter effects by plot
eyediagram(ReceivedSig(1:2000),numSamplespersymbol*2);