% Determine stability limits for 3d wave-propagation algorithms, as
% described in the paper
%  A wave propagation method for three-dimensional hyperbolic conservation laws
%  by Jan Olav Langseth and R. J. LeVeque,  submitted to J. Comput. Phys., 1999

% Results are output to the file stabXXX.dat where the three X's
% are replaced by the values of meth1, meth2, meth3 defining the method.
%
% Copyright Jan Olav Langseth and Randall J. LeVeque, 1999
%

meth1 = 1;
meth2 = 0;
meth3 = 0;

nckh = 8; % number of values of c*k/h
npts = 10;  % number of xi,eta,omega points to check for each ckh
	   % (checks on an  npts x npts x npts  array)

amplfac = nan*ones(nckh,5);
ckhv    = nan*ones(nckh);

% Array containing the Courant numbers
ckhv = [ 0.3 0.33 0.36 0.39 0.42 0.45 0.48 0.51];

fname = ['stab',num2str(meth1),num2str(meth2),num2str(meth3),'.dat'];
fid = fopen(fname,'w');

% Euler (primitive variables)
 neqn = 5;
 u = 0;
 v = 0;
 w = 0;
 r = 1;
 p = 1;
 g = 1;

 a = [u r 0 0 0; 0 u 0 0 1/r; 0 0 u 0 0; 0 0 0 u 0; 0 g*p 0 0 u];
 b = [v 0 r 0 0; 0 v 0 0 0; 0 0 v 0 1/r; 0 0 0 v 0; 0 0 g*p 0 v];
 c = [w 0 0 r 0; 0 w 0 0 0; 0 0 w 0 0; 0 0 0 w 1/r; 0 0 0 g*p w];

% acoustics:
% neqn = 4;
% a = [0 1 0 0; 1 0 0 0; 0 0 0 0; 0 0 0 0];
% b = [0 0 1 0; 0 0 0 0; 1 0 0 0; 0 0 0 0];
% c = [0 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 0];

% advection:
% neqn = 1;
% a = 1;
% b = 1;
% c = 1;
[r,lam] = eig(a);
lamp = max(lam,0);
lamm = min(lam,0);
ap = r*lamp*inv(r);
am = r*lamm*inv(r);
amax = max(max(abs(lam)));

[r,lam] = eig(b);
lamp = max(lam,0);
lamm = min(lam,0);
bp = r*lamp*inv(r);
bm = r*lamm*inv(r);
bmax = max(max(abs(lam)));

[r,lam] = eig(c);
lamp = max(lam,0);
lamm = min(lam,0);
cp = r*lamp*inv(r);
cm = r*lamm*inv(r);
cmax = max(max(abs(lam)));


sixth = 1/6;

% maximum of 1D Courant numbers
cc = max([amax bmax cmax]);

for jckh = 1:nckh
  ckh = ckhv(jckh)/cc;

  ckh*cc            % Just for output

  maxampl = 0;
  for jxi = 1:npts
    xi = 2*pi*jxi / npts;

    for jeta = 1:npts
      eta = 2*pi*jeta / npts;

      for jomega = 1:npts
      omega = 2*pi*jomega / npts;

      emxi = exp(-i*xi);            
      epxi = exp(i*xi);

      emeta = exp(-i*eta);            
      epeta = exp(i*eta);

      emomega = exp(-i*omega);            
      epomega = exp(i*omega);

      % donor cell:
      f = ap*emxi + am;
      % transverse corrections:
      if (meth2 > 0)
        % AB 
        f = f - 0.5*ckh * (ap*bp*(1-emeta)*emxi + ap*bm*(epeta-1)*emxi);
        f = f - 0.5*ckh * (am*bp*(1-emeta) + am*bm*(epeta-1));
        % AC
        f = f - 0.5*ckh * (ap*cp*(1-emomega)*emxi + ap*cm*(epomega-1)*emxi);
        f = f - 0.5*ckh * (am*cp*(1-emomega) + am*cm*(epomega-1));
      end

      if (meth3 > 0)
	% ACB + ABC
        f = f + sixth*ckh*ckh*(1-emomega)*(1-emeta)*emxi * ap*(cp*bp+bp*cp) ; 
        f = f + sixth*ckh*ckh*(epomega-1)*(1-emeta)*emxi * ap*(cm*bp+bp*cm) ;
        f = f + sixth*ckh*ckh*(1-emomega)*(epeta-1)*emxi * ap*(cp*bm+bm*cp) ; 
        f = f + sixth*ckh*ckh*(epomega-1)*(epeta-1)*emxi * ap*(cm*bm+bm*cm) ;
        f = f + sixth*ckh*ckh*(1-emomega)*(1-emeta) * am*(cp*bp+bp*cp) ; 
        f = f + sixth*ckh*ckh*(epomega-1)*(1-emeta) * am*(cm*bp+bp*cm) ; 
        f = f + sixth*ckh*ckh*(1-emomega)*(epeta-1) * am*(cp*bm+bm*cp) ; 
        f = f + sixth*ckh*ckh*(epomega-1)*(epeta-1) * am*(cm*bm+bm*cm) ; 
        end
      % second order terms:
      if (meth1 == 2) 
        cqxx = (ap-am)*(eye(neqn)-ckh*(ap-am))*(1-emxi);
        f = f + 0.5*cqxx;
        end
      
      % donor cell:
      g = bp*emeta + bm;
      % transverse corrections:
      if (meth2 > 0) 
        % BA
        g = g - 0.5*ckh * (bp*ap*(1-emxi)*emeta + bp*am*(epxi-1)*emeta);
        g = g - 0.5*ckh * (bm*ap*(1-emxi) + bm*am*(epxi-1));
        % BC
        g = g - 0.5*ckh * (bp*cp*(1-emomega)*emeta + bp*cm*(epomega-1)*emeta);
        g = g - 0.5*ckh * (bm*cp*(1-emomega) + bm*cm*(epomega-1));
      end

      if( meth3 > 0 )
	% BAC + BCA
        g = g + sixth*ckh*ckh*(1-emxi)*(1-emomega)*emeta * bp*(ap*cp+cp*ap) ; 
        g = g + sixth*ckh*ckh*(epxi-1)*(1-emomega)*emeta * bp*(am*cp+cp*am) ; 
        g = g + sixth*ckh*ckh*(1-emxi)*(epomega-1)*emeta * bp*(ap*cm+cm*ap) ; 
        g = g + sixth*ckh*ckh*(epxi-1)*(epomega-1)*emeta * bp*(am*cm+cm*am) ; 
        g = g + sixth*ckh*ckh*(1-emxi)*(1-emomega) * bm*(ap*cp+cp*ap) ; 
        g = g + sixth*ckh*ckh*(epxi-1)*(1-emomega) * bm*(am*cp+cp*am) ; 
        g = g + sixth*ckh*ckh*(1-emxi)*(epomega-1) * bm*(ap*cm+cm*ap) ; 
        g = g + sixth*ckh*ckh*(epxi-1)*(epomega-1) * bm*(am*cm+cm*am) ; 
      end

      % second order terms:
      if (meth1 == 2) 
        cqyy = (bp-bm)*(eye(neqn)-ckh*(bp-bm))*(1-emeta);
        g = g + 0.5*cqyy;
        end
      
      % donor cell:
      h = cp*emomega + cm;
      % transverse corrections:
      if (meth2 > 0)
        % CA 
        h = h - 0.5*ckh * (cp*ap*(1-emxi)*emomega + cp*am*(epxi-1)*emomega);
        h = h - 0.5*ckh * (cm*ap*(1-emxi) + cm*am*(epxi-1));
        % CB
        h = h - 0.5*ckh * (cp*bp*(1-emeta)*emomega + cp*bm*(epeta-1)*emomega);
        h = h - 0.5*ckh * (cm*bp*(1-emeta) + cm*bm*(epeta-1));
      end

      if (meth3 > 0)
	% CAB + CBA
        h = h + sixth*ckh*ckh*(1-emxi)*(1-emeta)*emomega * cp*(ap*bp+bp*ap) ;
        h = h + sixth*ckh*ckh*(epxi-1)*(1-emeta)*emomega * cp*(am*bp+bp*am) ;
        h = h + sixth*ckh*ckh*(1-emxi)*(epeta-1)*emomega * cp*(ap*bm+bm*ap) ; 
        h = h + sixth*ckh*ckh*(epxi-1)*(epeta-1)*emomega * cp*(am*bm+bm*am) ;
        h = h + sixth*ckh*ckh*(1-emxi)*(1-emeta) * cm*(ap*bp+bp*ap) ; 
        h = h + sixth*ckh*ckh*(epxi-1)*(1-emeta) * cm*(am*bp+bp*am) ;
        h = h + sixth*ckh*ckh*(1-emxi)*(epeta-1) * cm*(ap*bm+bm*ap) ; 
        h = h + sixth*ckh*ckh*(epxi-1)*(epeta-1) * cm*(am*bm+bm*am) ;
      end

      % second order terms:
      if (meth1 ==  2) 
        cqzz = (cp-cm)*(eye(neqn)-ckh*(cp-cm))*(1-emomega);
        h = h + 0.5*cqzz;
      end
      
      % transverse propagation of second order terms:
      if (meth2 == 2)
        % BAA
        g = g - (bm + bp*emeta)*0.5*ckh*cqxx * (epxi-1);
        % BCC
        g = g - (bm + bp*emeta)*0.5*ckh*cqzz * (epomega-1);
        % ABB
        f = f - (am + ap*emxi)*0.5*ckh*cqyy * (epeta-1);
        % ACC
        f = f - (am + ap*emxi)*0.5*ckh*cqzz * (epomega-1);
        % CBB
        h = h - (cm + cp*emomega)*0.5*ckh*cqyy * (epeta-1);
        % CAA
        h = h - (cm + cp*emomega)*0.5*ckh*cqxx * (epxi-1);
      end

      if (meth3 == 2) 
        % ACBB + ABCC
        f = f - 0.25*ckh*ckh*am*cp*cqyy*(1-epeta)*(1-emomega);
        f = f - 0.25*ckh*ckh*am*cm*cqyy*(1-epeta)*(epomega-1);
        f = f - 0.25*ckh*ckh*ap*cp*cqyy*(1-epeta)*(1-emomega)*emxi;
        f = f - 0.25*ckh*ckh*ap*cm*cqyy*(1-epeta)*(epomega-1)*emxi;
        f = f - 0.25*ckh*ckh*am*bp*cqzz*(1-epomega)*(1-emeta);
        f = f - 0.25*ckh*ckh*am*bm*cqzz*(1-epomega)*(epeta-1);
        f = f - 0.25*ckh*ckh*ap*bp*cqzz*(1-epomega)*(1-emeta)*emxi;
        f = f - 0.25*ckh*ckh*ap*bm*cqzz*(1-epomega)*(epeta-1)*emxi;

        % BACC + BCAA
        g = g - 0.25*ckh*ckh*bm*ap*cqzz*(1-epomega)*(1-emxi);
        g = g - 0.25*ckh*ckh*bm*am*cqzz*(1-epomega)*(epxi-1);
        g = g - 0.25*ckh*ckh*bp*ap*cqzz*(1-epomega)*(1-emxi)*emeta;
        g = g - 0.25*ckh*ckh*bp*am*cqzz*(1-epomega)*(epxi-1)*emeta;
        g = g - 0.25*ckh*ckh*bm*cp*cqxx*(1-epxi)*(1-emomega);
        g = g - 0.25*ckh*ckh*bm*cm*cqxx*(1-epxi)*(epomega-1);
        g = g - 0.25*ckh*ckh*bp*cp*cqxx*(1-epxi)*(1-emomega)*emeta;
        g = g - 0.25*ckh*ckh*bp*cm*cqxx*(1-epxi)*(epomega-1)*emeta;

        % CBAA + CABB
        h = h - 0.25*ckh*ckh*cm*bp*cqxx*(1-epxi)*(1-emeta);
        h = h - 0.25*ckh*ckh*cm*bm*cqxx*(1-epxi)*(epeta-1);
        h = h - 0.25*ckh*ckh*cp*bp*cqxx*(1-epxi)*(1-emeta)*emomega;
        h = h - 0.25*ckh*ckh*cp*bm*cqxx*(1-epxi)*(epeta-1)*emomega;
        h = h - 0.25*ckh*ckh*cm*ap*cqyy*(1-epeta)*(1-emxi);
        h = h - 0.25*ckh*ckh*cm*am*cqyy*(1-epeta)*(epxi-1);
        h = h - 0.25*ckh*ckh*cp*ap*cqyy*(1-epeta)*(1-emxi)*emomega;
        h = h - 0.25*ckh*ckh*cp*am*cqyy*(1-epeta)*(epxi-1)*emomega;

        end

      ampl = eye(neqn) - ckh * ((epxi-1)*f + (epeta-1)*g + (epomega-1)*h);

      amplfactor = max(abs(eig(ampl)));
 
      [R,D] = eig(ampl);

      if (rank(R) ~= neqn & norm(ampl^100) > 10 & amplfactor < 1+1.0e-8 )
           'stoped'
           pause      
      end

      if ( amplfactor > maxampl )
           maxampl = amplfactor;
           maxxi    = xi;
           maxeta   = eta;
           maxomega = omega;
      end

      end % loop on omega
      end % loop on eta
    end   % loop on xi


    amplfac(jckh,1) = cc * ckh;
    amplfac(jckh,2) = maxampl;
    amplfac(jckh,3) = maxxi;
    amplfac(jckh,4) = maxeta;
    amplfac(jckh,5) = maxomega;

  end % loop on ckh

fprintf(fid,'%6.2f %12.8f %12.8f %12.8f %12.8f\n',amplfac');