function [SHR_outwave] = shoulder_reflect(theta, elevation, fs, ITD_out_L, ITD_out_R)
% Calculates Shoulder reflections from the angle theta (horizontal) and angle of elevation as defined in the script script.m which calls this function.

% Inputs are:

% ITD_out_L - mono signal that has been processed with reverb
% ITD_out_R - mono signal that has been processed with reverb

% Variables 

% theta - angle on the horizontal plane, as declared in the script that
% calls this function
% elevation - angle of elevation, as declared in the script that calls this
% function
% fs - sampling frequency to calculate samples, which is read from the
% inputted waveform with wavread in the script that calls this function

% Outputs

% SHR_outwave - a stereo file which is the final output, and is played back
% and written to a wav file in the script that calls this function.

% Generic Example of how to call

% [SHR_outwave] = shoulder_reflect(-35, 45, 44100, input_l, input_r);

% Sally-anne Kellaway - skel2480@uni.sydney.edu.au // Based on script from
% DAFX - Rocchesso, D 2002, 'Spatial Effects', in U Zolzer (ed.)

% ----------------------------------------------------------------------


one = 180/(180+theta);                            % steps out the calculation of the shoulder reflection delay in to 4 sequential calculations
two = ((elevation-80)*one)^2;
three = 1-0.00004*two;
four = 1.2*((180-theta)/180);

sh_delay_ms = four*three;                         % arrives at the delay in miliseconds
sh_delay_sec = sh_delay_ms/1000;                  % converts to seconds

SHR_delay_samples = round(fs*sh_delay_sec);       % converts to samples

extra_samples = zeros(SHR_delay_samples, 1);      % Create a vector of zero samples the length as calculated, to use as a delay

SHR_L = [extra_samples; ITD_out_L];               % Create a delayed version by adding the zeros to the filtered head shadow 
SHR_R = [ITD_out_R; extra_samples];               % Add the delay to the end of the undelayed channel to ensure both files are the same length

% this section adjusts the gains of the files to prepare for summation of the ITD and SHR...
% 2 x files = halve the gain of each

gain = 0.5;
SHR_L = SHR_L*gain;
SHR_R = SHR_R*gain;
ITD_out_L = ITD_out_L*gain;
ITD_out_R = ITD_out_R*gain;


% this section assigns an overall value for the size/length of the ITD and
% SHR files
if size(ITD_out_L,1) == size(ITD_out_R,1);
    
    ITD_size = size(ITD_out_L,1);
 
end;

if size(SHR_L,1) == size(SHR_R,1);
    
    SHR_size = size(SHR_L,1);
    
end;

% this section tests whether the ITD or SHR results in the longer signal
% and adjusts the opposite with zeros on the end to ensure both are the
% same length for addition

if SHR_size > ITD_size

    add_ITD = SHR_size-ITD_size;
    ITD_adjusted_L = [ITD_out_L;zeros(add_ITD,1)];      % adjusts the signal with ITD to suit the length of delayed signal
    ITD_adjusted_R = [ITD_out_R;zeros(add_ITD,1)];
    
    SHR_out_L = SHR_L + ITD_adjusted_L;                 % add the ITD output signals and SHR output signals together
    SHR_out_R = SHR_R + ITD_adjusted_R;

else
    
    add_SHR = ITD_size-SHR_size;
    SHR_adjusted_L = [SHR_L;zeros(add_SHR,1)];          % adjusts the signal with SHR to suit the length of delayed signal
    SHR_adjusted_R = [SHR_R;zeros(add_SHR,1)];

    SHR_out_L = SHR_adjusted_L + ITD_out_L;             % add the ITD output signals and SHR output signals together
    SHR_out_R = SHR_adjusted_R + ITD_out_R;
    
end;

    SHR_out_L = 0.99.*(SHR_out_L./max(abs(SHR_out_L)));  % normalise
    SHR_out_R = 0.99.*(SHR_out_R./max(abs(SHR_out_R)));

SHR_outwave = [SHR_out_L, SHR_out_R];                   % provides a stereo file for testing 