function [R_stencil,L_stencil,P_stencil]=lfa(pdeflag,smoothingflag,restrictionflag,prolongationflag,coarsegridflag,plottingflag);
% Local Fourier analysis for hyperbolic operator with continuous bilinear
% finite element discretization.
if(pdeflag==0)
  disp('***');
  disp('usage: lfa(pdeflag,smoothingflag,restrictionflag,prolongationflag,plottingflag)');
  disp(' ');
  disp('           pdeflag:          1 = Poisson');
  disp('                             2 = Hyperbolic');
  disp('           smoothingflag:    1 = GS-LEX');
  disp('                             2 = w-JAC');
  disp('                             2 = w-GS-LEX');
  disp('           restrictionflag:  1 = Full Weighting');
  disp('                             2 = Injection');
  disp('           prolongationflag: 1 = standard');
  disp('           coarsegridflag:   1 = Fine Grid');
  disp('                             2 = Galerkin Formulation RAP');
  disp('           plottingflag:     0 = off');
  disp('                             1 = on');
  disp('***');
  break;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=100;
M=51;
[i1,i2]=meshgrid(-1:1:1);
h_all=2*pi/N;
h_high=(pi/2)/(N/4);
h_low=pi/M;
[T1_all,T2_all]=meshgrid(-pi:h_all:pi);
[T1_high_w,T2_high_w]=meshgrid(-pi:h_high:-pi/2,-pi:h_all:pi);
[T1_high_e,T2_high_e]=meshgrid(pi/2:h_high:pi,-pi:h_all:pi);
[T1_high_s,T2_high_s]=meshgrid(-pi:h_all:pi,-pi:h_high:-pi/2);
[T1_high_n,T2_high_n]=meshgrid(-pi:h_all:pi,pi/2:h_high:pi);
[T1_low,T2_low]=meshgrid(-pi/2:h_low:pi/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L stencil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(pdeflag==1)
  pdestring='Poisson';
  L_stencil=[0 -1 0;-1 4 -1;0 -1 0];
elseif(pdeflag==2)
  pdestring='Hyperbolic';
  L_stencil=[-5/12 -1/6 1/12;-1/6 4/3 -1/6;1/12 -1/6 -5/12];
  %L_stencil=[1/12 -1/6 -5/12;-1/6 4/3 -1/6;-5/12 -1/6 1/12];
end
syms T1 T2;           % Theta 1 and Theta 2
L_symbol=simplify(sum(sum(L_stencil.*(exp(i*i1*T1+i*i2*T2)))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smoothing stencil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L_plus=L_stencil;
L_minus=L_stencil;
if(smoothingflag==1)
  % GS-LEX
  smootherstring='GS-LEX';
  L_plus(1,:)=0;
  L_plus(2,3)=0;
  L_minus(3,:)=0;
  L_minus(2,1:2)=0;
elseif(smoothingflag==2)
  % w-JAC
  omega=2/3;
  smootherstring=['w-JAC with w=' num2str(omega)];
  L_plus=0*L_plus;
  L_plus(2,2)=L_minus(2,2);
  L_plus=(1/omega)*L_plus;
  L_minus=L_minus-L_plus;
elseif(smoothingflag==3)
  % w-GS-LEX
  omega=.95;
  smootherstring=['w-GS-LEX with w=' num2str(omega)];
  L_plus(1,:)=0;
  L_plus(2,3)=0;
  L_plus(2,2)=(1/omega)*L_plus(2,2);
  L_minus(3,:)=0;
  L_minus(2,1)=0;
  L_minus(2,2)=L_stencil(2,2)-L_plus(2,2);
end
L_plus_symbol=simplify(sum(sum(L_plus.*(exp(i*i1*T1+i*i2*T2)))));
L_minus_symbol=simplify(sum(sum(L_minus.*(exp(i*i1*T1+i*i2*T2)))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smoothing analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Amp_Factor=-L_minus_symbol/L_plus_symbol;
A_high_w=subs(Amp_Factor,{T1,T2},{T1_high_w,T2_high_w});
A_high_e=subs(Amp_Factor,{T1,T2},{T1_high_e,T2_high_e});
A_high_s=subs(Amp_Factor,{T1,T2},{T1_high_s,T2_high_s});
A_high_n=subs(Amp_Factor,{T1,T2},{T1_high_n,T2_high_n});
A_low=subs(Amp_Factor,{T1,T2},{T1_low,T2_low});
Smoothing_Factor=max(max([max(max(abs(A_high_w))),max(max(abs(A_high_e))),max(max(abs(A_high_s))),max(max(abs(A_high_n)))]));
F_high_w=subs(L_symbol,{T1,T2},{T1_high_w,T2_high_w});
F_high_e=subs(L_symbol,{T1,T2},{T1_high_e,T2_high_e});
F_high_s=subs(L_symbol,{T1,T2},{T1_high_s,T2_high_s});
F_high_n=subs(L_symbol,{T1,T2},{T1_high_n,T2_high_n});
F_low=subs(L_symbol,{T1,T2},{T1_low,T2_low});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Restriction and Prolongation stencils
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(restrictionflag==1)
  restrictionstring='Full Weighting';
  R_stencil=(1/4)*[1 2 1;2 4 2;1 2 1];
elseif(restrictionflag==2)
  restrictionstring='Full Weighting';
  R_stencil=zeros(3);
  R_stencil(2,2)=4;
end
if(prolongationflag==1)
  prolongationstring='Bilinear Interpolation';
  P_stencil=(1/16)*[1 2 1;2 4 2;1 2 1];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Coarse grid stencil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(coarsegridflag==1)
  coarsegridstring='Fine grid';
  LH_stencil=L_stencil;
elseif(coarsegridflag==1)
  coarsegridstring='Galerkin Formulation RAP';
  LH_stencil=R_stencil*L_stencil*P_stencil;
end
R_symbol=simplify(sum(sum(R_stencil.*(exp(i*i1*T1+i*i2*T2)))));
P_symbol=simplify(sum(sum(P_stencil.*(exp(i*i1*T1+i*i2*T2)))));
LH_symbol=simplify(sum(sum(LH_stencil.*(exp(i*i1*T1+i*i2*T2)))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% low freqs and their harmonics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T1_00=T1_low;
T2_00=T2_low;
T2_10=T2_low;
T1_01=T1_low;
for j1=1:size(T1_00,1)
  for j2=1:size(T1_00,2)
    if(T1_00(j1,j2) < 0)
      T1_10(j1,j2)=T1_00(j1,j2)+pi;
      T1_11(j1,j2)=T1_00(j1,j2)+pi;
    elseif(T1_00(j1,j2) >= 0)
      T1_10(j1,j2)=T1_00(j1,j2)-pi;
      T1_11(j1,j2)=T1_00(j1,j2)-pi;
    end
  end
end
for j1=1:size(T2_00,1)
  for j2=1:size(T2_00,2)
    if(T2_00(j1,j2) < 0)
      T2_01(j1,j2)=T2_00(j1,j2)+pi;
      T2_11(j1,j2)=T2_00(j1,j2)+pi;
    elseif(T2_00(j1,j2) >= 0)
      T2_01(j1,j2)=T2_00(j1,j2)-pi;
      T2_11(j1,j2)=T2_00(j1,j2)-pi;
    end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% building K
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lh_00=subs(L_symbol,{T1,T2},{T1_00,T2_00});
Lh_01=subs(L_symbol,{T1,T2},{T1_01,T2_10});
Lh_10=subs(L_symbol,{T1,T2},{T1_10,T2_01});
Lh_11=subs(L_symbol,{T1,T2},{T1_11,T2_11});
LH_00=subs(LH_symbol,{T1,T2},{T1_00,T2_00});
LH_01=subs(LH_symbol,{T1,T2},{T1_01,T2_10});
LH_10=subs(LH_symbol,{T1,T2},{T1_10,T2_01});
LH_11=subs(LH_symbol,{T1,T2},{T1_11,T2_11});
LH_00_inv=1./LH_00;
LH_01_inv=1./LH_01;
LH_10_inv=1./LH_10;
LH_11_inv=1./LH_11;
R_00=subs(R_symbol,{T1,T2},{T1_00,T2_00});
R_01=subs(R_symbol,{T1,T2},{T1_01,T2_10});
R_10=subs(R_symbol,{T1,T2},{T1_10,T2_01});
R_11=subs(R_symbol,{T1,T2},{T1_11,T2_11});
P_00=subs(P_symbol,{T1,T2},{T1_00,T2_00});
P_01=subs(P_symbol,{T1,T2},{T1_01,T2_10});
P_10=subs(P_symbol,{T1,T2},{T1_10,T2_01});
P_11=subs(P_symbol,{T1,T2},{T1_11,T2_11});
S_symbol=Amp_Factor;
S_00=subs(S_symbol,{T1,T2},{T1_00,T2_00});
S_01=subs(S_symbol,{T1,T2},{T1_01,T2_10});
S_10=subs(S_symbol,{T1,T2},{T1_10,T2_01});
S_11=subs(S_symbol,{T1,T2},{T1_11,T2_11});
I_00=ones(size(S_00));
I_01=ones(size(S_00));
I_10=ones(size(S_00));
I_11=ones(size(S_00));
COEF_00=P_00.*LH_00_inv.*R_00;
COEF_01=P_01.*LH_01_inv.*R_01;
COEF_10=P_10.*LH_10_inv.*R_10;
COEF_11=P_11.*LH_11_inv.*R_11;
COEF=COEF_00+COEF_01+COEF_10+COEF_11;
K_00=I_00-COEF.*Lh_00;
K_01=I_01-COEF.*Lh_01;
K_10=I_10-COEF.*Lh_10;
K_11=I_11-COEF.*Lh_11;
M_00=S_00.*K_00.*S_00;
M_01=S_01.*K_01.*S_01;
M_10=S_10.*K_10.*S_10;
M_11=S_11.*K_11.*S_11;

figure(2);clf;
plot([0 0],[-pi/2 pi/2],'k--',[-pi/2 pi/2],[0 0],'k--');
hold on;
for j1=1:size(M_00)
  for j2=1:size(M_00)
    M_temp=zeros(4,4);
    M_temp(1,1)=M_00(j1,j2);
    M_temp(2,2)=M_11(j1,j2);
    M_temp(3,3)=M_10(j1,j2);
    M_temp(4,4)=M_01(j1,j2);
    rho_all(j1,j2)=max(abs(eig(M_temp)));
    if(rho_all(j1,j2) < 4)
      %disp(['j1=' num2str(j1) ', j2=' num2str(j2) ', rho=' num2str(rho_all(j1,j2))]);
    else
      plot(T1_00(j1,j2),T2_00(j1,j2),'x');
      axis([-pi/2 pi/2 -pi/2 pi/2]);
      %disp(['j1=' num2str(j1) ', j2=' num2str(j2) ', rho=' num2str(rho_all(j1,j2)) '  XXXXXXXXXXXXX']);
    end
  end
end
rho=max(max(rho_all));



%
% Plotting
%
A_max=max([max(max(A_high_w)),max(max(A_high_e)),max(max(A_high_s)),max(max(A_high_n)),max(max(A_low))]);
A_min=min([min(min(A_high_w)),min(min(A_high_e)),min(min(A_high_s)),min(min(A_high_n)),min(min(A_low))]);
if(plottingflag==1)
  figure
  %subplot(2,1,1),
  hold on;axis([-pi pi -pi pi A_min A_max]);
  surf(T1_low,T2_low,abs(A_low))
  caxis([A_min,A_max]);
  shading interp;
  %subplot(2,1,2),
  hold on;axis([-pi pi -pi pi A_min A_max]);
  surf(T1_high_w,T2_high_w,abs(A_high_w))
  surf(T1_high_e,T2_high_e,abs(A_high_e))
  surf(T1_high_s,T2_high_s,abs(A_high_s))
  surf(T1_high_n,T2_high_n,abs(A_high_n))
  caxis([A_min,A_max]);
  shading interp;
end
disp('--------------------------------------------');
disp(['PDE:                  ' pdestring]);
disp(['Smoother:             ' smootherstring]);
disp(['Restriction:          ' restrictionstring]);
disp(['Prolongation:         ' prolongationstring]);
disp(['Coarse Grid Operator: ' coarsegridstring]);
disp('............................................');
disp(['Smoothing Factor = ' num2str(Smoothing_Factor)]);
disp('--------------------------------------------');