% thin_lens_matrix_root
%
% Trace rays through thin lens using ray matrices, and then find the 
%   focus via root-finder.

global gF

% parameters
f = 50; % focal length of thin lens in whatever units you like (let's say mm)

% we need to pass a parameters to our function via global variables;
%   prepend global names with a 'g' and keep them separate for clarity
gF = f;


% function to return ray height at the end of the imaging system
%   inputs:
%     d -> distance of propagation after the lens
%   outputs:
%     yout -> height of ray after propagation through lens and distance d
%
function yout = rayfunc(d)

  % get and rename global variables
  global gF
  f = gF;

  thin_lens_M = [1, 0; -1/f, 1];
  exterior_free_space_M = [1, d; 0, 1];

  % note order of multiplication!
  %   use '...' to continue an expression on the next line
  composite_M = ...
    exterior_free_space_M * ...
    thin_lens_M;

  % now trace an initially parallel ray, height = 1
  yin_vec = [1; 0];
  yout_vec = composite_M * yin_vec;
    
  % return the height of the final ray
  yout = yout_vec(1);
  
end %function rayfunc


% find distance where ray crosses the axis, using 'fzero' root finder
%   i.e., we want to find the distance d such that rayfunc(d) = 0
%   Note that we need to supply an initial guess, say 10f here, just
%   to make it not so trivial.
d_focus = fzero('rayfunc', 10*f);

% print result
printf('numerically computed focal length = %f mm\n', d_focus);

% compare to what you know it should be
printf('analytically computed focal length = %f mm\n', f);