% calculates the luminous intensity on a surface caused by some light
% emitters

% all units in mm

%% definitions
% Board dimensions
board_x = 100;
board_y = 160;

% distance bpard to emitters
distance = 100;

% emitter positions
% one emitter per row, each with x and y coordinates
emitters = [-20  10; ...
            -20 150; ...
            120 150; ...
            120  10];

% number of simulated board elements in each dimension
elements_x = board_x/1;
elements_y = board_y/1;

% emitter profile
% EDISON 3 Watt Emitter EDEV-SLC1-03
emitter_angle       = [   0,    5,   10,   15,   20,   25,   30,   35,   40,   45,   50,   55,   60,   65,   70,   75,   80,   85,   90];
emitter_intensity_x = [1.00, 0.99, 0.98, 0.97, 0.96, 0.95, 0.90, 0.85, 0.80, 0.75, 0.60, 0.40, 0.35, 0.20, 0.12, 0.06, 0.05, 0.04, 0.03];
emitter_intensity_y = [1.00, 1.00, 1.00, 0.98, 0.95, 0.93, 0.90, 0.86, 0.83, 0.79, 0.75, 0.72, 0.68, 0.63, 0.56, 0.63, 0.38, 0.27, 0.20];

%% simulation data
step_x = board_x/elements_x;
step_y = board_y/elements_y;
[X,Y]  = meshgrid(linspace(step_x/2, board_x-step_x/2, elements_x), linspace(step_y/2, board_y-step_y/2 ,elements_y));
intensity = zeros(size(X));


%% do the simulation for each emitter
for n_emitter = 1 : length(emitters)
    % radial distance from emitter
    deltaX = abs(X - emitters(n_emitter, 1));
    deltaY = abs(Y - emitters(n_emitter, 2));
    
    % intensity due to distance
    local_intensity = 1 ./ (deltaX.^2 + deltaY.^2 + distance^2);
    
    % angle from emitter to surface
    angle_x = atand(deltaX/distance); 
    angle_y = atand(deltaY/distance);
    
    % intensity due to angle
    local_intensity = local_intensity .* interp1(emitter_angle, emitter_intensity_x, angle_x) .* interp1(emitter_angle, emitter_intensity_y, angle_y);
    
    % add to gloabl intensity
    intensity = intensity + local_intensity;
end


%% output data

fprintf('Maximum Intensity: %f (100.0%%)\n', max(intensity(:)));
fprintf('Minimum Intensity: %f (%5.1f%%)\n', min(intensity(:)), min(intensity(:))/max(intensity(:)));

% normalize
intensity = intensity / max(intensity(:));

newplot 
surf(X, Y, intensity, 'EdgeColor','none');
% imagesc(intensity);