%  Modified Wave Number Examples  THP 10/17/2001
%  Plots Modified Wave number and an example of a 
%  'k' wave number function being propagated at speed
%  'a' to demonstate the effect of phase (convective)
%  and amplitude (diffusive) errors for various
%  finite difference schemes.

clear;
jmax = 50;
dx = 2*pi/jmax;
x = [0:jmax-1]*dx;
k_dx = [0:jmax/2]*dx;
On_e = ones(size(k_dx,1),1);
Zer_o = zeros(size(k_dx,1),1);
k = 1;
a = 1.0;
alpha = 3.0;
nmax = 25;
menu2 = 0;

while menu2 >= 0

  menu2 = menu('',...
      ' RUN ',...
      ' new wave number k',...
      ' longer time ',...
      ' Change # of points',...
      ' Plot all mod wave equations',...
      ' Exit');

  if menu2 == 2
    k = Enter('Wave # ', k);
  elseif menu2 == 3
    a = Enter('Wave speed', a);
    nmax = Enter('Number of its', nmax);
  elseif menu2 == 4
    jmax = Enter(' jmax ',jmax);
    dx = 2*pi/jmax;
    x = [0:jmax-1]*dx;
    k_dx = [0:jmax/2]*dx;
    On_e = ones(size(k_dx,1),1);
    Zer_o = zeros(size(k_dx,1),1);
    clear ks_c_p ps_d_p;
  elseif menu2 == 6
   break;
  end

  if k > 21 break; end
  

  if menu2 == 1
    menu1 = menu('Wave Eq. Example',...
	'2nd O Central',...
	'1st O backward',...
	'2nd O backward',...
	'4th O Central Pade',...
	'4th O Central',...
	'3rd O Opt Pade', ...
    ' 4th-6th Gen Pade');
  
    kdx = k*dx;
    if menu1 == 1
      ks_c =sin(kdx)/dx;
      ks_c_p =sin(k_dx);
      ks_d = 0.0;
      ks_d_p = Zer_o;
      lab1 = ' O(\Delta x^2) Central';
      n_2 = 1;
    elseif menu1 == 2
      ks_c =sin(kdx)/dx;
      ks_c_p =sin(k_dx);
      ks_d = (1-cos(kdx))/dx;
      ks_d_p = exp(-(1-cos(k_dx)));
      lab1 = ' O(\Delta x^2) Backward';
      n_2 = 2;
    elseif menu1 == 3
      ks_c =(4*sin(kdx)-sin(2*kdx))/dx*0.5;
      ks_c_p =(4*sin(k_dx)-sin(2*k_dx))*0.5;
      ks_d = (3-4*cos(kdx)+cos(2*kdx))/dx*0.5;
      ks_d_p = exp(-(3-4*cos(k_dx)+cos(2*k_dx))*0.5);
      lab1 = ' O(\Delta x^2) Backward';
      n_2 = 2;
    elseif menu1 == 4
      ks_c =3*sin(kdx)/(2 + cos(kdx))/dx;
      ks_c_p =3*sin(k_dx)./(2 + cos(k_dx));
      ks_d = 0.0;
      ks_d_p = Zer_o;
      lab1 = ' O(\Delta x^4) Pade';
      n_2 = 1;
    elseif menu1 == 5
      ks_c =(8*sin(kdx)-sin(2*kdx))/6/dx;
      ks_c_p =(8*sin(k_dx)-sin(2*k_dx))/6;
      ks_d = 0.0;
      ks_d_p = Zer_o;
      lab1 = ' O(\Delta x^4) central';
      n_2 = 1;
    elseif menu1 == 6
      kss = (2-2*cos(kdx) + 3*i*sin(kdx))/(2 + cos(kdx) - i*sin(kdx))/dx;
      ks_c = imag(kss);
      ks_d = real(kss);
      kss_p = (2-2*cos(k_dx) + 3*i*sin(k_dx))./(2 + cos(k_dx) - i*sin(k_dx));
      ks_c_p = imag(kss_p);
      ks_d_p = exp(-real(kss_p));
      lab1 = ' O(\Delta x^3) Optimized Pade';
      n_2 = 2;
    elseif menu1 == 7
      alpha = Enter('alpha for Pade Scheme',alpha);
      C = (4.0-alpha)/3.0;  B = (4.0*alpha+2.0)/3.0;
      ks_c =(B*sin(kdx) + C*0.5*sin(2*kdx))/(alpha+2*cos(kdx))/dx;
      ks_c_p =(B*sin(k_dx) + C*0.5*sin(2*k_dx))./(alpha+2*cos(k_dx));
      ks_d = 0.0;
      ks_d_p = Zer_o;
      lab1 = [ 'Pade alpha= ' num2str(alpha)];
      n_2 = 1;
    end 
    ks = ks_d + i*ks_c;
    
    figure(1);
    subplot(2,n_2,1);
    plot(k_dx,k_dx,'k',k_dx,ks_c_p,'r',kdx,kdx,'bo','Markersize',10);
    title('Modified Wave Number plot');
    xlabel('k dx');
    ylabel('k^* dx');
    axis square;
    if n_2 == 2
      subplot(2,n_2,2);
      plot(k_dx,On_e,'k',k_dx,ks_d_p,'r',kdx,1.0,'bo','Markersize',10);
      title('Modified Wave Number plot');
      xlabel('k dx');
      ylabel('Amplitude Error');
      axis square;
    end
    subplot(2,1,2);
    t = a/nmax/dx;
      % $$$     u =sin(k*(x-t));
      % $$$     un = sin(k*(x-ks_c*t));
    u = (exp(i*k*(x-t))-exp(-i*k*(x-t)))*0.5/i;
    un = exp(-ks_d*t)*(exp(i*ks_c*(x-t))-exp(-i*ks_c*(x-t)))*0.5/i;
    h=plot(x,u,'r-o',x,un,'k+--');
    set(gcf,'doublebuffer','on');
    for n = 1:nmax
      t = a*n/nmax/dx;
      % $$$     u =sin(k*(x-t));
      % $$$     un = sin(k*(x-ks_c*t));
      u = (exp(i*k*(x-t))-exp(-i*k*(x-t)))*0.5/i;
      un = exp(-ks_d*t)*(exp(i*ks_c*(x-t))-exp(-i*ks_c*(x-t)))*0.5/i;

      set(h,'xdata',x,x,'ydata',u,un);
      legend('U exact',lab1);
      title([ 'Wave Solution For K = ' num2str(k)]);
      xlabel('x');
      ylabel('sin(k x)');
      drawnow;
      pause(0.1);
    end
   
  end
  
  if menu2 == 5
      figure(2);
      ks_c_p =sin(k_dx);
      plot(k_dx,k_dx,'k',k_dx,ks_c_p,'r');
      hold on;
      ks_c_p =(4*sin(k_dx)-sin(2*k_dx))*0.5;
      plot(k_dx,ks_c_p,'g');
      ks_c_p =3*sin(k_dx)./(2 + cos(k_dx));
      plot(k_dx,ks_c_p,'b');
      ks_c_p =(8*sin(k_dx)-sin(2*k_dx))/6;
      plot(k_dx,ks_c_p,'c');
      
      C = (4.0-alpha)/3.0;  B = (4.0*alpha+2.0)/3.0;
      ks_c_p =(B*sin(k_dx) + C*0.5*sin(2*k_dx))./(alpha+2*cos(k_dx));
      plot(k_dx,ks_c_p,'m');
      legend('2^{nd} O Central','2^{nd} O Backward',...
	  '4^{th} O Pade','4^{th} O Central',...
	  [ 'Pade: \alpha = ' num2str(alpha)]);
      title('Modified Wave Numbers');
      xlabel('k dx');
      ylabel('k^* dx');
  end
end