% go_validate 
% 
% MATLAB script to find the spatial coherence of far-field 
% echos from beryllium and aluminum spheres. 
% This script is provided to validate the model implementation 
% against a number of published results. 
% Martin Anderson 
% Created: 9-2-98  Last modified: 10-29-99 
% 
% References are to Hickling (JASA 34: 1582-1592, 1962) 
% and Hampton and McKinney (JASA 33: 664-673, 1961). 
% =============================================== 

% Beryllium backscatter over spatial frequency: 
% Compare to Figure 3 in Hickling. 

tic 

k = .1:.1:30; 
a = 1; 
c1 = 12890000; 
c3 = 1480000; 
rho1 = 1.87; 
rho3 = 1; 
sigma = 0.05; 
n = 40; 
theta = pi; 

[fe,fr] = faran(k,a,c1,c3,rho1,rho3,sigma,theta,n); 

save beryl fe fr k a c1 c3 rho1 rho3 sigma theta n 

% =============================================== 

% Aluminum backscatter over spatial frequency: 
% Compare to Figures 7 and 8 in Hickling. 

k = .1:.1:30; 
a = 1; 
c1 = 6420000; 
c3 = 1480000; 
rho1 = 2.70; 
rho3 = 1; 
sigma = 0.355; 
n = 40; 
theta = pi; 

[fe,fr] = faran(k,a,c1,c3,rho1,rho3,sigma,theta,n); 

save alum fe fr k a c1 c3 rho1 rho3 sigma theta n 

toc 

% =============================================== 

load beryl 

clf 
subplot(2,1,1) 
plot(k*a, abs(fe), k*a, abs(fr), '--'); 
title('Far field magnitude function using elastic model for beryllium sphere') 
ylabel('|f_{\infty}|') 
ylabel('Magnitude') 
xlabel('ka') 
grid 

subplot(2,1,2) 
plot(k*a, unwrap(-angle(fe))./(k*a), k*a,  unwrap(-angle(fr))./(k*a), '--'); 
title('Far field phase function using elastic model for beryllium sphere') 
ylabel('-arg(f_{\infty})/ka') 
ylabel('-Phase/ka [rad]') 
xlabel('ka') 
grid 
orient tall 
print beryl.ps 

% =============================================== 

load alum 

clf 
subplot(2,1,1) 
plot(k*a, abs(fe), k*a, abs(fr), '--'); 
title('Far field magnitude function using elastic model for aluminum sphere') 
ylabel('|f_{\infty}|') 
ylabel('Magnitude') 
xlabel('ka') 
grid 

subplot(2,1,2) 
plot(k*a, unwrap(-angle(fe))./(k*a), k*a,  unwrap(-angle(fr))./(k*a), '--'); 
title('Far field phase function using elastic model for aluminum sphere') 
ylabel('-arg(f_{\infty})/ka') 
ylabel('-Phase/ka [rad]') 
xlabel('ka') 
grid 
orient tall 
print alum.ps 

% =============================================== 

% Aluminum scatter over angle, 1-inch and 3-inch spheres, 150 kHz: 
% Compare to Figure 10 in Hampton and McKinney. 

% Note that Hampton and McKinney place the sound source at 0 degrees 
% rather than at 180 degrees, which is the Faran convention. 
% Also note that Hampton's and McKinney's results include the 
% sound source response, and thus are only comparable for angles 
% within 120 degrees of backscatter. The description of their 
% apparatus suggests that some dropout may be present 
% in their measurements at the direct backscatter angle, 180 degrees 
% (0 degrees in their plots) due to overlap of the source and hydrophone. 
% Also note that the results presented here are monochromatic, 
% while Hampton and McKinney used a 500 microsecond gated pulse. 
% The dB reference used is also arbitrary to avoid negative 
% values on a rho-theta plot. 
% The rapid increase in oscillations in the spectra at higher 
% frequencies suggest that best agreement will be achieved at 
% relatively low ka values, such as for the 1-inch sphere here. 

tic 

a = (1/2)*25.6; 
c1 = 6420000; 
c3 = 1480000; 
k = 2*pi*150e3/c3; 
rho1 = 2.70; 
rho3 = 1; 
sigma = 0.355; 
n = 40; 

theta = [pi-120*pi/180:pi/80:pi+120*pi/180]; 

fe = zeros(size(theta)); 
fr = zeros(size(theta)); 

for ct = 1:length(theta) 
        [fe(ct),fr(ct)] = faran(k,a,c1,c3,rho1,rho3,sigma,theta(ct),n); 
end 
save alum_angle1 fe fr k a c1 c3 rho1 rho3 sigma theta n 

% =============================================== 
% Aluminum scatter over angle, 3-inch sphere, 150 kHz: 

a = (3/2)*25.6; 
c1 = 6420000; 
c3 = 1480000; 
k = 2*pi*150e3/c3; 
rho1 = 2.70; 
rho3 = 1; 
sigma = 0.355; 
n = 40; 

theta = [pi-120*pi/180:pi/80:pi+120*pi/180]; 

fe = zeros(size(theta)); 
fr = zeros(size(theta)); 

for ct = 1:length(theta) 
        [fe(ct),fr(ct)] = faran(k,a,c1,c3,rho1,rho3,sigma,theta(ct),n); 
end 
save alum_angle3 fe fr k a c1 c3 rho1 rho3 sigma theta n 

toc 

% =============================================== 

load alum_angle1 

clf 
subplot(2,1,1) 
db = 10*log10(fe.*conj(fe)); 
db = db-min(db); 
polar(theta, db); 
title('Far field magnitude function using elastic model for 1-inch aluminum sphere') 

load alum_angle3 

subplot(2,1,2) 
db = 10*log10(fe.*conj(fe)); 
db = db-min(db); 
polar(theta, db); 
title('Far field magnitude function using elastic model for 3-inch aluminum sphere') 

orient tall 
print alum_angle.ps 

