% Script to validate results produced by anderson.m MATLAB function (R.E. Baddour, 2004).
% Reproduce Figures 3 and 4 from Anderson, V.C., JASA 22:426-431, 1950.

theta = 0; % backscatter direction
n = 25;    % 25 terms of series solution
r = 100;   % observation point, far from sphere

% Fig. 3
k = linspace(0,16,2000);
a = 1; % to reproduce ka = 0 to 16 as in Fig. 3
[Power_A,p_A,R_A] = anderson_EB(k,a,0.8,1,1,1,theta,n,r); % g = 1.0, h = 0.8
[Power_A,p_B,R_B] = anderson_EB(k,a,1.2,1,1,1,theta,n,r); % g = 1.0, h = 1.2
figure;
plot(k*a,R_A,k*a,R_B);
legend('curve A - g = 1.0, h = 0.8', 'curve B - g = 1.0, h = 1.2');
xlabel('ka');
ylabel('R');
title('Fig. 3 in Anderson, V.C., JASA 22:426-431, 1950');

% Fig. 4
k = linspace(0,16,2000);
a = 1; % to reproduce ka = 0 to 6.5 as in Fig. 4
[Power_A,p_A,R_A] = anderson_EB(k,a,2,1,2,1,theta,n,r);       % g = 2.0, h = 2.0
[Power_B,p_B,R_B] = anderson_EB(k,a,10^8,1,10^8,1,theta,n,r); % g and h approach infinity (rigid sphere)
[Power_C,p_C,R_C] = anderson_EB(k,a,0.5,1,0.5,1,theta,n,r);   % g = 0.5, h = 0.5
figure;
plot(k*a,R_A,k*a,R_B,k*a,R_C);
legend('curve A - g = 2.0, h = 2.0', 'curve B - rigid sphere', 'curve C - g = 0.5, h = 0.5');
xlabel('ka');
ylabel('R');
title('Fig. 4 in Anderson, V.C., JASA 22:426-431, 1950');

figure;
plot(k*a,Power_A,k*a,Power_B,k*a,Power_C);
legend('curve A - g = 2.0, h = 2.0', 'curve B - rigid sphere', 'curve C - g = 0.5, h = 0.5');
xlabel('ka');
ylabel('Power');
title('Fig. 5 in Anderson, V.C., JASA 22:426-431, 1950');
