%n:\mat5work\sizescatter\abs_inversion\inversion_for_M_for_paper.m   15\10\13

%This software is in support of the paper
% 'An overview on the use of backscattered sound for measuring suspended particle size and concentration profiles
% in non-cohesive inorganic sediment transport studies'.
% by Peter D. Thorne and David Hurther in Continental Shelf Research

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


%1 GENERATES THE SUSPENSION FIELD. 
%2 CALCULATE AVERAGE SCATTERING CHARACTERISTICS 
%3 COMPUTES THE PRESSURE FIELD. 
%4 INVERSION TO OBTAIN CONCENTRATION

close all
clear all

%  ********1. CONSTRUCT CONCENTRATION AND PARTICLE SIZE SUSPENSION FEILD***********

% taken from n:\mat5work\sizescatter\abs_inversion\sus_c_a1.m


%range
deltar=0.005; 
r=deltar:deltar:1.05; 
%height
bed=1.05+deltar; z=bed-r; 
lr=length(r);


%underlying structure with temporal variability
nn=110;
xz=1:nn;
y=abs(0.8*(sin(2*pi*xz/20)+4*cos(2*pi*xz/15)+3*sin(2*pi*xz/40))+2*cos(2*pi*xz/17));
y=y/5;

rand('seed',0), %keeps field same from run to run

for vz=1:nn
%vertical concentration profile Rouse which near the bed becomes a power law
bp=1.0; %ws/(0.4*ustar);
za=deltar;
modelconc=(z/za).^(-bp);% relative profile in eqn 12a in the paper
modelconcb=(y(vz).^2)*modelconc;
%add reandom component
modelconc=modelconcb+0.5*rand(1,210).*modelconc;

%particle size
%power law
modelsize=(z/za).^(-.2); % relative profile in eqn 12b in the paper
modelsizeb=modelsize+y(vz)*modelsize;
%add reandom component
modelsize=modelsizeb+0.5*modelsize.*rand(1,210);

cnc(vz,:)=modelconc;
a_cnc(vz,:)=modelsize; 


end
clear modelconcb modelsizeb

%sets up absolute concentration and size, these are not reference values.
c0=13; a0=48.7e-6; 

%Smoothing
conc=medfilt2(cnc,[2 2]); conc=c0*conc(5:104,10:209);
a_size=medfilt2(a_cnc,[2 2]); a_size=a0*a_size(5:104,10:209);
z2=z(11:210);



%Reference values
MM1=min(conc(:,200)); MM2=mean(conc(:,200)); MM3=max(conc(:,200));
AA1=1e6*min(a_size(:,200)); AA2=1e6*mean(a_size(:,200)); AA3=1e6*max(a_size(:,200)); 
%from n(a) particle number radius convert to particle mass radius m(a) using (1/0.69)
D50=2*(1/.69)*1e6*max(a_size(:,200)); 


%fontsize
nz=15;

%SIZE radius

% a_size particle mean radius field contourf
figure(1), orient tall
axes('position',[.1 .65 .525 .23])
contourf([1:100],[0.005:0.005:1],1e6*(flipud(a_size')),10),hold on
text(105,1.06,'a_c (\mum)','fontsize',nz-2)
set(gca,'xticklabel',[])
set(gca,'fontsize',nz)
ylabel('z (m)','fontsize',nz)
text(95,0.9,'a','fontsize',nz,'color',[1 1 1])
axis([1 100 0 1])
title(['a_{rmin}=' num2str(AA1), ' ' 'a_{rmax}=' num2str(AA3), ' ' 'a_{rmean}=' num2str(AA2), ' ' 'd50=' num2str(D50)]) 

% make colour bar
x=(0:12:119);  onz=ones(2,1);
xonz=onz*x;
axes('position',[.68 .65 .03 .23])
contourf(1:2,0:12:119,xonz',10)
set(gca,'xticklabel',[])
set(gca,'fontsize',nz)


% mean a_size profile
axes('position',[.75 .55 .23 .25])
hp=loglog(1e6*mean(a_size),z2,'k'); set(hp,'linewidth',1.3)
axis([10 200 0.005 1])
set(gca,'yticklabel',[])
text(120,0.7,'b','fontsize',nz)
text(8,.98,'1','fontsize',nz)
text(8,.005,'0','fontsize',nz)
text(8,0.015,'z (m)','Rotation',90,'fontsize',nz)
xlabel('(\mum)','fontsize',nz)
set(gca,'fontsize',nz)

% a_size reference at z=0.005m
axes('position',[.1 .53 .525 .1])
hp=plot(1e6*a_size(:,200),'k'); set(hp,'linewidth',1.3)
text(-10,55,'a_r (\mum)','Rotation',90,'fontsize',nz)
text(95,110,'c','fontsize',nz)
set(gca,'fontsize',nz)
set(gca,'xticklabel',[])
axis([0 100 50 120])



% CONCENTRATION

% conc suspensded concentration contourf
axes('position',[.1 .2 .525 .23])
contourf([1:100],[0.005:0.005:1],log10(flipud(conc')),10), hold on
text(102, 1.06,'log(C) (kgm^{-3})','fontsize',nz-2)
set(gca,'xticklabel',[])
set(gca,'fontsize',nz)
ylabel('z (m)','fontsize',nz)
text(95,0.9,'d','fontsize',nz,'color',[1 1 1])
axis([1 100 0 1])
title(['C_{rmin}=' num2str(MM1), ' ' 'C_{rmax}=' num2str(MM3), ' ' 'C_{rmean}='  num2str(MM2)]) 

% make colour bar
x=(-4:.5:1.5);  onz=ones(2,1);
xonz=onz*x;
axes('position',[.68 .2 .03 .23])
contourf(1:2,-4:.5:1.5,xonz',10)
set(gca,'xticklabel',[])
set(gca,'fontsize',nz)


% mean conc profile
axes('position',[.75 .1 .23 .25])
hp=loglog(mean(conc),z2,'k'); set(hp,'linewidth',1.3)
axis([.01 20 0.005 1])
text(5e-3,.98,'1','fontsize',nz)
text(5e-3,.0,'0','fontsize',nz)
text(5e-3,0.015,'z (m)','Rotation',90,'fontsize',nz)
set(gca,'yticklabel',[])
text(5,0.7,'e','fontsize',nz)
set(gca,'fontsize',nz)
xlabel('(kgm^{-3})','fontsize',nz)

% conc reference concentration at z=0.005m
axes('position',[.1 .08 .525 .1])
hp=semilogy(conc(:,200),'k'); set(hp,'linewidth',1.3)
axis([0 100 1 20])
ylabel('C_r (kgm^{-3})','fontsize',nz)
xlabel('Time','fontsize',nz)
text(95,15,'f','fontsize',nz)
set(gca,'fontsize',nz)

save acfield_paper a_size conc

clear all


  %******** 2. CALCILATE AVERAGED FORM FUNCTION AND NORMALISED TOTAL SCATTERING CROSS-SECTION***********


%n:\mat5work\sizescatter\abs_inversion\av_chi_fm_lognorm1.m  

clear all

x=0.001:0.001:200; %x=ka where a is the particle size and k the acousitc wavenumber in the fluid

%Moate B.D. and Thorne P.D. 2012. Interpreting acoustic backscatter from suspended sediments
%of different and mixed mineralogical composition. Continental Shelf Research, 46, 67-82.

%calculate the intrinsic scattering expressions eqn's 3a and 3b in the paper

% Expression for f/sqrt(rho) 
term1=( 1 - (0.25*exp( -1*(((x-1.5)/0.35).^2) )) );
term2=( 1 + (0.6*exp( -1*(((x-2.9)/1.15).^2) )) );
E=term1.*term2;
fe=(E.*(x.^2))./( 42 + (28*(x.^2)) );
% Expression for chi/rho 
chi=( 0.09*(x.^4) )./( 1380 + (560*(x.^2)) + (150*(x.^4)) );


%calculates ensemble scattering characteristics
nz=1000;
dsigma=0.4;, dx=x(2)-x(1);, xx=logspace(-2,1.7,nz);


%Lognormal eqn (9b) in the paper 
for jj=1:nz;
    xo=xx(jj);
sigma=dsigma*xo;
mu=log((xo.^2)./sqrt(xo.^2+sigma.^2));
sigman=sqrt(log((sigma/xo).^2+1));
plognorm=(1./(x*sigman*sqrt(2*pi))).*exp(-((log(x)-mu).^2)/(2*sigman.^2));

axlng(jj)=dx*sum(x.*plognorm);   %eqn 2c in the paper
ax3=dx*sum((x.^3).*plognorm);
axf=dx*sum(((x.*fe).^2).*plognorm);
axchi=dx*sum((x.^2).*chi.*plognorm);
aflng(jj)=sqrt((axlng(jj).*axf)./ax3); %eqn 2a in the paper
achilng(jj)=(axlng(jj).*axchi)./ax3; %eqn 2b in the paper
end

figure(2), orient tall
subplot(2,1,1)
loglog(axlng,aflng,'k'),hold on 
plot(x,fe,'r')
xlabel('x,x_o','fontsize',15), ylabel('f/\rho^{1/2}','fontsize',15)
title(['dsigma=' num2str(dsigma)])
axis([1e-2 1e1 1e-4 1e-1])
set(gca,'fontsize',15)
text(8,1e-2,'a','fontsize',15)
legend('ensemble scattering', 'intrinsic scattering')

subplot(2,1,2)
loglog(axlng,achilng,'k'), hold on
plot(x,chi,'r')
xlabel('x,x_o','fontsize',15), ylabel('\chi/\rho','fontsize',15)
axis([1e-1 1e1 1e-6 1e-3])
set(gca,'fontsize',15)
text(8,1e-4,'b','fontsize',15)


save avchifm_paper axlng achilng aflng

clear all


%            ************3  CALCULATION OF THE BACKSCATTERED FIELD*******


% n:\mat5work\sizescatter\abs_inversion\press2MHzb.m    21-11-11

%calculates bs at 2MHz for the inversion
%mod of press2MHz.m using seperating errors for scattering and calibration


clear all

tic

%n:\mat5work\sizescatter\abs_inversion\sus_c_a1.m
load acfield_paper
%n:\mat5work\sizescatter\abs_inversion\av_chi_fm_lognorm1.m   
load avchifm_paper

%forward problem calculate the backscattered signal

%range
deltar=0.005; 
r=deltar:deltar:1.0; 
%height
bed=1+deltar; z=bed-r; 
lr=length(r);

%THIS IS WHERE UNCERTAINTIES ARE INTRODUCED INTO THE BACKSCATTERED SIGNAL
%theta=1.0  and drand=0.0 no errors
%typical error would be between 0-20%
theta=1.0; %systematic (epsilon the paper)
drand=0.0;  %random   (gammo_o in the paper)


% ************    calculate the backscattered signal     ************

%sound vel
c=1500; 
%frequency
f=2.0*1e6; 
%transducer radius
at=5.0*1e-3; 
%water attanuation
alw=f.*f*31e-15; % alpha_w in the paper at 15C for distilled water
%system constant from measurement given by R in the paper
meankt=0.0175;

zxo=rand(1);
meankt=meankt*theta; %systematic error in the paper
meankt=meankt.*(1+drand*(2*zxo-1)); % random error in the paper

cnc=conc; 
a=a_size; 


% near field correction
%------------------------------------------------
k=2*pi*f/c;

	Rstar=r./((k/2).*at.^2);
	psi = 1./(1 - 1./(1+1.35.*Rstar+(2.5.*Rstar).^3.2))';
   
   
 % attenuation constants for water and sediment
 %-----------------------------------------------

for jj=1:100 
 
x=k*a(jj,:);

%form function 
fm=interp1(axlng,aflng,x);
ao3=a(jj,:);
ks=fm./sqrt(ao3); % sediment backscattering properties K in paper

ks=ks*theta; %systematic error in the paper
zx=rand(1,200);
ks(1,:)=ks(1,:).*(1+drand*(2*zx-1));  % random error in the paper


%attenuation
chi=interp1(axlng,achilng,x);
eta=3*chi./(4*ao3); % attenuation scattering properties greek-Xi in the paper

eta=eta*theta; %systematic error in the paper
eta(1,:)=eta(1,:).*(1+drand*(2*zx-1)); % random error in the paper

%sediment attenuation alpha_s in the paper
als1=cumsum(eta(1,:).*cnc(jj,:))./[1:lr];
att1=alw(1)+als1;


%backscatter signal profiles, eqn A1.4 in the paper
v1(jj,:)=((ks(1,:).*meankt(1).*sqrt(cnc(jj,:)))./(r.*psi(1,:))).*exp(-2*r.*att1);	

%add random component to backscattering signal due to Rayleigh noise 
v1(jj,:)=v1(jj,:).*(1+drand*(2*rand(1,200)-1));

end


save press2MHz_paper v1 r z lr k psi c f meankt alw theta drand

figure(3)
semilogy(r,v1')
title(['press2Mz-paper' '  freq=' num2str(f'/1e6) '  theta=' num2str(theta) '  drand=' num2str(drand)])
axis([0 1.01 1e-4 1e-1])
ylabel('Backscattered signal (volts)' )
xlabel('Range (m)')

clear all



%            ************4  CALCULATION FOR M*******


% taken from n:\mat5work\sizescatter\abs_inversion\press_inv_c1.m  
clear all, 

%n:\mat5work\sizescatter\abs_inversion\av_chi_fm_lognorm1.m   21/11/11
load avchifm_paper   % dsigma=0.4 gives ensemble scattering
% n:\mat5work\sizescatter\abs_inversion\press2MHzb.m    21-11-11
load press2MHz_paper
%n:\mat5work\sizescatter\abs_inversion\sus_c_a1.m
load acfield_paper

M1=zeros(100,lr); 
Z1=zeros(100,lr);

for jj=1:100 

atten=zeros(1,lr); eta=zeros(1,lr);

%select size for the inversion

%inversion check
ao=a_size(jj,:);  %particle size exact 

%using d_50 of the bed a_n=0.69*d_50/2 d_50=320
%ao=(1/0.69)*110e-6*ones(1,lr);  % an

%uses ao of the bed
%ao=110e-6*ones(1,lr);  %ao of the bed aprox max of a_r=max(a_cnc(:,200))

%uses ao of the bed with a profile where size reduces above the bed in the
%mean same way as the suspension field
%ao=110e-6*(z/z(200)).^(-.2);

a=ao;  
x=k.*a; 
    
   
%form function 
fm=interp1(axlng,aflng,x); ks=fm./sqrt(a);

%attenuation
chi=interp1(axlng,achilng,x); eta=3*chi./(4*a); 

%to make atten negligible choose this option
%eta=eta*1e-10; 
%to remove nearfield correction choose this option
%psi=psi.\psi;


%initial mass eqn (4) in the paper with the nearfield correction included
M1a=((v1(jj,:).*r.*psi(1,:).*exp(2*r*(alw(1))))./(ks(1,:)*meankt(1))).^2;

% iterations to account for attenuation equations (6) and (7) in the paper
for n=1:lr
    
    M1(jj,n)=M1a(n);
   
for loop=1:50
Mtemp1=M1(jj,n);
atten(1,n)=sum(eta(1,1:n).*M1(jj,1:n))/n; 
M1(jj,n)=M1a(n).*exp(4*r(n)*atten(1,n));
if abs((M1(jj,n)-Mtemp1)/M1(jj,n))<0.0001, break, end
Z1(jj,n)=loop;
end %loop
 
end  %n

end %jj


figure(4), orient tall, mn=15;
title('ration of M/C')
subplot(3,1,1) 
semilogy(r,M1./conc,'k') 
ylabel('M/C','fontsize',mn)
title(['ratio of the M and C fields' '   ' '\epsilon=' num2str(theta) '  ' '\gamma_o=' num2str(drand)])
axis([0 1 0.1 10])


subplot(3,1,2), 
semilogy(conc(:,200),'k'), hold on
semilogy(conc(:,20),'k')
plot(M1(:,200),'.r')
plot(M1(:,20),'.r')
xlabel('Time'), ylabel('(kgm^{-3})')
title('compare M with C at z=0.1 & 1 m')
axis([0 100 .001 200])


subplot(3,1,3), 
semilogy(r,mean(conc),'k'), hold on
plot(r,mean(M1),'.r')
xlabel( 'r (m)'), ylabel('C  (kgm^{-3})')
title('mean concentration profiles for M  & C')
axis([0 1 .001 20])
legend('C','M')

save tempdata1 M1


