How to plot multiple curves in one figure for changable values of phi and psi?

I have phi & psi and each one has values from 0 to 30 as shown below:

What I am trying to do is plotting the received optical based on the above values.

For example, let's consider the first case below: phi = 0 & psi = 0,5,10,15,20,25,30

phi = 0  
psi = 0
      5
      10
      15
      20
      25
      30 

I drew them like this below:

The above figure was only for one case.

What I want now is to draw the received optical power for the values as shown below and put them all in one figure:

Let's name them as

case_1: phi = 0 & psi = 0,5,10,15,20,25,30

case_2: phi = 5 & psi = 0,5,10,15,20,25,30

case_3: phi = 10 & psi = 0,5,10,15,20,25,30

case_4: phi = 15 & psi = 0,5,10,15,20,25,30

case_5: phi = 20 & psi = 0,5,10,15,20,25,30

case_6: phi = 25 & psi = 0,5,10,15,20,25,30

case_7: phi = 30 & psi = 0,5,10,15,20,25,30

I need all above cases in one figure, on the same x-axis and y-axis.

I tried to draw them with this below piece of code:,

figure
hold on;
for i = 1 : 7
    for j = 1 : 7        
        plot(X_r,P_rec_dbm(:,21,1,i,j)); 
    end
end
xlabel('X (m)');
ylabel('Received power (dBm)');
grid on;

but it only shows me one curve.

Any assistance please?

Below my code

===========

close all;
clear variables;
clc;

%% Simulation Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    Main Simulation Parameters                     %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------%
% NUMBER OF LIGHT SOURCES %
%-------------------------%
N_t = 1; % Number of light sources
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                          AP Parameters                            %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------%
% LIGHT SOURCES GEOMETRY %
%----------%
L=20; W=20; H=3; % Length, width and height of the room (m)
theta_half = 30;
m = -log(2)./log(cosd(theta_half)); % Lambertian order of emission
coord_t = [0 0 0]; % Positions of the light sources
n_t_LED = [0, 0, -1]; n_t_LED = n_t_LED/norm(n_t_LED); % Normalized normal vector of each light source
n_t = repmat(n_t_LED, N_t, 1); % Normalized normal vectors of the light sources

%-------------------------------------%
% LIGHT SOURCES ELECTRICAL PARAMETERS %
%-------------------------------------%
P_LED = 2.84; % Average electrical power consumed by each light source (W)

param_t = {coord_t, n_t, P_LED, m};
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                          Rx Parameters                            %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------%
% PHOTODETECTOR PARAMETERS %
%--------------------------%
A_det = 1e-4;        % Photoreceiver sensitive area (m²)
FOV=60*pi/180;     % Field-of-view of the photoreceiver
T_s = 1;                      % Gain of the optical filter (ignore if not used)
index=1.5;                  % Refractive index of the Rx concentrator/lens (ignore if not used)
G_Con = (index^2)/sin(FOV);   % Gain of an optical concentrator; ignore if no lens is used
n_r = [0, 0, 1]; % Normal vector of the photoreceiver
n_r = n_r/norm(n_r); % Normal vector of the photoreceiver (normalized)

%---------------------------%
% RECEIVER PLANE PARAMETERS %
%---------------------------%
step = 0.5; % Distance between each receiving point (m)
X_r = -L/2:step:L/2; % Range of Rx points along x axis
Y_r = -W/2:step:W/2; % Range of Rx points along y axis
N_rx = length(X_r); N_ry = length(Y_r); % Number of reception points simulated along the x and y axis
z_ref = 0.85; % Height of the receiver plane from the ground (m)
z = z_ref-H; %z=-1.65; % Height of the Rx points ("-" because coordinates system origin at the center of the ceiling)
if ( abs(z) > H )
    fprintf( 'ERROR: The receiver plane is out of the room.\n' );
    return
end
param_r = {A_det, n_r, FOV}; % Vector of the Rx parameters used for channel simulation

%% LOS received optical power calculation

phi = 0:5:30;
psi = (0:5:30).';
N_phi = numel(phi);
N_psi = numel(psi);
H0_LOS = zeros(N_rx,N_ry,N_t,N_phi,N_psi);
P_rec_dbm = zeros(N_rx,N_ry,N_t,N_phi,N_psi);
T = param_t{1}(1,:);
P_t = param_t{3};

for r_x = 1:N_rx
    for r_y = 1:N_ry
        for i_t = 1:N_t
            x = X_r(r_x); y = Y_r(r_y);
            R = [x,y,z];
            v_tr = (R-T)./norm(R-T);
            d_tr = sqrt(dot(R-T,R-T));          
            H0_LOS(r_x,r_y,i_t,:,:) = (m+1)*A_det/(2*pi*d_tr^2)*cosd(phi).^m*cosd(psi); 
        end
    end
end

P_r_LOS = P_t.*H0_LOS.*T_s.*G_Con;
P_rec_dbm = 10*log10(P_r_LOS*1000);

%% Plotting

figure
[x,y] = meshgrid(X_r,Y_r);
P_rec_dbm_extracted = P_rec_dbm(:,:,:,3,5);
size(x);
size(y);
size(P_rec_dbm_extracted);
meshc(X_r,Y_r,P_rec_dbm_extracted);
colorbar
xlabel( 'X(m)' );
ylabel( 'Y(m)' );
zlabel( 'Received power (dBm)' );
axis([-L/2 L/2 -W/2 W/2 min(min(P_rec_dbm_extracted)) max(max(P_rec_dbm_extracted))]);


figure
plot(X_r,P_rec_dbm_extracted)
xlabel( 'X(m)' );
ylabel( 'Received power (dBm)' );

figure
hold on;
for i = 1 : 7
    for j = 1 : 7        
        plot(X_r,P_rec_dbm(:,21,1,i,j)); 
    end
end
xlabel('X (m)');
ylabel('Received power (dBm)');
grid on;

[dpb]

figure
hold on;
for i = 1 : 7
    for j = 1 : 7        
        plot(X_r,P_rec_dbm(:,21,1,i,j)); 

would plot 49 lines on the one axes; if they don't show up as separate lines it would imply the data aren't different for the varioul i,j combinations in array P_rec_dbm

Use the debugger and see; that's a pretty complex storage arrangement with a 5D array, I'd guess something went wrong in building it.

[Haitham AL Satai]

@dpb I got your point dear. Thank you very much for your comment.

[Torsten]

Why do you modify the code at the wrong places ?

H0_LOS(r_x,r_y,i_t,:,:) = (m+1)*A_det/(2*pi*d_tr^2)*cosd(phi).^m.*cosd(psi);

instead of

H0_LOS(r_x,r_y,i_t,:,:) = (m+1)*A_det/(2*pi*d_tr^2)*cosd(phi).^m*cosd(psi);

%% Simulation Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    Main Simulation Parameters                     %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------%
% NUMBER OF LIGHT SOURCES %
%-------------------------%
N_t = 1; % Number of light sources
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                          AP Parameters                            %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------%
% LIGHT SOURCES GEOMETRY %
%----------%
L=20; W=20; H=3; % Length, width and height of the room (m)
theta_half = 30;
m = -log(2)./log(cosd(theta_half)); % Lambertian order of emission
coord_t = [0 0 0]; % Positions of the light sources
n_t_LED = [0, 0, -1]; n_t_LED = n_t_LED/norm(n_t_LED); % Normalized normal vector of each light source
n_t = repmat(n_t_LED, N_t, 1); % Normalized normal vectors of the light sources
%-------------------------------------%
% LIGHT SOURCES ELECTRICAL PARAMETERS %
%-------------------------------------%
P_LED = 2.84; % Average electrical power consumed by each light source (W)
param_t = {coord_t, n_t, P_LED, m};
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                          Rx Parameters                            %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------%
% PHOTODETECTOR PARAMETERS %
%--------------------------%
A_det = 1e-4;        % Photoreceiver sensitive area (m²)
FOV=60*pi/180;     % Field-of-view of the photoreceiver
T_s = 1;                      % Gain of the optical filter (ignore if not used)
index=1.5;                  % Refractive index of the Rx concentrator/lens (ignore if not used)
G_Con = (index^2)/sin(FOV);   % Gain of an optical concentrator; ignore if no lens is used
n_r = [0, 0, 1]; % Normal vector of the photoreceiver
n_r = n_r/norm(n_r); % Normal vector of the photoreceiver (normalized)
%---------------------------%
% RECEIVER PLANE PARAMETERS %
%---------------------------%
step = 0.5; % Distance between each receiving point (m)
X_r = -L/2:step:L/2; % Range of Rx points along x axis
Y_r = -W/2:step:W/2; % Range of Rx points along y axis
N_rx = length(X_r); N_ry = length(Y_r); % Number of reception points simulated along the x and y axis
z_ref = 0.85; % Height of the receiver plane from the ground (m)
z = z_ref-H; %z=-1.65; % Height of the Rx points ("-" because coordinates system origin at the center of the ceiling)
if ( abs(z) > H )
    fprintf( 'ERROR: The receiver plane is out of the room.\n' );
    return
end
param_r = {A_det, n_r, FOV}; % Vector of the Rx parameters used for channel simulation
%% LOS received optical power calculation
phi = 0:5:30;
psi = (0:5:30).';
N_phi = numel(phi);
N_psi = numel(psi);
H0_LOS = zeros(N_rx,N_ry,N_t,N_phi,N_psi);
P_rec_dbm = zeros(N_rx,N_ry,N_t,N_phi,N_psi);
T = param_t{1}(1,:);
P_t = param_t{3};
for r_x = 1:N_rx
    for r_y = 1:N_ry
        for i_t = 1:N_t
            x = X_r(r_x); y = Y_r(r_y);
            R = [x,y,z];
            v_tr = (R-T)./norm(R-T);
            d_tr = sqrt(dot(R-T,R-T));          
            H0_LOS(r_x,r_y,i_t,:,:) = (m+1)*A_det/(2*pi*d_tr^2)*cosd(phi).^m.*cosd(psi); 
        end
    end
end
P_r_LOS = P_t.*H0_LOS.*T_s.*G_Con;
P_rec_dbm = 10*log10(P_r_LOS*1000);
%% Plotting
figure
[x,y] = meshgrid(X_r,Y_r);
P_rec_dbm_extracted = P_rec_dbm(:,:,:,3,5);
size(x);
size(y);
size(P_rec_dbm_extracted);
meshc(X_r,Y_r,P_rec_dbm_extracted);
colorbar
xlabel( 'X(m)' );
ylabel( 'Y(m)' );
zlabel( 'Received power (dBm)' );
axis([-L/2 L/2 -W/2 W/2 min(min(P_rec_dbm_extracted)) max(max(P_rec_dbm_extracted))]);

figure
plot(X_r,P_rec_dbm_extracted)
xlabel( 'X(m)' );
ylabel( 'Received power (dBm)' );

figure
hold on;
for i = 1 : N_phi
    for j = 1 : N_psi        
        plot(X_r,P_rec_dbm(:,21,1,i,j));
        legendInfo{(i-1)*N_psi+j} = ['\Phi =',num2str((i-1)*(phi(2)-phi(1))),' ,\Psi =',num2str((j-1)*(psi(2)-psi(1)))];
        legend(legendInfo);
    end
end
xlabel('X (m)');
ylabel('Received power (dBm)');
grid on;

[Haitham AL Satai]

@Torsten Oops. I think I mistook when I rewrite it. Thanks a lot dear sir. Is it possible to add the legends on the results like below, but for the above picture? I mean mentioning the status of every phi & psi. For example:

________ phi = 0 & psi = 0

________ phi = 0 & psi = 5

________ phi = 0 & psi = 10

....

________ phi = 5 & psi = 0

________ phi = 5 & psi = 5

...

...

________ phi = 30 & psi = 30

[Torsten]

See above.

[Haitham AL Satai]

@Torsten Thanks a lot dear sir. May you explain the idea of this line, if you do not mind?

legendInfo{(i-1)*N_psi+j} = [ '\Phi =' ,num2str((i-1)*(phi(2)-phi(1))), ',\Psi =' ,num2str( (j-1)*(psi(2)-psi(1)))];

What is the idea in this part?

legendInfo{(i-1)*N_psi+j}

[Torsten]

It defines the legend that should be printed if phi = phi(i) and psi = psi(j).

[Haitham AL Satai]

@Torsten It means store what is resulted from here

[ '\Phi =' ,num2str((i-1)*(phi(2)-phi(1))), ',\Psi =' ,num2str( (j-1)*(psi(2)-psi(1)))];

and store it here

legendInfo

I understood the above ideas, just this one (i-1)*N_psi+j

[Torsten]

What do you get for (i-1)*N_psi + j for i = 1,...,N_phi and j = 1,...,N_psi ?

1,2,3,...,N_psi, N_psi+1,Npsi+2,...,2*N_psi, ... ,(Nphi-1)*N_psi+1,(N_phi-1)*N_psi+2,...,N_phi*N_psi.

So it transforms the 2d-array (i,j) (1<=i<=N_phi, 1<=j<=N_psi) to an 1d-array in the above order.

[Haitham AL Satai]

@Torsten Thanks a lot for your patience, cooperation, and your assistance.

[Torsten]

You're welcome.