% hornrefl - inverse solution for horn diameter from reflectance
%
function hornrefl
global sr fmax xmax D0 T crt A_lim B_lim r_lim wd
nfrq=2^16;	% number of frequencies
nhrn=3;		% number of horn shapes
nfig=5;		% number of figures
sr=1e6;		% sampling rate
fmax=3e4;	% maximum frequency plotted
xmax=12;	% maximum x
x1=10;		% x at diameter match
D0=1;       % diameter at x=0
T=20; 		% air temperature (deg C)
crt=0;      % compare time-domain reflectance
wd=0;       % write data ?
%
c=34723;
m=((2*sqrt(x1/pi)/D0)-1)/x1;
r_lim=[-1 0.1]*(m*c/1000);
A_lim=[D0-0.1 4.1];
B_lim=[0 2]*m;
%
for k=1:nfig
   figure(k);clf
end
for k=1:nhrn
   [A,B,x,Yo,Yr,f,r1,t,color]=horn_admit(1,k,nfrq,x1);
   [R,f]=horn_refl_fd(2,k,color,Yo,Yr,f);
   [r,t]=horn_refl_td(3,k,color,R,f,r1,t);
   rate=max(f)*2; nt=(length(f)-1)*2; t=t+nt/rate/2;
   [A,B]=horn_inverse(4,k,color,A,B,x,r,t,Yo);
end
return

% horn admittance
function [A,B,x,Yo,Yr,f,r,t,color]=horn_admit(fig_num,horn_type,nfrq,x1)
global sr xmax D0 T A_lim wd fmax horn
r0=D0/2;
A0=pi*r0^2;
A1=x1;
r1=sqrt(A1/pi);
rr=r1/r0;
rho = 1.1769e-3 * (1 - 0.00335 * (T - 26.85));
c   = 3.4723e4  * (1 + 0.00166 * (T - 26.85));
dx=c/sr;
nx=ceil(xmax/dx);
x=(0:nx)'*dx;
%
switch(horn_type)
case 1				% exponential horn
   m=log(rr)/x1;
   r=r0*exp(m*(x));
   B=m*ones(size(x));
   color='b';
   horn='exponential';
case 2				% conical horn
   m=(rr-1)/x1;
   r=r0*(1+m*x);
   B=m./(1+m*x);
   color='g';
   horn='conical';
case 3				% parabolic horn
   m=(rr^2-1)/x1;
   r=r0*sqrt(1+m*x);
   B=(m/2)./(1+m*x);
   color='r';
   horn='parabolic';
end
A=pi*r.^2;
Yo=A(1)/(rho*c);
mc=m*c;
nt=9000;
t=(0:(nt-1))'/sr;
t(1)=1e-12; % avoid division by zero
f=(sr/2)*(0:nfrq)'/nfrq;
f(1)=1e-12; % avoid division by zero
w=2*pi*f;
s=i*w;
switch(horn_type)
case 1				% exponential horn
   Yr=Yo*(sqrt(s.^2+mc^2)+mc)./s;
   r=-besselj(1,mc*t)./t;
case 2				% conical horn
   Yr=Yo*(1+mc./s);
   r=-(mc/2)*exp(-mc*t/2);
case 3				% parabolic horn
   a=-w/mc;
   Yr=i*Yo*besselh(1,a)./besselh(0,a);
   r=-(mc/4)*exp(-mc*t./(2+mc*t./(4+mc*t/16)));
end

%
if (1)
   yo=real(Yr(end)/Yo);
   ys=real(s(end)*(Yr(end)/Yo-1)/(m*c));
   fprintf('r0=%4.2f m=%5.3f mc=%5.0f ',r0,m,mc);
   fprintf('surge/Yo=%5.3f step/mc=%5.3f horn=%s\n',yo,ys,horn);
end
%
figure(fig_num)
D=2*sqrt(A/pi);
plot(x,D,color);
axis([0 xmax A_lim]);
ylabel('diameter (cm)')
xlabel('axial distance (cm)')
hold on
if (1)
   xt=9;
   yt=1+horn_type*0.2;
   plot(xt-[1.2 0.2],[yt yt],color)
   text(xt,yt,horn)
end
drawnow
if (wd)
    write_data(sprintf('hornrefl1%s.txt',color),[x D]);
end
return

% horn reflectance in frequency domain
function [R,f]=horn_refl_fd(fig_num,horn_type,color,Yo,Yr,f)
global fmax wd horn
f_lim=[fmax/1000 fmax];
R=(Yo-Yr)./(Yo+Yr);
if (max(abs(R))<1e-4) return; end;
db=20*log10(abs(R));
ph=unwrap(angle(R))/(2*pi);
figure(fig_num)
subplot(3,1,1)
semilogx(f,db,color)
axis([f_lim -50 10])
title('frequency-domain reflectance')
ylabel('magnitude (dB)')
hold on
subplot(3,1,2)
semilogx(f,ph,color)
axis([f_lim 0.20 0.55])
ylabel('phase (cyc)')
hold on
subplot(3,1,3)
gd=-diff(ph)./diff(f);
ff=(f(1:(end-1))+f(2:end))/2;
semilogx(ff,gd*1000,color)
axis([f_lim -0.1 1.1])
xlabel('frequency (Hz)')
ylabel('delay (ms)')
hold on
if (1)
   xt=10^4;
   yt=0.1+horn_type*0.2;
   plot(xt-[5*10^3 10^3],[yt yt],color)
   text(xt,yt,horn)
end
drawnow
if (wd)
    % magnitude
    [f1, db1] = plt_data(f, db, f_lim(1), f_lim(2));
    filename = sprintf('hornrefl2m%s',color);
    write_data([filename,'.txt'],[f1 db1]);
    % phase
    [f1, ph1] = plt_data(f, ph, f_lim(1), f_lim(2));
    filename = sprintf('hornrefl2p%s',color);
    write_data([filename,'.txt'],[f1 ph1]);
    % group delay
    [ff1, gd1] = plt_data(ff, gd, f_lim(1), f_lim(2));
    filename = sprintf('hornrefl2gd%s',color);
    write_data([filename,'.txt'],[ff1 gd1*1000]);
end
return

% horn reflectance in time domain
function [r,t]=horn_refl_td(fig_num,horn_type,color,R,f,r,t);
global sr D0 r_lim crt wd horn
figure(fig_num)
if (crt)
   line=strcat(color,'--');
   plot(t*1000,r/1000,line)
   hold on
end
t0=t; r0=r;
%
rate=max(f)*2;
nt=(length(f)-1)*2;
t=(0:(nt-1))'/rate;
t_lim = [-0.1 1.1];
r=ffs(R);
if (max(abs(r))<1e-4), return; end
r=fftshift(r);
t=t-nt/rate/2;
plot(t*1000,r*sr/1000,color)
t1=t; r1=r;
r=fftshift(r);
axis([t_lim r_lim])
xlabel('time (ms)')
ylabel('(kHz)')
title('time-domain reflectance / {\Delta}t')
hold on
if (1)
   xt=0.8;
   yt=-8+horn_type*0.5;
   plot(xt-[0.2 0.05],[yt yt],color)
   text(xt,yt,horn)
end
drawnow
if (wd)
    % approximation
    [t3, r3] = plt_data(t0, r0, t_lim(1)/1000, t_lim(2)/1000);
    filename = sprintf('hornrefl3ap%s',color);
    write_data([filename,'.txt'],[t3*1000 r3/1000]);
    % inverse fft
    [t2, r2] = plt_data(t1, r1, t_lim(1)/1000, t_lim(2)/1000);
    filename = sprintf('hornrefl3fft%s',color);
    write_data([filename,'.txt'],[t2*1000 r2*sr/1000]);
end
return

% horn inverse solution
function [A,B]=horn_inverse(fig_num,horn_type,color,A,B,x,r,t,Yo)
global sr xmax T A_lim B_lim wd horn
figure(fig_num)
line=strcat(color,'--');
plot(x,2*sqrt(A/pi),line);
axis([0 xmax A_lim]);
hold on
figure(fig_num+1)
plot(x,B,line);
hold on

x1 = x; A1 = A; B1 = B;
%
figure(fig_num)
zo=1/Yo;
[A,B,x]=refl_inverse(r,zo,sr,xmax,T);
plot(x,2*sqrt(A/pi),color);
axis([0 xmax A_lim]);
ylabel('diameter (cm)')
xlabel('axial distance (cm)')
hold on
if (1)
   xt=9;
   yt=1+horn_type*0.2;
   plot(xt-[1.2 0.2],[yt yt],color)
   text(xt,yt,horn)
end

figure(fig_num+1)
plot(x,B,color);
axis([-1 xmax B_lim]);
title('inertance / \rho')
ylabel('(cm^{-1})')
xlabel('axial distance (cm)')
hold on
drawnow


if (wd)
    % D(x) - true value
    [x2, A2] = plt_data(x1, A1, 0, xmax);
    filename = sprintf('hornrefl4Dtv%s',color);
    write_data([filename,'.txt'],[x2 2*sqrt(A2/pi)]);
    % B(x) - true value
    [x2, B2] = plt_data(x1, B1, 0, xmax);
    filename = sprintf('hornrefl5Btv%s',color);
    write_data([filename,'.txt'],[x2 B2]);
    % D(x) - inverse sloution
    [x3, A3] = plt_data(x', A', 0, xmax);
    filename = sprintf('hornrefl4Dis%s',color);
    write_data([filename,'.txt'],[x3 2*sqrt(A3/pi)]);
    % B(x) - inverse sloution
    [x3, B3] = plt_data(x', B', 0, xmax);
    filename = sprintf('hornrefl5Bis%s',color);
    write_data([filename,'.txt'],[x3 B3]);
end
return

% inverse solution for diameter from reflectance
function [A,B,x]=refl_inverse(r,zo,rate,xmax,T);
rho = 1.1769e-3 * (1 - 0.00335 * (T - 26.85));
c   = 3.4723e4  * (1 + 0.00166 * (T - 26.85));
dt=1/rate;
dx=c*dt;
Ao=rho*c/zo;
m=ceil(xmax/dx);
if (length(r)<(2*m))
   m=floor(length(r)/2);
   xmax=m*dx;
   fprintf('reduced xmax=%5.2f\n',xmax);
end
x=(0:(m-1))*dx;
epmx=1/dx;						% maximum epsilon
etmx=1/dx;						% maximum eta
pfmx=1.01; 						% maximum power flow
pfmn=-0.01;						% minimum power flow
%
n=2*m;
pm=r(1:n);
pp=zeros(n,1);
pp(1)=1;
pm(1)=0;
A=zeros(1,m);
B=zeros(1,m);
A(1)=Ao;
for k=1:(m-1)
   p=pp+pm;						% pressure
   q=cumsum(p)*dt;              % integrated pressure
   u=pp-pm+c*B(k)*q;			% velocity / Yo
   a=A(k)/Ao;					% relative area
   pf=p(k)*conj(u(k))*a;        % power flow
   if ((pf>pfmx)|(pf<pfmn)) break; end	% check #1
   j1=k;
   j2=k+1;
   bk=B(k)/2;
   v1=bk*(p(j1)+u(j1))-pm(j1)/dx;
   v2=bk*(p(j2)+u(j2))-pm(j2)/dx;
   dd=u(j1)*q(j2)-q(j1)*u(j2);
   if (dd==0), break; end					% check #2
   ep=(q(j2)*v1-q(j1)*v2)/dd;
   et=(u(j1)*v2-u(j2)*v1)/dd;
   if ((ep>epmx)), break; end				% check #3
   if ((et>etmx)&(bk<-1)), break; end        % check #4
   npp=zeros(n,1);
   npm=zeros(n,1);
   for j=(k+1):(n-k)
      jp=j+1;
      jm=j-1;
      vp=bk*p(jp)-((ep-bk)*u(jp)+et*q(jp));
      vm=bk*p(jm)+((ep-bk)*u(jm)+et*q(jm));
      npp(j)=pp(jm)-vm*dx;
      npm(j)=pm(jp)-vp*dx;
   end
   pp=npp;
   pm=npm;
   A(k+1)=A(k)*exp(2*ep*dx);
   B(k+1)=B(k)+2*et*dt;
end
if (k<(m-1))
   fprintf('x=%5.2f A=%6.3g B=%6.3g ',x(k),A(k),B(k));
   fprintf('ep=%6.3g et=%6.3g pf=%6.3g\n',ep,et,pf);
end;
return

% fast Fourier synthesize real signal
function h=ffs(H)
m=length(H);
n=2*(m-1);
H(1,:)=real(H(1,:));
H(m,:)=real(H(m,:));
H((m+1):n,:)=conj(H((m-1):-1:2,:));
h=real(ifft(H));
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function write_data(fn,data)
[nr,nc] = size(data);
fp=fopen(fn,'wt');
fprintf(fp,'; %s\n', fn);
for i=1:nr
    for j=1:nc;
        fprintf(fp,' %14.5g',data(i,j));
    end
    fprintf(fp,'\n');
end
fclose(fp);


function [x1, y1]=plt_data(x, y, xmin, xmax)
j1 = getSampleNumber(x,xmin);
j2 = getSampleNumber(x,xmax);
x1 = x(j1:j2);
y1 = y(j1:j2);
return


function sampleNum = getSampleNumber(fun, value)
for i = 1:length(fun)-1
    if fun(i) > value
        sampleNum = i-1;
        if sampleNum > 0
            return
        else sampleNum = 1;
            return
        end
    else
        sampleNum = length(fun);
    end
end