%% Total sound field of a difraction problem at a rectagular wedge using UTD

% Author: Jonas Stienen


%% Scene

n1  = [ 0  1  0 ]; % main face 1
n2  = [ 1  0  0 ]; % opposite face 2 / facade
n3  = [ 0  1  0 ]; % ground plane / face 3
loc = [  0  0  0 ];
w = itaInfiniteWedge( n1, n2, loc );

r    = [  2  -2  0 ];
r2   = [ -2  -2  0 ]; % image receiver 2
r3   = [  2  -4  0 ]; % image receiver 3
r23  = [ -2  -4  0 ]; % image receiver 23 /  32
s1  = 5 * [ 0  -1  0 ];

a = w.apex_point( s1, r );
ad = w.edge_direction;
c = 344; % speed of sound

f = [ 20, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000 ]';
f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 Hz', '800 Hz', ...
    '1.6 kHz', '3.2 kHz', '6.4 kHz', '12.8 kHz', '24 kHz' };

%f = [ 20, 50, 100, 200, 400 ]';
%f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 Hz' };

num_angles = 1000;
alpha_d_start = 0;
alpha_d_end = pi;
alpha_d = linspace( alpha_d_start, alpha_d_end, num_angles );

% Set different receiver positions rotated around the aperture
S = norm( s1 ) * [ - cos( alpha_d ); sin( alpha_d ); zeros( 1, numel( alpha_d ) ) ]';

% wavenumber
k = 2 * pi * f ./ c;

% store results in ita audio format
res_data_template = itaResult();
res_data_template.resultType = 'signal';
res_data_template.freqVector = f;


%% Sound field simulation

N = size( S, 1 );
for n = 1 : N
    
    s = S( n, 1:3 );
    
    % Zero order
    
    r_dir = norm( s - r );
    H_i = 1 ./ r_dir .* exp( -1i .* k .* r_dir ); % direct sound component
    
    r_1_dir = norm( s - r2 );
    H_r = 1 ./ r_1_dir .* exp( -1i .* k .* r_1_dir ); % reflected sound component

    H_d = ita_diffraction_utd( w, s, r, f, c ); % diffracted sound component
    
    % Ground reflection image
    
    r_dir_1 = norm( s - r3 );
    H_i_1 = 1 ./ r_dir_1 .* exp( -1i .* k .* r_dir_1 ); % direct sound component
    
    r_1_dir_1 = norm( s - r23 );
    H_r_1 = 1 ./ r_1_dir_1 .* exp( -1i .* k .* r_1_dir_1 ); % reflected sound component

    [ H_d_1 , D_half ] = ita_diffraction_utd( w, s, r3, f, c ); % diffracted sound component
    
    % zero and first order components
    
    H = H_d; % Total sound field ... diffracted component always present
    if ~ita_diffraction_shadow_zone( w, s, r )
       H = H + H_i;
    end
    if ita_diffraction_reflection_zone( w, s, r, false ) % opposite face reflection
        H = H + H_r;
    end
    
    H_1 = H_d_1;
    if ~ita_diffraction_shadow_zone( w, s, r3 )
        H_1 = H_1 + H_i_1;
    end
    if ita_diffraction_reflection_zone( w, s, r3, false ) % opposite face reflection
        H_1 = H_1 + H_r_1;
    end
    
    H_1_half = zeros( numel( H_d_1 ), 1 );
    if ita_diffraction_shadow_zone( w, s, r3 )
        H_1_half = H_d_1; % Shadow zone yes diffraction only
    else
        H_1_half = H_i_1; % Direct sound yes
        if ~ita_diffraction_reflection_zone( w, s, r3, false )
            if ita_diffraction_utd_illumination_sign( w, s, r3, false ) == -1
                sf = ita_diffraction_utd_illumination_fadeout( w, s, r3, false );
                H_1_half = H_1_half + sf .* H_d_1; % closer to shadow zone diffraction yes
            end
        else
            humbug = 1;
        end
    end
    
    att_sum = res_data_template;
    att_sum.freqData = ( H_1_half ) ./ H_i; % Insertion loss
        
    if n == 1
        att_sum_all = att_sum;
    else
        att_sum_all = ita_merge( att_sum_all, att_sum );
    end
    
end


%% Plot

fg = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );

h = plot( rad2deg( alpha_d ), flip( att_sum_all.freqData_dB', 2 ) );

title( 'Transfer function nomalized to free-field transmission' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )

lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;

xlim( [ -5 185 ] );
ylim( [ -35 10 ] );
grid on;

lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
    h( i ).LineWidth = lw( i );
end


%% Export

set( fg, 'papersize', [ 16*1.5 9*1.5 ], 'units', 'centimeters' )
print( 'ita_diffraction_plot_utd_total_field_stienen_diss_with_ground', '-fillpage', '-dpdf' )