%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% freeform_circular.m
%
% Simple freeform lens design for uniform circular illumination.
%
% Source is assumed to be point-like and Lambertian.
% Distance unit is mm.
%
% May 2025
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;
% Inputs
d = 5000;   % distance between source and target
n = 1.4935; % typical index of refraction of PMMA
A = 90/2;   % half emitting angle of lens in degrees

% initial calculations
R = d*tand(A); % radius of target
Ri= 4;         % inner radius of first surface
z0= 8;         % outer radius of second surface
N = 50;        % number of points

t = zeros(N,1); % angle from the source
r = zeros(N,1); % radius on target

dFlux = Integrate(0,pi/2)/N; % flux in each segment

% compute angle and radii
% Angle theta (t) is only for lambertian.
for i = 2:N
  t(i) = asin( sqrt(2*dFlux+sin(t(i-1))^2) );
  r(i) = R*sqrt(i/N);
end

% points on inner surface of the lens (circle)
Fx = Ri*sin(t);
Fz = Ri*cos(t);
% points on target (hit positions on target)
Qx = r;
Qz = d*ones(N,1);
% points on outer surface of the lens
Px = zeros(N,1);
Pz = z0*ones(N,1);
% unit vectors
L0 = [0 0 0];
nv = [0 0 1];
P  = [0 0 z0];

for i=2:N
  % unit vector in the direction of incident ray
  J  = [sin(t(i)) 0 cos(t(i))];

  % hit position on the target
  Q  = [r(i) 0 d];

  % intersection of line and plane
  s  = dot(P-L0,nv) / dot(J,nv);

  % line equation
  P  = L0 + s*J;

  % unit vector in the direction of outgoing ray
  O  = (Q-P)/norm(Q-P);

  % normal (unit) vector to the surface
  nv = (O - n*J) / sqrt(1+n*n-2*n*dot(O,J));

  % location on the freeform lens surface
  Px(i) = P(1);
  Pz(i) = P(3);
end

% draw lens profile
plot(+Fx,Fz,'b.')
hold on
plot(-Fx,Fz,'b.')
plot(+Px,Pz,'b.');
plot(-Px,Pz,'b.');
scale = max(Px)*1.1;
axis([-scale scale 0 2*scale])

title('Freeform Lens Profile')
xlabel('Distance (mm)')
ylabel('Distance (mm)')
grid on
hold off

disp('ALL OK.')


%************************************%
%         Helper functions           %
%                                    %
%             /\_/\                  %
%            ( o.o )                 %
%             > ^ <                  %
%************************************%

%----------------------------------------------------
% light intensity distribution of the LED source, Lambertian
function y = I(t)
   y = 1 * cos(t);
end
%----------------------------------------------------
% integrand
function y = f(t)
   y = I(t) * sin(t);
end
%----------------------------------------------------
% Integral of f() above over the range a,b.
function etf = Integrate(a,b)
  n = 2000;
  h = (b-a)/n;
  etf = 0;
  for i=1:n-1
    etf = etf + f(a+i*h);
  end
  etf = etf * h;
end

