function l1d04(tmax, iplot, Delay, iN, nL, iL, Z, boundary_type, ...
               pulse_type, rt)
%
  global d u;

% $$$   tmax  = 350;
% $$$   tmaxi = input(['tmax (', int2str(tmax), ') = ']);
% $$$   if ~isempty(tmaxi), tmax = tmaxi;  end;
% $$$ 
% $$$   iplot  = 1;
% $$$   iploti = input(['iplot (', int2str(iplot), ') = ']);
% $$$   if ~isempty(iploti), iplot = iploti;  end;
% $$$ 
% $$$   Delay = 0;
% $$$   Delayi = input(['Time delay (', int2str(Delay), ') = ']);
% $$$   if ~isempty(Delayi), Delay = Delayi;  end;
% $$$   %  printf("Time delay (msec) = ");
% $$$ 
% $$$   iN = 200;
% $$$   iNi = input(['iN (', int2str(iN), ') = ']);
% $$$   if ~isempty(iNi), iN = iNi;  end;
% $$$ 
% $$$   nL = 2;
% $$$   nLi = input(['number of layers (', int2str(nL), ') = ']);
% $$$   if ~isempty(nLi), nL = nLi;  end;
% $$$ 
% $$$   iL = 100;
% $$$   iLi = input(['thickness of layers(1:nL-1) (', int2str(iL), ') = ']);
% $$$   if ~isempty(iLi), iL = iLi;  end;
% $$$ 
% $$$   Z = [1 2];
% $$$   Zi = input('impedance of layers(1:nL) ([1 2]) = ');
% $$$   if ~isempty(iLi), Z = Zi;  end;
% $$$ 
% $$$   boundary_type = 1;
% $$$   bti = input(['boundary (refl.=1, abs.=2, zero stress=3, periodic=4) (', ...
% $$$                 int2str(boundary_type), ') = ']);
% $$$   if ~isempty(bti), boundary_type = bti;  end;
% $$$ 
% $$$   pulse_type = 2;
% $$$   pti = input(['input pulse (gauss=1, gradient=2, sinus=3) (', ...
% $$$                 int2str(pulse_type), ') = ']);
% $$$   if ~isempty(pti), pulse_type = pti;  end;
% $$$ 
% $$$   rt = 1;
% $$$   rti = input(['right (0,1) (', int2str(rt), ') = ']);
% $$$   if ~isempty(rti), rt = rti;  end;
  lt = 1 - rt;

  h  = zeros(1,iN);
  d  = h; u  = h;
  d1 = d; d2 = d;

  t1 = ones(1,iN);
  t2 = t1;
  t1(iN) = 0;
  t2(1)  = 0;

%  t1(1)  = 2;   % boundary - vacuuum
%  t2(iN) = 2;

  if (iL > 1)
    i0 = 1;
    for i = 1:nL-1
      i0 = i0 + iL(i);
      zz = Z(i) + Z(i+1);
      t1(i0) = 2 * Z(i+1) ./ zz;
      t2(i0) = 2 * Z(i) ./ zz;
    end;
  end;

  switch pulse_type
    case 1,  gauss;
    case 2,  grad ;
    case 3,  sinus;
  end;

  u = d;
  d1(1:end-1) = d(2:end);
  d2(2:end)   = d(1:end-1);
  h = (rt .* d1 + lt .* d2) ./ (rt + lt);

  plot_signal_init;

  t0 = clock;
  for it=0:tmax-1
    d1(1:end-1) = d(2:end);
    d1(end)     = d1(end-1);
    d2(2:end)   = d(1:end-1);
    d2(1)       = d2(2);
    u = t1 .* d1 + t2 .* d2 - h;

    switch boundary_type
      case 1,  % vacuum - reflection
	      u(1)  = t1(1) * d(2) - h(1);
	      u(iN) =  t2(iN) * d(iN - 1) - h(iN);

      case 2, % absorption
	      u(1)  = d(2);
	      u(iN) = d(iN - 1);

      case 3, % zero stress
	      % the same as boundary_type = 1 and t1(1) = 2; t2(iN) = 2;
	      u(1)  = u(2);
	      u(iN) = u(iN - 1);
      case 4, % periodic
	      u(1)  = u(iN-1);
	      u(iN) = u(2);
    end;  % switch

    if isequal(it, floor(it/iplot) * iplot)
      if isinf(Delay),  pause; else  pause(Delay); end
      plot_signal(it)
    end;

    h = d;
    d = u;

  end  % for it=0:tmax-1
  fprintf('elapsed time = %6.2f \n', etime(clock,t0))
  
%
% end l1d

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gauss()
%
  global d;

  ampl = 1;
  ai = input(['amplitude (', int2str(ampl), ') = ']);
  if ~isempty(ai), ampl = ai;  end;

  i0 = 25;
  i0i = input(['center (', int2str(i0), ') = ']);
  if ~isempty(i0i), i0 = i0i;  end;

  sigma = 7;
  sigmai = input(['sigma (', int2str(sigma), ') = ']);
  if ~isempty(sigmai), sigma = sigmai;  end;

%  d = ampl/ (2 .* sigma.^2) .* exp(-((1:length(d))-i0).^2);
  d = ampl .* exp(-0.5.*((1:length(d))-i0).^2 ./ sigma.^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function grad()
%
  global d;

  ampl = 1;
  ai = input(['amplitude (', int2str(ampl), ') = ']);
  if ~isempty(ai), ampl = ai;  end;

  i1 = 30;
  i1i = input(['imin (', int2str(i1), ') = ']);
  if ~isempty(i1i), i1 = i1i;  end;

  i2 = 50;
  i2i = input(['imax (', int2str(i2), ') = ']);
  if ~isempty(i2i), i2 = i2i;  end;

  d(i1:i2) = ampl.* ones(1,i2-i1+1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sinus()
%
  global d;

  ampl = 1;
  ai = input(['amplitude (', int2str(ampl), ') = ']);
  if ~isempty(ai), ampl = ai;  end;

  i1 = 20;
  i1i = input(['imin (', int2str(i1), ') = ']);
  if ~isempty(i1i), i1 = i1i;  end;

  i2 = 30;
  i2i = input(['imax (', int2str(i2), ') = ']);
  if ~isempty(i2i), i2 = i2i;  end;

  phase = 1;
  phasei = input(['phase (', int2str(phase), ') = ']);
  if ~isempty(phasei), phase = phasei;  end;

  d(i1:i2) = ampl .* sin(2 .* pi .*((i1:i2)-i1) ./(i2-i1) + phase);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_signal(time)
%
  global u h_line h_line_sep;

%  plot(u, 'LineWidth', 2)
  set(h_line, 'ydata', u)
  axs = axis;
  axs(3) = 1.1*min(u(:));
  axs(4) = 1.1*max(u(:));
  axis([axs(1) axs(2) axs(3) axs(4)]);
  set(h_line_sep, 'ydata', [axs(3) axs(4)])
%  text(length(u)/10, 0.9*(max(u) - min(u)), ['time = ', num2str(time)])
  title(['time = ', num2str(time)])
  drawnow;
%  pause;
%  pause(n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_signal_init()
%
  global u h_line h_line_sep;

  cla;
  h_line = line('xdata', 1:length(u), 'ydata', u, 'LineWidth',2, 'color','b');
  axs = axis;
  axs(3) = 1.1*min(u(:));
  axs(4) = 1.1*max(u(:));
  axis([axs(1) axs(2) axs(3) axs(4)]);
  h_line_sep=line('XData',[120 120], 'YData', [axs(3) axs(4)], 'Color', 'r');
  text(107, 0.05, 'Z=1')
  text(125, 0.05, 'Z=2')
  pause;
% ======== EOF ===========
