%attnorm.m
a=.5e-3;		%particle radius
mu=1e-3;		%liquid dynamic viscosity
rhol=1e3;		%liquid density
kinvisc=mu/rhol;
rhos=2.5e3;		%solid density
Kl=2.25e9;		%bulk modulus of water
Ks=2.e10;		%bulk modulus of glass
%f=1e5;
f=0.5e6;
%f=4.77e5;
w=2*pi*f;
i=sqrt(-1);
delta=sqrt(2*kinvisc/w);
AA=(9*mu/(2*a*a))*(1+a/delta);
DD=rhos/rhol;
nu=.02:.02:.62;
dnu=(ones(nu)-nu);
C11=0.5*(ones(nu)+3.32*nu);
%C00=0.5;
%Cnu=0.5*dnu;
BB1=rhol*(C11+(9*delta/(4*a))*ones(nu));
%BB0=rhol*(C00+(9*delta/(4*a))*ones(nu));
%BB=rhol*(Cnu+(9*delta/(4*a))*ones(nu));
rhostar=dnu+DD*nu;
rhohat=nu+DD*dnu;
const1=((1/Ks)*nu)+((1/Kl)*dnu);
%c0=ones(nu)./(rhol*(const1.*rhostar));
%c0=sqrt(c0);
%k2num=((AA*ones(nu)+(i*w)*BB).*(rhol*rhostar))+(i*w*rhos*rhol)*dnu;
%k2den=AA+((i*w*rhol)*(dnu.*rhohat))+i*w*BB;
%k2=((w^2*const1).*k2num)./k2den;
%k2num0=((AA*ones(nu)+(i*w)*BB0).*(rhol*rhostar))+(i*w*rhos*rhol)*dnu;
%k2den0=AA+((i*w*rhol)*(dnu.*rhohat))+i*w*BB0;
%k20=((w^2*const1).*k2num0)./k2den0;
k2num1=((AA*ones(nu)+(i*w)*BB1).*(rhol*rhostar))+(i*w*rhos*rhol)*dnu;
k2den1=AA+((i*w*rhol)*(dnu.*rhohat))+i*w*BB1;
k21=((w^2*const1).*k2num1)./k2den1;
%
%c=w*(ones(nu)./(real(sqrt(k2))));
%alphaa=abs(a*imag(sqrt(k2)));
%c2=w*(ones(nu)./(real(sqrt(k20))));
%alphaa2=abs(a*imag(sqrt(k20)));
c1=w*(ones(nu)./(real(sqrt(k21))));
alphaa1=abs(a*imag(sqrt(k21)));
alphaa1=alphaa1./(sqrt(mu*w)*ones(alphaa1));
clear Kl Ks rhol rhos kinvisc f delta AA DD dnu Cnu BB rhostar rhohat
clear const1 k2 k2num k2den a w i
%plot(nu,c,'-r',nu,c0,'-g',nu,c2,'-b',nu,c1,'--g')
%pause
%
%ATKINSON ATTENUATION DATA
%
xmin=0;
delx=0.65;
ymin=0.;
dely=0.005;
xwidth=6.41;		%inches
ywidth=6.36;		%inches
xtriang=[.43
.816
1.323
1.6
1.925
2.18
2.425
2.716
3.022
3.4
3.834
4.507
4.941
5.47
5.822
6.129
6.129];
ytriang=[.397
.507
.907
.9925
1.067
1.234
1.329
1.155
1.227
1.236
1.411
1.4
1.376
1.358
1.17
.86
1.29];
xsquare=[.43
.816
1.323
1.6
1.925
2.18
2.425
2.716
3.022
3.4
3.834
4.507
4.941
5.47
5.822
6.129
6.129];
ysquare=[.358
.502
.665
.750
.946
1.144
1.14
1.375
1.146
1.143
1.242
1.347
1.594
1.574
1.323
1.313
1.09];
xcirc=[.43
.816
1.323
1.6
1.925
2.18
2.425
2.716
3.022
3.4
3.834
4.507
4.941
5.47
5.822
6.129
6.129];
ycirc=[1.171
1.58
2.07
2.295
2.294
2.4
2.52
2.4
2.386
2.492
2.6
2.67
2.42
2.266
1.77
1.75
1.92];
xstar=[.43
.816
1.323
1.6
1.925
2.18
2.425
2.716
3.022
3.4
3.834
4.507
4.941
5.47
5.822
6.129];
ystar=[1.728
2.325
2.942
3.343
3.594
4.031
4.222
4.021
4.176
4.273
4.243
4.292
4.228
3.905
3.08
3.02];
xbsqu=[.43
.816
1.323
1.6
1.925
2.18
2.425
2.716
3.022
3.4
3.834
4.507
4.941
5.47
5.822
6.129];
ybsqu=[2.409
3.203
4.136
4.38
4.78
4.99
5.54
5.54
5.65
5.81
5.76
5.65
5.53
5.35
4.25
4.84];
xsquare=xmin+(delx/xwidth)*xsquare; 
ysquare=(1/sqrt(2*pi*mu*.2e6))*(ymin*ones(ysquare)+(dely/ywidth)*ysquare);
xtriang=xmin+(delx/xwidth)*xtriang; 
ytriang=(1/sqrt(2*pi*.25e6*mu))*(ymin*ones(ytriang)+(dely/ywidth)*ytriang);
xcirc=xmin+(delx/xwidth)*xcirc; 
ycirc=(1/sqrt(2*pi*.3e6*mu))*(ymin*ones(ycirc)+(dely/ywidth)*ycirc);
xstar=xmin+(delx/xwidth)*xstar; 
ystar=(1/sqrt(2*pi*.35e6*mu))*(ymin*ones(ystar)+(dely/ywidth)*ystar);
xbsqu=xmin+(delx/xwidth)*xbsqu; 
ybsqu=(1/sqrt(2*pi*.4e6*mu))*(ymin*ones(ybsqu)+(dely/ywidth)*ybsqu);
clear xmin delx  ymin dely xwidth ywidth
%
%URICK DATA
%
xmin=0;
delx=0.65;
ymin=0.;
dely=220;
xwidth=6.41;		%inches
ywidth=6.37;		%inches
xurick=[.276
.769
.783
1.054
1.737
1.833
2.406
2.667
2.985
3.854];
yurick=[1.88
3.54
3.744
5.67
6.017
5.85
5.932
4.56
5.163
3.621];
xurick=xmin+(delx/xwidth)*xurick; 
yurick=(1/sqrt(2*pi*mu*1e6))*(ymin*ones(yurick)+(dely/ywidth)*yurick);
a=0.5e-6;
const=10*log10(exp(1));
yurick=(a/const)*yurick;
clear xmin delx ymin dely xwidth ywidth a const
%
%HAMPTON DATA
%
xmin=0;
delx=0.65;
ymin=0;
dely=26;
xwidth=6.4;		%inches
ywidth=6.37;		%inches
xhampton=[.1345
.354
.463
.659
.643
1.115
1.228
1.494
1.506
1.6885
2.022
2.14
2.515
2.615
3.016
3.016
3.7];
yhampton=[.798
.901
1.625
1.715
3.182
3.2
4.573
4.07
3.18
4.54
3.69
2.55
3.19
2.95
3.3
2.56
.976];
xhampton=xmin+(delx/xwidth)*xhampton; 
yhampton=(1/sqrt(2*pi*mu*.1e6))*(ymin*ones(yhampton)+(dely/ywidth)*yhampton);
a=1e-6;
const=10*log10(exp(1));
yhampton=(a/const)*yhampton;
clear xmin delx ymin dely xwidth ywidth a const
semilogy(xsquare,ysquare,'+r',xtriang,ytriang,'xg',xcirc,ycirc,'ob',xstar,ystar,'*',xbsqu,ybsqu,'+',nu,alphaa1,'-g',xhampton,yhampton,'or',xurick,yurick,'*b')
%pause
%semilogy(nu,abs((real(k2))./(imag(k2))),'r',nu,abs((real(sqrt(k2)))./(imag(sqrt(k2)))),'g')