%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% source_lambertian_random.m
%
% This program generates a Lambertian light source file to be used in
% NonSequeantial mode of Zemax OpticStudio. The source is assumed to be
% point-like located at (0,0,0). The ray directions (angles) are selected 
% randomly from the non-uniform Lambertian Distribution Function 
% defined by:
%
%    I(theta) = I0 * cos(theta)
%
% Output file format:
% ------------------------------------------
% number_of_rays dimension_flag
% x1 y1 z1 l1 m1 n1 intensity1 wavelength1
% x2 y2 z2 l2 m2 n2 intensity2 wavelength2
% ...
% ------------------------------------------
% dimension_flag is 0 for meters,
%                   1 for inches,
%                   2 for centimeters,
%                   3 for feet, and
%                   4 for millimeters
% xyz        are initial coordinates of the ray
% lmn        are direction cosines in 3D
% intensity  is 1.0 [you can change as desired weights between (0,1)]
% wavelength is in um
%
% Jan 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;

% total number of rays requared
n  = 100;

% output file name
filename = 'lambertian_rnd.dat';

% MC simulation of angles
x     = rand(n,1);       % set of random number (0,1)
phi   = 2*pi*rand(n,1);  % azimuthal angle (0,2pi)
Tz    = asind(sqrt(x));  % zenith angle (0,pi/2)

% direction cosines
cosTx = sind(Tz).*cos(phi);
cosTy = sind(Tz).*sin(phi);
cosTz = cosd(Tz).*ones(n,1);
waveLength = 0.55*ones(n,1); % um

% to print file
fileID = fopen(filename,'w');
fprintf(fileID,'%d 4\n',n);

for i = 1:n
  fprintf(fileID,'%d %d %d %8.5f %8.5f %8.5f %d %6.4f\n', ...
                  0,0,0, cosTx(i),cosTy(i),cosTz(i), 1, waveLength(i));
end

fclose(fileID);
fprintf('Data written to file: %s\n',filename)

