%% Define a beam and visualize it propagating % Define "optical system" using matrices. For this first example, we will just use a free space propagation matrix, and let the beam propagate a distance w0 = 1E-3; % 1mm beam waist lam = 421E-9; % wavelength zR = Zr(w0, lam); % Rayleigh range in m z0 = 0; % location of waist in m syms d; % Define 'd' as a symbolic variable M = propagator(d); [R, w] = q1_inv_func(0, w0, lam, M); % Do all the ABCD and q-parameter math, and return the waist and radius of curvature functions plotResult(w, d, linspace(0,10)) %% Add a lens to the system w0 = 1E-3; % 1mm beam waist lam = 421E-9; % wavelength zR = Zr(w0, lam); % Rayleigh range in m z0 = 0; % location of waist in m syms d; % Define 'd' as a symbolic variable M = matmultiplier(propagator(d), lens(.5), propagator(1)); [R, w] = q1_inv_func(0, w0, lam, M); plotResult(w, d, linspace(0,1,50)) %% Add another to make a beam expander, that enlarges the beam waist by a factor 5 and then collimates it w0 = 1E-3; % 1mm beam waist lam = 421E-9; % wavelength zR = Zr(w0, lam); % Rayleigh range in m z0 = 0; % location of waist in m expansion_factor = 5; syms d1 d2 d3 f1 f2 M = matmultiplier(propagator(d3), lens(f2), propagator(d2), lens(f1), propagator(d1)); [R, w] = q1_inv_func(0, w0, lam, M); % Substitute values into w w = subs(w, {d1, d3}, {1, 0}); % Solve for f1 f1_eq = solve(w - expansion_factor*w0 == 0, f1); f1_eq_subs = subs(f1_eq, d2, 1); % Convert solutions to numeric values f1_val_numeric = double(f1_eq_subs); fprintf('f1 = %.2f m or %.2f m, for a lens separation of 1 meter\n', f1_val_numeric(1), f1_val_numeric(2)) for i = 1:length(f1_val_numeric) f1_val = f1_val_numeric(i); R = subs(R, {d1, d2, d3, f1}, {1, 1, 0, f1_val}); f2_val = solve(1/R == 0, f2); fprintf('f2 = %.2f, for a collimated beam, 5x the original waist, after propagating 1m to the first lens of f1 = %.2f, and propagating another 1m to the second lens\n', double(f2_val), double(f1_val)) end chosen_f1_val = f1_val_numeric(1); fprintf('Chosen f1 = %.2f m \n', chosen_f1_val) M = matmultiplier(propagator(d3), lens(f2_val), propagator(1), lens(chosen_f1_val), propagator(1)); [R, w] = q1_inv_func(0, w0, lam, M); plotResult(w,d3) estimated_expansion_factor = subs(w, d3, 0)/ w0; fprintf('beam is w = %.2f x w0\n', double(estimated_expansion_factor)) beam_size_change = ((subs(w,d3,10) - subs(w,d3,0)) / subs(w,d3,0)) * 100; fprintf('Over 10 m after second lens, beam changes by %.5f percent\n', double(beam_size_change)) %% Function definitions % Input Ray parameter, i.e. height and angle function ret = ray(y, theta) % Parameters % ---------- % y : float or integer in meters % The vertical height of a ray. % theta : float or integer in radians % The angle of divergence of the ray. % % Returns % ------- % mat : 2x1 matrix % [ % [y], % [theta] % ] ret = [y; theta]; end % Ray Transfer Matrix for ideal lens with focal length f function ret = lens(f) % Parameters % ---------- % f : float or integer in meters % Thin lens focal length in meters % % Returns % ------- % mat : 2x2 matrix % [ % [ 1, 0], % [-1/f, 1] % ] ret = [1, 0; -1/f, 1]; end % Ray Transfer Matrix for propagation over a distance d function ret = propagator(d) % Parameters % ---------- % d : float or integer % Distance light is propagating along the z-axis. % % Returns % ------- % mat: 2x2 matrix % [ % [1, d], % [0, 1] % ] ret = [1, d; 0, 1]; end % Multiplying the matrices together. mat1 is the last matrix the light interacts with. function ret = matmultiplier(mat, varargin) % Parameters % ---------- % mat1 : 2x2 ABCD matrix % Last matrix light interacts with. % varargin : 2x2 ABCD matrices % From left to right, the matrices should be entered such that the leftmost matrix interacts % with light temporally after the rightmost matrix. % % Returns % ------- % Mat : 2x2 matrix % The ABCD matrix describing the whole optical system. ret = mat; for i = 1:length(varargin) ret = ret * varargin{i}; end end % Adding Gaussian beam parameters function zr = Zr(wo, lam) % Parameters % ---------- % wo : float or integer % Beam waist radius in meters. % lam : float or integer % Wavelength of light in meters. % % Returns % ------- % zr : float or integer % Rayleigh range for given beam waist and wavelength. zr = pi * wo^2 / lam; end % Calculate beam waist radius function w0 = W0(zr, lam) % Parameters % ---------- % zr : float or integer % Rayleigh range in meters. % lam : float or integer % Wavelength of light in meters. % % Returns % ------- % w0 : float or integer % Beam waist radius in meters. w0 = sqrt(lam * zr / pi); end function [z_out, zr_out] = q1_func(z, w0, lam, mat) % Parameters % ---------- % z : float or integer % Position of the beam waist in meters. % w0 : float or integer % Radial waist size in meters (of the embedded Gaussian, i.e. W0/M). % lam : float or integer % Wavelength of light in meters. % mat : 2x2 matrix % The ABCD matrix describing the optical system. % % Returns % ------- % z_out : float or integer % Position of the beam waist after the optical system. % zr_out : float or integer % Rayleigh range of the beam after the optical system. A = mat(1, 1); B = mat(1, 2); C = mat(2, 1); D = mat(2, 2); % Calculate Rayleigh range for the given beam waist and wavelength zr = Zr(w0, lam); % Calculate real and imaginary parts real_part = (A * C * (z^2 + zr^2) + z * (A * D + B * C) + B * D) / (C^2 * (z^2 + zr^2) + 2 * C * D * z + D^2); imag_part = (zr * (A * D - B * C)) / (C^2 * (z^2 + zr^2) + 2 * C * D * z + D^2); % Output the new position and Rayleigh range z_out = real_part; zr_out = imag_part; end function [R, w] = q1_inv_func(z, w0, lam, mat) % Parameters % ---------- % z : float or integer % Position of the beam waist in meters. % w0 : float or integer % Radial waist size in meters (of the embedded Gaussian, i.e. W0/M). % lam : float or integer % Wavelength of light in meters. % mat : 2x2 matrix % The ABCD matrix describing the optical system. % % Returns % ------- % R : float or integer % Radius of curvature of the wavefront in meters. % w : float or integer % Radius of the beam in meters. A = mat(1, 1); B = mat(1, 2); C = mat(2, 1); D = mat(2, 2); % Calculate Rayleigh range for the given beam waist and wavelength zr = Zr(w0, lam); % Calculate real and imaginary parts real_part = (A * C * (z^2 + zr^2) + z * (A * D + B * C) + B * D) / (A^2 * (z^2 + zr^2) + 2 * A * B * z + B^2); imag_part = -zr * ((A * D) - (B * C)) / (A^2 * (z^2 + zr^2) + (2 * A * B * z) + B^2); % Calculate radius of curvature and beam radius R = 1/real_part; w = sqrt(-lam / imag_part / pi); end function plotResult(func, var, rang) % Parameters % ---------- % func : symfun (symbolic function of one variable) % Symbolic function defining the beam width after the last optical element. % var : symbolic variable % Variable in func that will be plotted. % rang : array (optional) % Array of values along the optical axis to be plotted. Default range is 0 to 3 with step size 0.01. % Set default range if not provided if nargin < 3 rang = 0:0.01:3; end % Create a numeric function from the symbolic function func_handle = matlabFunction(func, 'Vars', var); % Compute the values of the function and its negative y_vals = func_handle(rang) * 1E3; % in mm y_neg_vals = -y_vals; % Create a figure figure; set(gcf, 'Position', [100 100 950 750]); % Plot the function and its negative plot(rang, y_vals, 'b'); hold on; plot(rang, y_neg_vals, 'b'); % Fill the area between the curves with translucent blue color fill([rang, fliplr(rang)], [y_vals, fliplr(y_neg_vals)], 'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none'); % Add grid and labels grid on; xlabel('Optical Axis (m)', 'FontSize', 16); ylabel('Beam size (mm)', 'FontSize', 16); end