% Berechnung der Relativen Lichtintensität eines Planarstrahlers auf eine 
% Parallele Platte mit Variablem Abstand
% Geschrieben am 20.04.2011 von Marc Seiffert
% Nutzung, Kommerziell sowie Privat gestattet
% Veröffentlichen von geänderten Versionen nur mit ausdrücklichem Hinweis
% auf Veränderung !
%
%
% Version 1.0



clc
clear all

% Alle Einheiten in mm

Feldbreite              = 300;
Feldlaenge              = 300;

Oberflaeche             = zeros(Feldbreite+1,Feldlaenge+1);  % Oberfläche
Normierte_Oberflaeche   = zeros(Feldbreite+1,Feldlaenge+1);  % Oberfläche
Emitter                 = zeros(Feldbreite+1,Feldlaenge+1);  % Emittiende Oberfläche


% LED-Kennlinien von 0 bis 90° in 5°-Schritten
% Muss für jede LED individuell angegeben werden
% Beispielwerte für EDISON 3 Watt Emitter EDEV-SLC1-03
% Wenn nicht benötigt die beiden oberen Zeilen auskommentieren und die
% Zeilen darunter aktivieren

                              %0     5     10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90                                 
LED_Intensitaetskennlinie_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];
LED_Intensitaetskennlinie_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];

%LED_Intensitaetskennlinie_X = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00];
%LED_Intensitaetskennlinie_X = [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00];

% Abstand der Bestrahlten Oberfläche zum Strahler, mm
Abstand                  =   100;

% Hier die Emitter Platzieren. Es können beliebige Planare Emitterformen
% verwendet werden, allerdings wird die Rechnung SEHR langsam bei
% verwendung von Punkt oder Flächenemittern oder vielen Emittern

Emitter(80,80)           = 100;
Emitter(220,80)          = 100;
Emitter(80,220)          = 100;
Emitter(220,220)         = 100;

% Es wird die Emitierende Platte Punkt für Punkt abgefahren, und zu jedem
% Punkt die Rel. Lichtintensität auf der Oberfläche berechnet

% Abfahren des Emitters
for Emitter_X_Achse = 1 : Feldbreite+1
    for Emitter_Y_Achse = 1 : Feldlaenge+1
        % Strahlintensität am Aktuellen Emitterpunkt
        Emitterintensitaet = Emitter(Emitter_X_Achse, Emitter_Y_Achse);
        % Abfahren der Oberfläche
        if (Emitterintensitaet > 0)
            for Oberflaeche_X_Achse = 1 : Feldbreite+1
                for Oberflaeche_Y_Achse = 1 : Feldlaenge+1
                    % Berechnen des Abstandes zwischen den entspr. Punkten
                    deltaX = abs(Oberflaeche_X_Achse - Emitter_X_Achse);
                    deltaY = abs(Oberflaeche_Y_Achse - Emitter_Y_Achse);
                    % Berechnen der Wegstrecke (=Kugelradius)
                    Wegstrecke = (deltaX^2+deltaY^2+Abstand^2)^(0.5);
                    % Berechnen der Intensität, bedingt durch Entfernung
                    RelIntensitaet = (Abstand / Wegstrecke)^2;
                    % Berechnen des Aktuellen Winkels zwischen Emitter und
                    % Oberflächenpunkt
                    X_Winkel = uint8(atand(deltaX/Abstand)/5); 
                    Y_Winkel = uint8(atand(deltaY/Abstand)/5);
                    % Berechnen der Strahlungsstärke des Emitters für den
                    % aktuellen Winkel
                    LED_X_Intensitaet = LED_Intensitaetskennlinie_X(X_Winkel+1);
                    LED_Y_Intensitaet = LED_Intensitaetskennlinie_Y(Y_Winkel+1);
                    % Berechnen der Gesammtintensität
                    Intensitaet = RelIntensitaet * Emitterintensitaet * LED_X_Intensitaet * LED_Y_Intensitaet;
                    % Hinzufügen der Intensität des Emitterpunktes zum
                    % Oberflächenpunkt
                    Oberflaeche(Oberflaeche_X_Achse, Oberflaeche_Y_Achse) = Oberflaeche(Oberflaeche_X_Achse, Oberflaeche_Y_Achse) + Intensitaet;
                end
            end
        end
    end
end


% Normieren und Analysieren

Maximum = max(max(Oberflaeche));
Minimum = min(min(Oberflaeche));

for Oberflaeche_X_Achse = 1 : Feldbreite+1
    for Oberflaeche_Y_Achse = 1 : Feldlaenge+1
        Normierte_Oberflaeche(Oberflaeche_X_Achse, Oberflaeche_Y_Achse) = Oberflaeche(Oberflaeche_X_Achse, Oberflaeche_Y_Achse) / Maximum;
    end
end

% Auf der Oberfläche eine €-Platine plazieren
Platine = zeros(160,100);
for Platine_X_Achse = 1 : 100
    for Platine_Y_Achse = 1 : 160
            Platine(Platine_Y_Achse, Platine_X_Achse) = Normierte_Oberflaeche(Platine_X_Achse+100,Platine_Y_Achse+70);
    end
end

% Bestimmen des Minimums & des Maximums auf der Platine
Maximum = max(max(Platine))
Minimum = min(min(Platine))
% Die Homogenität auf der Oberfläche der €-Platine
Homogenitaet = (Minimum / Maximum)

newplot 
surface(Platine);
%surface(Normierte_Oberflaeche); % auskommentieren für Anzeige der
%Gesammten Berechneten Bestrahlten Oberfläche (300x300mm)
