  echo off
j=sqrt(-1);
f = logspace(log10(10),3,50);
n = max(size(f));
w = 2*pi*f;
c=344; r0=1.18;
%------------------ start user Input
%input voltage
e=2.83;
%box volume
Vol = 0.39e-3;
%distance 1m
r=1;
%driver parameters (i.e. FRS 8)
Re = 7.2;
Le= 0.85e-3;
fs=120;
Bl = 3.2;
Sm = 31e-4;
Vas = 0.91e-3;
Qt = 1.04;
Qe = 1.32;
Qm = 4.85;
%----------------------- end user input
Cm = Vas/(r0*c^2*Sm^2);
Mm=1/((2*pi*fs)^2*Cm);
Rm = 1/((2*pi*fs)*Qm*Cm);
a=sqrt(Sm/pi);
ws=2*pi*fs;
pout = ones(1,n);
pin = ones(1,n);
pnf = ones(1,n);
xuk=ones(1,n);
zul=ones(1,n);
%radiation impedance constants
ra1 = 0.1404*r0*c/a^2;
ra2 = r0*c/(pi*a^2);
ma1 = 8*r0/(3*pi^2*a);
ca1 = 5.94*a^3/(r0*c^2);
% box constants
Cb=Vol/(r0*c^2); %box compliance
Mao=0.85*a*r0/Sm; %membrane air mass loading
for k=1:n
 ww = w(k); jw = j*ww;
 % radiation impedance
 z1 = jw*ma1; z2 =1/(jw*ca1);
 z3=z2*ra1/(ra1+z2); z3= z3+ra2;
 Zar=1/(1/z1+1/z3);

 %box analogous circuit
 Ze= Re + jw*Le;
 Zm= Rm + jw*Mm +1/(jw*Cm);
 A11 = Ze*Sm/Bl;
 A12 = Ze*Zm/(Sm*Bl) + Bl/Sm;
 A21 = Sm/Bl;
 A22 = Zm/(Sm*Bl);
 Zak =Zar + jw*Mao + 1/(jw*Cb);
 zul(k)= (A11*Zak+A12)/(A21*Zak+A22);
%volume velocity
 U1 = e/(A11*Zak+A12);
 %displacement
 xuk(k) = (U1/Sm)/jw;
 %presure at distance r (infinite baffle mounted)
 pout(k) = jw*r0*U1/(2*pi*r);
%pressure inside the box
 pin(k) = U1/(jw*Cb);
 %near field pressure
 kk = ww*a/2/c;
 pnf(k)= jw*r0*(U1/(a*pi))*sin(kk)/kk;

end
%level in far filed (db/1m)
pdb= 20*log10(abs(pout)/2e-5);
%level in the box
pudb= 20*log10(abs(pin)/2e-5);
%level in the near field
pnfdb= 20*log10(abs(pnf)/2e-5);
semilogx(f,pudb),title('SPL (dB) inside the box'),grid,pause
semilogx(f,pnfdb),title('SPL (dB) in near field'),grid,pause
subplot(221);
echo off;
semilogx(f,pdb),title('SPL (dB/2pi)'),grid;
xlabel('f (Hz)');
subplot(222);
semilogx(f,abs(xuk)*1000),title('Displacement (mm) '),grid;
xlabel('f (Hz)');
subplot(223);
semilogx(f,angle(pout)*180/pi),title('Phase (deg) '),grid;
xlabel('f (Hz)');
subplot(224);
semilogx(f,abs(zul)),title('Impedance (ohm)'),grid;
xlabel('f (Hz)'), pause;
subplot;  
