% 1DFDTD     
% Demonstrates 1-D FDTD initial state formation. 
% Please edit the FLAGS for demonstration of different cases.
% BASED ON "1-D Digital Waveguide Modeling for Improved Sound Synthesis"
% Editer: Chun-Feng Lai

% FLAG DESCRIPTION
  % FORCEFLAG: Choice between initial state formation (0) and external force (1).
  % EXC_TYPE: 
  %  1: Traveling impulse given by Eqs(5-8)
  %  2: Boxcar given by Eqs (9-10)
  % SIMULATION_MODE
  %  0: Static (such as Fig.6)
  %  1: Dynamic (animation). 
  
  % FLAGS
 FORCEFLAG = 0;       % <0,1>
 EXC_TYPE = 1;        % <1,2,3>
 SIMULATION_MODE = 1; % <0,1>
 
 % Screen update in simulation mode.
 update = 10;   % 10 for Static, 1 for Dynamic recommended
 MAXSTEP = 1000;% 100 for Static, 1000 for Dynamic recommended
 
% General simulation parameters
 fs = 44100;
 ff = 441;
 n = fix(fs/ff);
 c = 340;
 
% LOSSES
 g = 0.9998;a = -g^2;
  
% Force excitation parameters
 PLUCKPOINT = 10/50; % Relative ration between <0-1> 
 FW = 5; % Force width [odd samples]
 FW2 = fix(FW/2);
 AMP = .025/FW;
 
 Plp = fix(n*PLUCKPOINT);
 
 % Prepare the force density matrix
 F = zeros(MAXSTEP,n);
  
% INITIALLY RELAXED STATES
 y_init2 = zeros(1,n);
 y_init1 = zeros(1,n);

 y = zeros(1,n);
 
if FORCEFLAG
  F(:,Plp-FW2:Plp+FW2) = AMP*(hanning(FW)*ones(1,MAXSTEP))'; 
else
 switch EXC_TYPE
 case(1);
  % EXCITATION 1
  % ORIGINAL FORWARD TRAVELING PLUCK Eq.(5,6)
  y_init1(n/2) = y_init1(n/2)+0.5;
  y_init2(n/2-1) = y_init2(n/2-1) + 0.5;
  
  % BACKWARD TRAVELING PLUCK Eq.(7,8)
  y_init1(n/2) = y_init1(n/2)+0.5;
  y_init2(n/2+1) = y_init2(n/2+1)+0.5;
  
 case(2)
  % ORIGINAL EXCITATION 2 Eq (9,10)
  y_init1(n/2) = y_init1(n/2)+0.5;
  y_init1(n/2+1) = y_init1(n/2+1)+0.5;
 end
end

figure(1);clf; 

 % Comment the following line out if your system does not support
 % OpenGL rendering for flickerless simulation
 set(gcf,'Renderer','OpenGL')

% SIMULATION
for i=0:MAXSTEP-1;
  i
  y(2:n-1) = g*(y_init1(1:n-2) + y_init1(3:n)) + a*y_init2(2:n-1)+...
     FORCEFLAG*F(i+1,2:n-1);
  
  y_init2(2:n-1) = y_init1(2:n-1);
  y_init1(2:n-1) = y(2:n-1);
  if (rem(i,update)==0)
    switch SIMULATION_MODE
     case 0
      plot(y+1*i/update); hold on;
      h=text(62+3*i/update,max(y+1*i/update+.125),['n = ' num2str(i)]);
      set(h,'FontName','Times')
     case 1
      plot(y(2:n-1));axis([2 n-1 -1 1]);drawnow;pause(.5);
    end
  end
end

set(gca,'FontName','Times');
set(gca,'FontSize',16);
xlabel('k');
ylabel('y');