%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% achromate_ynu.m
% Relatively advanced achromatic lens design
% using ynu ray tracing. Radius of curvatures
% are found by Monte Carlo Method
%
% Inputs:
%         f  Focal length of system
%     index  Refractive indices (FdC) of glasses
%     t1,t2  center thickness of glasses
% Outputs:
%         R  Radii of curvatures
%
% Nov 2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;
% system focal length
fsystem = 300;
R1= fsystem/2;
R2= inf;
R3= inf;
t1 = 4;
t2 = 3;
index = [1.5224 1.5168 1.5143; ... % first  lens
         1.6875 1.6727 1.6667];    % second lens

mbest = 1e9; % best value of the merit function
N = 1e7;     % number of MC trials

R2 = -rand(N,1)*1000; % select second and
R3 = -rand(N,1)*1000; % third surface randomly

% *** MC loop ***
for i = 1:N
  R = [inf R1 R2(i) R3(i) inf];
  t = [1 t1 t2];

    %------- for each wavelength
    for j=1:3
      y(1) = 1;
      u(1) = 0;
      n = [1  index(1,j) index(2,j) 1];
      % apply y-nu ray tracing method
      for  k = 1:3
        y(k+1) = y(k) + u(k)*t(k);
        u(k+1) = (n(k)*u(k) - (n(k+1)-n(k))*y(k+1)/R(k+1))/n(k+1);
      end
      EFFL(j) = -y(1)/u(4);
    end
    %-------

    % define merit function to be minimized and weights
    w1 = 0.5;
    w2 = 1.0;
    meritf = (w1*(EFFL(2)-fsystem)^2 + w2*(EFFL(1)-EFFL(3))^2)/(w1+w2);
    if  meritf < mbest
      mbest = meritf;
      fprintf('MF=%f R2=%9.3f R3=%9.3f => f(b)=%9.3f | f(y)=%9.3f | f(r)=%9.3f\n',...
              mbest,R2(i),R3(i),EFFL(1),EFFL(2),EFFL(3));
      Radius = [R1 R2(i) R2(i) R3(i)];
    end
end
% *** End of MC loop ***

fprintf('Target Focal length for d-line = %9.3f\n')
fprintf('Center Thickness of Lenses:\n')
fprintf('%9.3f\n',t1)
fprintf('%9.3f\n',t2)
fprintf('Radius of Curvatures:\n')
for i=1:length(Radius)
    fprintf('%9.3f\n',Radius(i))
end

%//////////////////////////////////////////////
% Example output
%
% Target Focal length for d-line = 300.000
% Center Thickness of Lenses:
%    4.000
%    3.000
% Radius of Curvatures:
%   150.000
%  -154.467
%  -154.467
%  -747.311
%
%/////////////////////////////////////////////
