% Continuous wave from fokused line-transducer
%24.01.01  Hans Torp
%Modified to include apodisation 24.01.01  Hans Torp

f=2e6;%frekvens
c=1540;%lydhastighet
w=2*pi*f;%
k=w/c;%bølgetall
lambda=c/f%bølgelengde
%Transducer geometri
a=22e-3;%aperture
R=40e-3;%radius of curvature
dx0=0.7*lambda%avstand mellom punktkilder på transducer


xmax=25e-3;% +/- størrelse i xretning
zmax=50e-3;%størrelse i zretning
Nx=300;%antall punkter i xretning
Nz=400;%antall punkter i zretning

xaxis=linspace(-xmax,xmax,Nx);
zaxis=linspace(0,zmax,Nz)';
x=ones(Nz,1)*xaxis;
z=zaxis*ones(1,Nx);

p0=0*x;%Akustisk felt ved tidspunkt t=0
figure(1);
set(gcf,'Doublebuffer','on');%prevent flickering
image(xaxis,zaxis,real(p0));
colormap(gray);axis('image');
dt=1/f/10;T=1/f;%for animation

wAp=boxcar(round(a/dx0));%apodisation 
%wAp=hamming(round(a/dx0));%apodisation 
wAp=wAp(end/2:end);%cut away left half

animate=1;

p0=0*x;%setter akustisk felt til 0
n=1;
for x0=dx0/2:dx0:a/2,
   z0=R-sqrt(R^2-x0^2);
   r=sqrt((x-x0).^2+(z-z0).^2);%positiv x0
  if x0>a/10, p0=p0+wAp(n)*exp(i*k*r)./r;end;
   r=sqrt((x+x0).^2+(z-z0).^2);%negativ x0
   %p0=p0+wAp(n)*exp(i*k*r)./r;
   n=n+1;
   crossover=(2*x0)^2/lambda;
   if animate,%animation of wave propagation
       T=1/f/4;
       for t=0:dt:T,
           p=p0*exp(-i*w*t);
           image(xaxis,zaxis,0.05*real(p));
           axis('image');pause(0.001);
       end;
   end;
end;

image(xaxis,zaxis,0.09*real(p0));
colormap(gray);axis('image');
if animate,%animation of wave propagation
    dt=1/f/10;T=8/f;
    for t=0:dt:T,
        p=p0*exp(-i*w*t);
        image(xaxis,zaxis,0.05*real(p));
        axis('image');pause(0.001);
    end;
end;
return;

%stråleprofil
% Prøv å sammenlign med de analytiske uttrykkene i N. Wright's kompendium!
z0=1;
pmax=max(max(abs(p0(100:end,:))));
while z0>0,
    disp('Stråleprofil: Klikk i bildet');
    figure(1);
    [x0,z0]=ginput(1);
    if z0>0,
        dz=zaxis(2)-zaxis(1);
        z0ind=1+round(z0/dz);
        dBprofil=20*log10(abs(p0(z0ind,:)/pmax));
        figure(2);plot(xaxis,dBprofil);axis([-xmax,xmax,-50,0]);
        grid;
    end;
end;

