function [  ] = CompareRAPDFs(  )
% Comparison of Reflection Angle Probability Density Functions (RAPDFs)

% This GUI allows the computation of the Reflection Angle Probability
% Density Functions (RAPDFs)for both the Vecor Mixing (VM) and the Hybrid
% Reflectance Model (HRM)

    %% Predefines (angles in degrees)
    thetaIn      =   -20;
    s            =   0.7;
    thetaOutList =   -89.997:0.1:89.5;
    Tol          =   1e-8;

    BackgroundColor = [0.7 0.7 0.7];
    ForegroundColor = [0.8 0.8 0.8];
    HighlightColor = [0.6 0.6 0.8];

    %% GUI
    h_Figure=figure('Position',[100 100 600 400],'NumberTitle','off','Name','filename','Color',BackgroundColor);

    % axes
    h_AnglePlot     = axes('Position',[0.05  .25  .9  .7]);

    % select thetaIn
                   uicontrol('Position', [ 10 25 100 20],  'Style', 'text',        'String', 'Incidence Angle [°]',  'BackgroundColor',  BackgroundColor);
    g_thetaIn =    uicontrol('Position', [110 25 25 20],   'Style', 'edit',        'String', num2str(thetaIn)     ,  'BackgroundColor',  HighlightColor);
    
    % select s
                uicontrol('Position', [135 25 100 20],  'Style', 'text',        'String', 'Scat. Coeff.',        'BackgroundColor',  BackgroundColor);
    g_s   =     uicontrol('Position', [235 25 25 20],   'Style', 'edit',        'String', num2str(s),          'BackgroundColor',  HighlightColor);
    
    % select tolerance
                uicontrol('Position', [260 25 50 20],   'Style', 'text',        'String', 'Tolerance',       	 'BackgroundColor',  BackgroundColor);
    g_Tol =     uicontrol('Position', [320 25 50 20],   'Style', 'edit',        'String', num2str(Tol),         'BackgroundColor',  HighlightColor);
    
    % compute
                uicontrol('Position', [385 25 100 20],  'Style','pushbutton',   'String',   'Compute'    ,     'Callback', @compute,       'BackgroundColor',  ForegroundColor);
                
    % hold
    h_hold =    uicontrol('Position', [500 25 45 20],  'Style', 'checkbox',     'String','Hold',                                            'BackgroundColor',  BackgroundColor);

    function[] = compute(hObj,event)
        % Clean graph
        if (get(h_hold,'value')==0)
            % clear Axes
            cla(h_AnglePlot);
        end
        
        % hold
        hold on;
        
        % read config
        Tol     = (str2num(get(g_Tol,   'String')));
        s       = (str2num(get(g_s,  'String')));
        thetaIn = (str2num(get(g_thetaIn,  'String')));
        
       
        %% VM
        % normalization constant
        R0_VM = 1./RAPVM(-pi/2,pi/2,thetaIn*pi/180,s,1,Tol);
        
        % compute Distribution
        RAPDF_VM = RAPDFVM(thetaOutList*pi/180, thetaIn*pi/180, s, R0_VM);
        
        % plot
        plot(thetaOutList,RAPDF_VM,'color','blue','LineWidth',2.0);
        
        
        %% HRM
        % normalization constant
        R0_HRM = 1./RAPL(-pi/2,pi/2,1,Tol);
        
        % compute distribution
        RAPDF_HRM = s * RAPDFL(thetaOutList*pi/180,R0_HRM);
        
        % plot
        plot(thetaOutList,RAPDF_HRM,[-thetaIn -thetaIn],[0 (1-s)],'color','red','LineWidth',2.0);
        legend('Vector Mixing (VM)', 'Hybrid Reflectance Model (HRM)');
        
        xlim([-95 95]);   
    end

    function[RI] = RAPL(theta_min, theta_max, R_0, Tol)  
        % Integration of distribution for normalization (HRM, Lambert part)
        
        RI=quad(@RAPDFL, theta_min, theta_max,Tol,[],R_0);
    end

    function[R] = RAPDFL(theta, R_0)
        % Compute distibution (HRM, Lambert part)
        
        R = R_0 .* cos(theta);
    end

function[RI] = RAPVM(thetaOut_min, thetaOut_max, thetaIn, s, R_0, Tol) 
        %  Integration of distribution for normalization (VM)
        
        % left most and right most vector ends possible
        xMin = -sin(thetaIn)*(1-s) - s;
        xMax = -sin(thetaIn)*(1-s) + s;
        yMinMax =  cos(thetaIn)*(1-s);
        
        % corresponding angles
        StartAngle = atan(xMin/yMinMax);
        EndAngle = atan(xMax/yMinMax);
        
        % adjustion of angles in case of small scattering values
        if(s<0.5)
            
            % specular vector end
            xM = sin(thetaIn)*(1-s);
            yM = yMinMax;
        
            a = (xM^2*yM^2-(yM^2-s^2)*(xM^2-s^2))/(xM^2-s^2)^2;
            m1 = -xM*yM/(xM^2-s^2)+sqrt(a);
            m2 = -xM*yM/(xM^2-s^2)-sqrt(a);

            L1 = atan(1/m1);
            L2 = atan(1/m2);

            if(L1>0&&StartAngle>L1)
                StartAngle=L1;
            elseif(L2<0&&EndAngle<L2)
                EndAngle=L2;
            end
        end
        
        RI=quad(@RAPDFVM, max(thetaOut_min,StartAngle), min(thetaOut_max,EndAngle),Tol,[],thetaIn,s,R_0);
end

function[R] = RAPDFVM(thetaOut, thetaIn, s, R_0)
     % Compute distibution (VM)
     
     % Find possible Lambert angles
    [T1, T2] = ComputeLambertFromOut(thetaIn,thetaOut,s);

    % For each theta
    for i=1:length(thetaOut)
        % If only T2 valid --> result=T2
        if(isnan(T1(i))&&isfinite(T2(i)))
            R(i)=R_0*cos(T2(i));
        % If only T1 valid --> result=T1
        elseif(isfinite(T1(i))&&isnan(T2(i)))
            R(i)=R_0*cos(T1(i));
        % If both T valid, add up result    
        elseif(isfinite(T1(i))&&isfinite(T2(i)))
            R(i)=R_0*cos(T1(i))+R_0*cos(T2(i));
        % If no T valid --> result=0
        elseif(isnan(T1(i))&&isnan(T2(i)))
            R(i)=0;
        else
            disp('strange behaviour')
        end
    end
end

function [ o_ThetaLambert1, o_ThetaLambert2] = ComputeLambertFromOut( i_ThetaIn, i_ThetaResult, i_S)
    % Find nescessary Lambert angle, such that resulting angle is computed
    % (up to two results possible)
    p_ThetaOut      = - i_ThetaIn;

    p_a = (1-i_S)/i_S;

    % temp values
    p_sTO = cos(p_ThetaOut);
    p_cTO = -sin(p_ThetaOut);

    p_tA = -cot(i_ThetaResult);
    
    % for each resulting angle
    for i=1:length(i_ThetaResult)
        
        % undefined, y=0
        if(abs(i_ThetaResult(i))<1e-6)
            k = -p_cTO*p_a;
            if(abs(k)>=1)
                o_ThetaLambert1(i) = NaN;
                o_ThetaLambert2(i) = NaN;
            else
                o_ThetaLambert1(i) = acos(k)-pi/2;
                o_ThetaLambert2(i) = acos(k)-pi/2;
            end

        else
            p_Value1 =cot(i_ThetaResult(i))^2*p_a^2*cos(p_ThetaOut)^2+cot(i_ThetaResult(i))^2-p_a^2*cot(i_ThetaResult(i))^2+1-p_a^2+2*p_a^2*cos(p_ThetaOut)*cot(i_ThetaResult(i))*sin(p_ThetaOut)+sin(p_ThetaOut)^2*p_a^2;

            if(p_Value1>0)
                p_L1 = atan2(-p_tA(i)*(p_tA(i)^2*p_a*p_cTO-p_tA(i)*p_a*p_sTO-sqrt(p_Value1))/(p_tA(i)^2+1)-p_a*p_sTO+p_a*p_tA(i)*p_cTO, -(p_tA(i)^2*p_a*p_cTO-p_tA(i)*p_a*p_sTO-sqrt(p_Value1))/(p_tA(i)^2+1))-pi/2;
                p_L2 = atan2(-p_tA(i)*(p_tA(i)^2*p_a*p_cTO-p_tA(i)*p_a*p_sTO+sqrt(p_Value1))/(p_tA(i)^2+1)-p_a*p_sTO+p_a*p_tA(i)*p_cTO, -(p_tA(i)^2*p_a*p_cTO-p_tA(i)*p_a*p_sTO+sqrt(p_Value1))/(p_tA(i)^2+1))-pi/2;
                if (p_L1<-pi/2||p_L1>pi/2)
                    o_ThetaLambert1(i) = NaN;
                else
                    o_ThetaLambert1(i) = p_L1;
                end
                if (p_L2<-pi/2||p_L2>pi/2)
                    o_ThetaLambert2(i) = NaN;
                else
                    o_ThetaLambert2(i) = p_L2;  
                end

            else
                o_ThetaLambert1(i)=NaN;
                o_ThetaLambert2(i)=NaN;
            end
        end
    end
    
end

    
end

