Complex Number Ignition Error: Ray transfer matrix analysis

I'm trying to simulate beam propagation with Ray transfer matrix analysis.

But, the result graphs are not that I expected.(Graph of Beam radius is just flat)

I hope someone could find my mistake.

clear,clc

clear;

clc;

% Parameters

lambda = 632.8e-6; % Wavelength in millimeters (632.8 nm)

w0 = 0.24; % Initial beam waist (0.24 mm)

zR = pi * w0^2 / lambda; % Rayleigh range in mm

z1 = 10; % Position of the first lens (10 mm)

z2 = 250; % Position of the second lens (250 mm)

f1 = -25; % Focal length of the first lens (-25 mm)

f2 = 250; % Focal length of the second lens (250 mm)

% Initial complex beam parameter q0

q0 = vpa(1i * zR);

% Discretize the propagation distance for plotting

z = linspace(0, 500, 1000); % Range from 0 to 500 mm

% Initialize arrays for q, beam radius, and curvature radius

q = zeros(size(z));

wz = zeros(size(z));

Rz = zeros(size(z));

% Calculate q, beam radius, and curvature radius

for i = 1:length(z)

if z(i) <= z1

% Propagation from 0 to 10 mm

ABCD = [1, z(i); 0, 1];

elseif z(i) <= z2

% Propagation from 10 mm to 250 mm through first lens

ABCD = [1, z1; 0, 1] * [1, 0; -1/f1, 1] * [1, z(i) - z1; 0, 1];

else

% Propagation from 250 mm to 500 mm through second lens

ABCD = [1, z1; 0, 1] * [1, 0; -1/f1, 1] * [1, z2 - z1; 0, 1] * [1, 0; -1/f2, 1] * [1, z(i) - z2; 0, 1];

end

q(i) = (ABCD(1,1) * q0 + ABCD(1,2)) / (ABCD(2,1) * q0 + ABCD(2,2));

wz(i) = sqrt(lambda / pi * imag(1/q(i)));

Rz(i) = 1 / real(1/q(i));

end

% Plotting the beam radius along the propagation

figure;

plot(z, wz); % Already in mm

xlabel('Propagation distance z (mm)');

ylabel('Beam radius w(z) (mm)');

title('Gaussian Beam Propagation Through Two Lenses');

grid on;

% Plotting the curvature radius along the propagation

figure;

plot(z, Rz); % Already in mm

xlabel('Propagation distance z (mm)');

ylabel('Curvature radius R(z) (mm)');

title('Curvature Radius of Gaussian Beam Through Two Lenses');

grid on;

% Display final parameters at specific positions

fprintf('Beam waist at 500 mm: %.3f mm\n', wz(end));

fprintf('Curvature radius at 500 mm: %.3f mm\n', Rz(end));

Answer 1

[Torsten]

imag(1./q)*lambda/pi

is negative, thus

sqrt(imag(1./q)*lambda/pi)

is complex.

MATLAB does not know how you want it to plot complex numbers (and I don't know either).