function standingwave(f,mywait,decayk,tubetype,truncinput)
% STANDINGWAVE Demo of development of standing waves in acoustic tubes
% function standingwave(f,mywait,decayk,tubetype,truncinput)
% standingwave: Version 17.11.2014
%   Syntax
%       tubetype: 1 = closed at glottis, open at lips; 2=closed at both ends; 3 = open at both ends
%       decayk: number of passes through vocal tract required for initial
%           amplitude to decay to 1/e. So larger numbers mean less decay, and
%           final amplitude of standing wave will be larger.
%           A decay value of about 10 makes it possible to distinguish input and
%           reflected waves reasonably quickly, whereas decay of 100 mean they
%           stay virtually identical for a long time.
%       f: input frequency, in vocal tract units, i.e use values like 1/4,
%           1/2, 3/4 etc.
%       mywait: delay (in s) between each calculation. Use 0 to quickly see
%           a fully-developed standing wave. Use e.g about 0.1 to follow
%           individual reflections
%       truncinput:
%           If true, set input to zero after one pass from glottis to lips.
%           To see a fully developed standing wave needs to be false (0).
%           But it is easier to see the reflections at each end if input is
%           truncated.
%       Other settings (hardwired):
%           autoscale: Determines whether scaling is adjusted as
%           oscillation grows.
%           Input (i.e first pass through the tube) is shown in red.
%           Result (i.e sum of input and all reflections) is shown in
%           black.
%           Reflections from right to left shown in green.
%           Reflections from left to right shown in mauve.
%           The leading edge of the pressure disturbance is shown be a blue
%           circle.
%
%   See Also
%       A good animation of the physics behind the polarity switch of the reflection at the
%       open end of a tube can be found here:
%       www.phys.unsw.edu.au/jw/flutes.v.clarinets.html

vtp=60;
vtp1=vtp+1;
vtpm1=vtp-1;
nr=100;

vtlen=17;
c=34000;        %in cm/s
timeinc=(vtlen/c)/vtp;
timeincms=timeinc*1000;

autoscale=1;
vstate=ones(nr,vtp1)*NaN;

%decayk=2;   %decay constant, in vocal tract units

%decay constant in effect determines the final amplitude reached
% (assuming input not truncated)
myylim=decayk+1;        
mydecay=exp(-(1:(nr*vtp1))/(decayk*vtp));
decaylim=0.01;
vd=find(mydecay<decaylim);
mydecay(vd)=NaN;

mydecay=reshape(mydecay,[vtp1 nr]);
mydecay=mydecay';
%flip every other line

mydecay(2:2:end,:)=fliplr(mydecay(2:2:end,:));

bign=vtp*nr;
endhard=0;
starthard=1;
if tubetype==2
    endhard=1;
starthard=1;
end;

if tubetype==3
    endhard=0;
starthard=0;
end;

xvec=(0:(vtp1-1))*(vtlen/vtp);

hl=plot(repmat(xvec',[1 nr]),vstate');
set(hl(1),'color','r','linewidth',2);
set(hl(2:2:end),'color','g');
set(hl(3:2:end),'color','m');

hlsw=line(xvec,ones(1,vtp1)*NaN,'color','k','linewidth',2);

if ~autoscale
ylim([-myylim myylim]);
else
    ylim([-2 2]);
end;

xlim([0 max(xvec)]);

if starthard
    line([0 0],get(gca,'ylim'),'linewidth',3,'color','k','clipping','off');
end;
if endhard
    line([vtlen vtlen],get(gca,'ylim'),'linewidth',3,'color','k','clipping','off');
end;


hlcp=line(0,0,'marker','o','markersize',12);

set(gcf,'doublebuffer','on')
maxbuf=zeros(bign,1);

set(gca,'fontsize',12,'fontweight','bold');
ylabel('Sound pressure');
xlabel('Distance from glottis (cm)');
htit=title('0');
currow=1;
curcol=1;
curdir=1;

set(gcf,'userdata',0);
set(gcf,'windowbuttondownfcn',{@buttoncb});
disp('Click in figure to pause; then hit any key to resume');
for ipi=1:bign

buttonpressed=get(gcf,'userdata');
if buttonpressed
    drawnow;
    pause;
    set(gcf,'userdata',0)
end;

    newv=cos(2*pi*f*((ipi-1)/vtp));
if truncinput
    if ipi>(vtp*(1/f))
    newv=0;
end;
end;
    
    vstate(1,2:vtp1)=vstate(1,1:vtp);
    vstate(1,1)=newv;
    
    %even lines reflect back, odd lines reflect forwards
    for ir=2:nr
        if mod(ir,2)==0
            vstate(ir,1:vtp)=vstate(ir,2:vtp1);
            if endhard 
                vstate(ir,end)=vstate(ir-1,end);
            else
                vstate(ir,end)=-vstate(ir-1,end);
            end;
        else
            vstate(ir,2:vtp1)=vstate(ir,1:vtp);
            if starthard 
                vstate(ir,1)=vstate(ir-1,1);
            else
                vstate(ir,1)=-vstate(ir-1,1);
            end;
        end;
    end;
    sw=nansum(vstate.*mydecay);
    maxbuf(ipi)=max(abs(sw));
    allmax=max(maxbuf);
    pbuf=(vstate.*mydecay)';
    curamp=pbuf(curcol,currow);        %rows and columns swapped
    set(hlcp,'ydata',curamp,'xdata',xvec(curcol));
    
    set(hl,'visible','off');
    lv=find(~all(isnan(pbuf)));
    nl=length(lv);
    for il=1:nl
        myline=lv(il);
        set(hl(myline),'ydata',pbuf(:,myline),'visible','on');
    end;
    if nl>1
        set(hlsw,'ydata',sw);
    end;

    set(htit,'string',['t (ms) : ' num2str(round(timeincms*ipi*10)/10)]);
    
    doscale=0;
    curcol=curcol+curdir;
    if curdir==1
        if curcol>vtp1
            currow=currow+1;
            curdir=-1;
            curcol=vtp;
            doscale=1;
        end;
    else
        if curcol<1
            currow=currow+1;
            curdir=1;
            curcol=2;
            doscale=1;
        end;
    end;
    
    if autoscale
        if doscale
        mm=round(allmax)+1;
%        disp(mm);
    oldlim=get(gca,'ylim');
    if mm>oldlim(2)
        ylim([-mm mm]);
    end;
    
    end;
    
    end;
            
    drawnow;
    if isempty(mywait)
        pause
    else
        
    pause(mywait);
end;

end;
keyboard;
function buttoncb(obj,data);
set(obj,'userdata',1);
