%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% source_gaussian_random.m
%
% This program generates a Gaussian 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 Gaussian Distribution Function
% defined by:
%
%    I(x) = exp(-x^2/(2*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 = 1
% wavelength is in um
%
% Jan 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;

% number of rays required
n  = 100;

% RMS value of the angle in degrees
theta = 5;

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

% MC simulation
phi   = 2*pi*rand(n,1);             % azimuth angle (0,2pi)
Tz    = abs(normrnd(0, theta,n,1)); % zenith angle

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

% to print
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)
