%n:\mat5work\sgn\glass_sphere_impact_paper.m   02/07/14

%sound generated by impacting spheres

%generates figure 4 in the paper

%This software is in support of the paper
%An overview of underwater sound generated by inter-particle collisions 
%and its application to the measurements of coarse sediment bedload transport
% by Peter D. Thorne


%1 The time domain signal for impactor, impactee and combined signal are
%calculated
%2 The Spectrum is calulated analytically with B=0, full analytically and
%numerically from the time series signal


% The program is for educational use and not for comercial purposes.
% The code is unsupported.


% set for B=0 and have particles of different size

close all
clear all


a1o=[5 15 50]*1e-3;  %impactee
a2o=[15 50 5]*1e-3;   %impactor

figure, orient tall, mn=15

laz='abcdef'

for jj=1:3

% 1- impactee   2- impactor

c=1500;  %vel sound in water
rho=1000; % density of the fluid
rhos=2500; %sphere density  stainless steel
a1=a1o(jj) %35e-3;  %sphere radius 1
a2=a2o(jj)    %35e-3;  %sphere radius 2
U=0.3; % impact velocity
theta=0;  %angle between line of motion and direction of field point
r=0.3; %range to filed point
e=7e10; sigma=0.2;%elasticity (youngs modulus) and poisson's ratio for stainless steel

m1=(4/3)*pi*rhos*a1.^3;
m2=(4/3)*pi*rhos*a2.^3;

Am1=2*m2/(m1+m2); Am2=-(m1/m2)*Am1;

to=4.53*((2*(m1*m2/(m1+m2))*(1-sigma^2)/(pi*e)).^0.4)*((a1+a2)/(a1*a2*U)).^0.2;  

eta1=c*to/(pi*a1);
eta2=c*to/(pi*a2);


td=a1*(1+pi/2)/c; %Time delay
%td=2*a1*cos(theta)/c;   %theta>0
f=(0:1:(2^19)-1);  %high freq required to get good time signal
fo=1/(2*to);
fn=f/fo;
omega=2*pi*f;
omegao=2*pi*fo;

%Frequency domain analytical full soln for any impacting particles
Po1=Am1*a1*(rho*c*U*cos(theta)/(2*r^2*omegao))*(1+exp(-i*pi*fn))./(1-fn.^2);% signal amplitude

%Impactee
xx1=(eta1*a1+i*r*fn);
xx2=2*eta1^2-fn.^2+2*i*eta1*fn;
Afa=Po1.*xx1./xx2;

%impactor
Po2=Am2*a2*(rho*c*U*cos(theta)/(2*r^2*omegao))*(1+exp(-i*pi*fn))./(1-fn.^2);% signal amplitude
xx1b=(eta2*a2+i*r*fn);
xx2b=2*eta2^2-fn.^2+2*i*eta2*fn;
Afb=Po2.*xx1b./xx2b;

Af=Afa+Afb.*exp(-i*omega*td);
magAf=abs(Af);

subplot(3,2,2*jj)
semilogx(f*1e-3,magAf*1e3,'k','linewidth',1); hold on
xlabel('f (kHz)','fontsize',mn)
ylabel('P (mPaHz^{-1})','fontsize',mn)
axis([.1 100 0 1.2*max(magAf)*1e3])
set(gca,'fontsize',mn)
text(60,1.1*max(magAf)*1e3,laz(jj+3),'fontsize',mn-2)



%Time domain analytic soln for B=0 
t=((0:2^19)-1)/1e7;
deltat=t(2)-t(1);
tlen=length(t);
hh=find(abs(t-to)==min(abs(t-to)));

ta=t(1:hh);
tna=ta/to;

Po1=Am1*rho*c*U*a1*cos(theta)/(2*r*(4*eta1^4+1)); % signal amplitude

s1=2*eta1^2-1;
s2=2*eta1^2+1;

p1a=s1.*cos(pi*tna)+2*eta1.*sin(pi.*tna);
p2a=(s2.*sin(pi.*eta1.*tna)+s1.*cos(pi.*eta1.*tna)).*exp(-pi.*eta1.*tna);
pa=p1a-p2a;

tb=t(hh+1:tlen);
tnb=tb/to;

tz=pi*eta1*(tnb-1);
p1b=(-s1*cos(tz)-s2*sin(tz)).*exp(-tz);
p2b=(s2*sin(pi*eta1*tnb)+s1*cos(pi*eta1*tnb)).*exp(-pi*eta1*tnb);
pb=p1b-p2b;

p1=Po1*[pa pb]; 


Po1=Am2*rho*c*U*a2*cos(theta)/(2*r*(4*eta2^4+1)); % signal amplitude

s1=2*eta2^2-1;
s2=2*eta2^2+1;

p1a=s1.*cos(pi*tna)+2*eta2.*sin(pi.*tna);
p2a=(s2.*sin(pi.*eta2.*tna)+s1.*cos(pi.*eta2.*tna)).*exp(-pi.*eta2.*tna);
pa=p1a-p2a;

tb=t(hh+1:tlen);
tnb=tb/to;

tz=pi*eta2*(tnb-1);
p1b=(-s1*cos(tz)-s2*sin(tz)).*exp(-tz);
p2b=(s2*sin(pi*eta2*tnb)+s1*cos(pi*eta2*tnb)).*exp(-pi*eta2*tnb);
pb=p1b-p2b;

p2=Po1*[pa pb]; 


hh=find(abs(t-td)==min(abs(t-td)));
pp=zeros(1,tlen);
pp(hh+1:tlen)=p2(1:(tlen-hh));, 
p=p1+pp;

subplot(3,2,1+2*(jj-1))
plot(t*1e3,p*1e-3,'k','linewidth',1), hold on
%plot(t*1e3,p10*1e-3,':k','linewidth',1.0)
xlabel('t (ms)','fontsize',mn)
ylabel('P (kPa)','fontsize',mn)
axis([0 0.6 -1.4*max(p)*1e-3 1.4*max(p)*1e-3])
set(gca,'fontsize',mn)
text(0.25,1.2*max(p)*1e-3,['a1=' num2str(a1) ' (m)'],'fontsize',mn-4) 
text(0.25,1.0*max(p)*1e-3,['a2=' num2str(a2) ' (m)'],'fontsize',mn-4) 
text(0.55,1.2*max(p)*1e-3,laz(jj),'fontsize',mn-2)


%Numerical comparisions


subplot(3,2,1+2*(jj-1))
%Time domain by ifft from full freq soln for any impacting particles
delt=1/(2*max(f));
tlen=length(Af);

%obtained by comparing FFT of a time series of an impact with 
%Frequency domain analytical full soln for any impacting particles
xz=[Af(1:tlen) conj(Af(tlen:-1:2))]; 

pt=ifft(xz)/delt;
tt=[0:length(pt)-1]*delt;
plot(tt*1e3,pt*1e-3,'--g')



subplot(3,2,2*jj)
%numerical from power density spectrum   Pfft=2*fft^2/T;
[pespec,fs]=periodogram(p,[],[],1/deltat);
plot(fs*1e-3,sqrt(max(t)*(pespec)/2)*1e3,'--g')

end

