%% Schaerfentieferechner
% Erstellt wird ein einfaches Diagramm, in dem die Schaerfentiefe fuer
% verschiedene Brennweiten und Entfernungen abgelesen werden kann.
% Ergebnis ist ein druckbares PNG-Bild
%
% Verwendung findet die einfache Schaerfentiefeformel.
%
% Getestet mit Matlab R2014b und Octave 4.0.0
%
% 10.01.2016 DL1DOW


clear all
close all


%% Einstellungen
% Im Struct "Cfg" koennen alle Einstellungen gemaess der verwendeten 
% Kamera-Objektiv-Kombination geaendert werden

switch 'Canon_APS-C'
    case 'Sensor'

        % Berechnung sinnvoller Zerstreuungskreis aus 
        % physikalischer Sensoraufloesung
        spx = [5472 3648]; % [Pixel]
        smm = [22.5 15]; % [mm]
        ressens = max(spx);
        Z = sqrt(smm(:)'*smm(:))/ressens;
        
    case 'SensorFausformel'

        % Berechnung sinnvoller Zerstreuungskreis aus 
        % Sensordiagonale/1500
        smm = [22.5 15]; % [mm]
        Z = sqrt(smm(:)'*smm(:))/1500;

    case 'Canon_APS-C'
        Z = 0.018;

end
        
% Erlaubter Zerstreuungskreis [mm]
% z.Z. 0.018mm fuer APS-C-Sensor
Cfg.circleOfConfusion = Z;

% Brennweite-Blenden-Kombinationen 
% In Struct "Lens" ist fuer jede gewuenschte Brennweite ein Eintrag. In 
% jedem Brennweiteneintrag ist die gewuenschte Blendenreihe
% Nur ganze Blendenwerte (k = [1 1.2 2 2.4 4 5.6 8 11 16 22]) werden im
% Diagramm farbkodiert.

switch '120-400mm'
    case '18-55mm'
        % Name fuer Diagrammueberschrift
        Cfg.name = 'EF-S 18-55mm f/3.5-5.6';


        % Brennweite in [mm], 
        Cfg.Lens(1).focalLen = 18; % [mm]
        Cfg.Lens(2).focalLen = 55; % [mm]

        % Blende in F-Stops
        Cfg.Lens(1).fStop = [4.5 5.6 8 11 16 22];
        Cfg.Lens(2).fStop = [5.6 8 11 16 22];

        % Naheinstellgrenze [mm]
        Cfg.Lens(1).minFocus = 390;
        Cfg.Lens(2).minFocus = 390;

    case '120-400mm'
        % Name fuer Diagrammueberschrift
        Cfg.name = '120-400mm f/3.5-5.6';


        % Brennweite in [mm], 
        Cfg.Lens(1).focalLen = 120; % [mm]
        Cfg.Lens(2).focalLen = 400; % [mm]

        % Blende in F-Stops
        Cfg.Lens(1).fStop = [4.5 5.6 8 11 16 22];
        Cfg.Lens(2).fStop = [5.6 8 11 16 22];

        % Naheinstellgrenze [mm]
        Cfg.Lens(1).minFocus = 1500;
        Cfg.Lens(2).minFocus = 1500;
end        





%% Einstellungen pruefen und uebernehmen
assert(ischar(Cfg.name),' ''Cfg.name'' must be of class ''char'' ');


for flen = 1:numel(Cfg.Lens)
    
    % Zerstreuungskreis
    Z = Cfg.circleOfConfusion;
    assert(isfloat(Z) && numel(Z) == 1,...
        ' ''circleOfConfusion'' must be scalar numerical value');
    if Z > 1
        warning(' ''circleOfConfusion'' is larger than 1 mm ! ');
    end
    
    % Brennweite
    f = Cfg.Lens(flen).focalLen;
    assert(isfloat(f) && numel(f) == 1,...
        ' ''focalLen'' must be scalar numerical value');
    if f <= 0
        warning(' ''focalLen'' must be positive scalar value ');
    end
    
    % Blendenreihe
    % liegender Vektor
    k = Cfg.Lens(flen).fStop(:)';
    assert(isfloat(k),...
        ' ''fStop'' must be numerical value');
    
    
    % Naheinstellgrenze
    d_np = Cfg.Lens(flen).minFocus;
    

    %% Berechnungen
    
    % Hyperfokale Entfernung [mm]
    % Zeilenvektor: d_h fuer alle k
    d_h = f.^2./(k * Z)+f;

    
    % Berechnete Gegenstandsweiten sind logarithmisch zwischen der
    % Naheinstellgrenze und der hyperfokaten Distanz verteilt.
    %
    % G: alle berechneten Gegenstandsweiten fuer alle Blendenzahlen
    % Zeilen: Verschiedene Gegenstandsweiten (Schrittweite in Abhaengigheit
    %          von der hyperfokaten Distanz)
    % Spalten: Verschiedene Blendenzahlen
    nDp = 100; % Anzahl der Datenpunkte zwischen Naheinstellgrenze und
               % hyperfokalen Entfernung
    G = zeros(nDp,numel(k));  
    for i = 1:numel(k)

        d_hi = d_h(i);

        % Gegenstandsweiten [mm]
        g = logspace(log10(d_np),log10(d_hi),nDp)';
        G(:,i) = g;
    end
    
    
    % Hilfswert: Einsen zum f�llen, so dass gilt: size(G) == size(onef*k)
    onef = 0*G(:,1)+1;


    % Nahpunkt [mm]
    % Zeilen: unterschiedliche Gegenstandsweiten
    % Spalten: unterschiedliche Blendenzahlen
    % d_n = g*(d_h - f)/(d_h +g -2*f) 
    d_n = ( G.*(onef*d_h -f ) )./( onef*d_h + G -2*f);



    % Fernpunkt [mm]
    % d_f = g*(d_h-f)/( d_h - g)
    % d_f = inf, wenn g > d_h
    % da die letzte Zeile von G die Hyperfokaldistanz fuer die jeweilige
    % Brennweite enthaelt, muss die letzte Zeile von d_f zwangsweise
    % Unendliche sein.
    d_f = ( G.*(onef*d_h -f ) )./( onef*d_h - G);
    d_f( end,:) = Inf;


    % Schaerfentiefe [mm]
    dd = d_f-d_n;



    %% Ploteinstellungen
    % Skalierungs des Plots, [mm] -> [m]
    pltscale = 1000;

    % Darstellung der einzelnen Blendenzahlen, kann beliebig ergaenzt werden
    Kc = { 1,'y','-';
         1.2,'y','--';
           2,'g','-';
         2.4,'g','--';
           4,'b','-';
         5.6,'b','--';  
           8,'r','-';
          11,'r','--';
          16,'k','-';
          22,'k','--';
         };

    % Stiftbreite
    linewidth = 1.5;
 
    %% Grafische Darstellung/Plotten: Elemente fuer jede Brennweite

    
    if flen==1
        % Dummy-Plot
        dummy = loglog(1,1);        
        hold on
        
        % Hintergrundlinien fuer mm, cm, m, km (mit Dummywerten)
        Hgl(6) = plot([1 1],[1 1]*1e0, 'k:','LineWidth',linewidth);
        Hgl(5) = plot([1 1],[1 1]*1e3, 'k:','LineWidth',linewidth);
        Hgl(4) = plot([1 1],[1 1]*1e-3,'k:','LineWidth',linewidth);
        Hgl(3) = plot([1 1]*1e0, [1 1],'k:','LineWidth',linewidth);
        Hgl(2) = plot([1 1]*1e3, [1 1],'k:','LineWidth',linewidth);
        Hgl(1) = plot([1 1]*1e-3,[1 1],'k:','LineWidth',linewidth);     
    end
    
    % Darstellung der Schaerfentiefe-Verlaeufe
    % Sowohl G als auch dd enthalten die Verlaeufe in der Reihenfolge wie k
    H = loglog(G(:,:)/pltscale,dd(:,:)/pltscale,...
        'k-','LineWidth',linewidth);


    % linke Plotgrenze: Auf ganze Zehnerpotenz unterhalb Naheinstellgrenze
    % runden
    % rechte Plotgrenze: 1 km
    xmin = 10^floor(log10(d_np/pltscale));
    xmax = 1000;
    axis([xmin xmax 5e-4 5e3]);


    hold on
    grid on
    objname = Cfg.name;
    title(['Schaerfentiefe ' objname '   Z=' num2str(Z) 'mm'],'FontSize',14);
    xlabel('Gegenstandsweite g [m]');
    ylabel('Schaerfentiefe [m]');

    
    % Legende
    %legend(strcat('F/',num2str(k')),'Location','NorthWest');

    % Beschriftung: Alle Blendenzahlen direkt an Plot
    Kt = strcat('F/',num2str(k'),' |'); % Beschriftungstext
    o = 0*k+1;                          % Hilfsvektor mit Einsen
    tx = d_np/pltscale*o;               % x-Positionen
    ty = dd(1,:)/pltscale;              % y-Positionen
    T = text(tx,ty,Kt,...
        'HorizontalAlignment','right');


    
    % Textfeld fuer Brennweitenbeschriftung
    xlim = get(gca,'XLim');
    bx = xlim(1)*1.2;
    by = [min(dd(1,:)) max(dd(1,:))]/pltscale;
    txt = [num2str(f) ' mm'];
    text(bx,min(by),txt,...
        'HorizontalAlignment','left','VerticalAlignment','top',...
        'Rotation',90,'BackgroundColor','w','FontSize',12,...
        'EdgeColor','k');

    
    % Plots fuer die einzelnen Blendenzahlen entsprechen der obigen Tabelle
    % formatieren
    for i = 1:size(Kc,1)
        ki = Kc{i,1};
        ci = Kc{i,2};
        li = Kc{i,3};
        idx = (k==ki);
        set(H(idx),'Color',ci,'LineStyle',li);
        set(T(idx),'Color',ci);
    end

    
    % Hilfslinien fuer mm, cm, m, km
    text(1e3,1e0 ,' [m]');
    text(1e3,1e3 ,' [km]');
    text(1e3,1e-2,' [cm]');
    text(1e3,1e-3,' [mm]');


    % Hintergrundlinien korrekt anordnen
    xlim = get(gca,'XLim');
    ylim = get(gca,'YLim');
    set(Hgl(6),'XData',xlim);
    set(Hgl(5),'XData',xlim);
    set(Hgl(4),'XData',xlim);
    set(Hgl(3),'YData',ylim);
    set(Hgl(2),'YData',ylim);
    set(Hgl(1),'YData',ylim);
    
    
    
end % Ende Schleife ueber alle Brennweiten    
    


% Verlauf drucken
filename = [Cfg.name '.png'];
notallowed = '\/?".*:|';
for i = 1:numel(notallowed)
    filename(filename == notallowed(i)) = ' ';
end
disp(filename)
print(filename,'-dpng','-r300');


