How to draw multiple surface plots on the same mesh?

I use this code to draw a Gaussian normal distribution. I want to plot several distributions at equal distance from the original one (like the figure).

How can I do so?

lambda = 1;               % wavelength

E0 = 1;                   % max value of E-field amplitude

% ==================================================================
% prepare the coordinates mesh
% ----------------------------

zLin = linspace(-201,201,400) * lambda;     % z-coordinate (horizontal)

rLin = linspace(-65,65,100) * lambda;       % radius-coordinate (vertical)

[z,r] = meshgrid(zLin,rLin);                % create mesh


% ===================================================================

wFactor = 4;                 
    
w0 = wFactor * lambda;                  % minimal waist

zR = pi* w0^2 / lambda;                 % Rayleigh length

w = w0 * sqrt( 1 + (z / zR) .^2 );      % Beam radius at any arbitrary distance (equ. 3)

% Irradiance
eta = 377;                              % for free space
I0 = abs(E0)^2/eta;

I = I0 .* (w0 ./ w).^2 .* exp( -2*(r./w).^2 );

surf(z,r,I); shading interp

[VBBV]

Use contourf function as

%if true
  % code
% end
v2 = 0:0.00009:0.0025;
contourf(z,r,I,v2,'edgecolor','none');shading interp;

Answer 1

[Durganshu]

Did you try hold on?

figure;
surf(z1,r1,I1);
hold on
zlin2= linspace(-101,101,300) * lambda; %whatever shifted values you want
rLin2 = linspace(-65,65,100) * lambda; %whatever shifted values you want
[z2,r2] = meshgrid(zLin2,rLin2);
% similarly define I2 and then plot
surf(z2,r2,I2);

[Muhamed Sewidan]

Thank you @Durganshu for your fast response.

In fact, I tried this approach but the two distributions are overlapped on each other. How can I add the shift?

[Durganshu]

Can you share the snippet of code that you're using for shifting?

[Muhamed Sewidan]

hold on

zlin2= linspace(-101,101,300) * lambda; %whatever shifted values you want

rLin2 = linspace(-65,65,100) * lambda; %whatever shifted values you want

[z2,r2] = meshgrid(zlin2,rLin2);

w2 = w0 * sqrt( 1 + (z2 / zR) .^2 );

% similarly define I2 and then plot
I2 = I0 .* (w0 ./ w2).^2 .* exp( -2*(r2./w2).^2 );

surf(z2,r2,I2);

[Durganshu]

Sorry for a bit delay.

The code is working fine. Please elaborate on what you mean by shift? If you'll plot the first surface within -201 to 201 and the another one from -101 to 101, then definitely it'll be plotted inside -201 to 201. If you want it to be plotted outside the range of the first plot then make sure its range shouldn't be the subset of the first range.

For example, linspace (201, 301) will not overlap with the plot from (-201, 201).

[Muhamed Sewidan]

Never mind, thank you for your interest.

Should I change rlin2 also the same way I change zlin2?

I did as you said, but it didn't plot I2, it is just extended the mesh of the first plot. Please have a look at my code after adding your suggestion.

lambda = 1;               % wavelength

E0 = 1;                   % max value of E-field amplitude

% ==================================================================
% prepare the coordinates mesh
% ----------------------------

zLin = linspace(-201,201,400) * lambda;     % z-coordinate (horizontal)

rLin = linspace(-65,65,100) * lambda;       % radius-coordinate (vertical)

[z,r] = meshgrid(zLin,rLin);                % create mesh

% ===================================================================
% Applying the code:
% ------------------

wFactor = 4;                 
    
w0 = wFactor * lambda;                  % minimal waist

zR = pi* w0^2 / lambda;                 % Rayleigh length

w = w0 * sqrt( 1 + (z / zR) .^2 );      % Beam radius at any arbitrary distance (equ. 3)

% Irradiance
eta = 377;                              % for free space
I0 = abs(E0)^2/eta;
I = I0 .* (w0 ./ w).^2 .* exp( -2*(r./w).^2 );
surf(z,r,I)
shading interp

%% Second source:

hold on

zlin2= linspace(201,602,400) * lambda; %whatever shifted values you want

rLin2 = linspace(-65,65,100) * lambda; %whatever shifted values you want

[z2,r2] = meshgrid(zlin2,rLin2);

w2 = w0 * sqrt( 1 + (z2 / zR) .^2 );

% similarly define I2 and then plot
I2 = I0 .* (w0 ./ w2).^2 .* exp( -2*(r2./w2).^2 );

surf(z2,r2,I2);

[Durganshu]

Again, your code is absolutely okay. Let's see the differences between I and I2:

1. I includes w and I2 includes w2.
2. I includes r and I2 includes r2.

Now, w2 and w depend on z2 and z. Looking at the values of z2 and z (that you get in the MATLAB workspace), it is clearly visible that there must be a continuos boundary between the curves I and I2. That is also visible in the curve.

Now, r2 and r are the same.

COCLUSION: I2 is just a natural extension of I. I2 starts from where the range of I ends. Basically, by writing this code, you have extended the range of I to 602. It could be written that range(I') = range(I) + range(I2). That's why you're getting that kind of plot.

My answers were to guide you and not give you the final answer. You never specified anything about the inputs and what exactly is going in the problem. The code is alright. You just need to be clear about what exactly curves do you want. If you want a different curve, then you need to put something else in the equation for I2.

How can I tell you if you need to change r2 or not? It depends on your original problem statement.

So, please check on the exact requirements and change the equation accordingly.

If you're still stuck somewhere, then post what kind of output/shifts you want with detailed explanation.

[Muhamed Sewidan]

Thank you for your valuable assist.

I wanted basically to draw another gaussian normal distribution next to the original one. So, should I add some offset value to I2 equation?

[Durganshu]

Yes, that should work. Did you try it? All you're supposed to do is tweaking the I2 function in order to find the appropriate fit.

[Muhamed Sewidan]

Thank you I will try it out.