%
% @file fourier_smooth.m
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
%   Brute force computation of \mu = MAX |G| for
%   damped Jacobi and Gauss-Seidel relaxation
%
%   Find optimal damping factor \omega for Jacobi relaxations
%
if ~exist('tau','var'), tau=0; end
if ~exist('alpha','var'), alpha=1.0; end
dth=0.01;

theta1=-pi+dth:dth:pi;
theta2=-pi+dth:dth:pi;
nth1=length(theta1);
nth2=length(theta2);

sint12=sin(theta1./2).^2;
sint22=sin(theta2./2).^2;

sint1=sin(theta1);
sint2=sin(theta2);
[S1,S2]=meshgrid(sint1,sint2);
ctau=0.25*alpha*tau;

omega=0.5:0.002:1;
n=length(omega);
G = zeros(nth1,nth2);
for i=1:n
  c=2*omega(i)/(1+alpha^2);
  for i1=1:nth1
      for i2=1:nth2
          if or(abs(theta1(i1))>= pi/2, abs(theta2(i2)) >= pi/2);
              G(i1,i2) = abs(1-c*(sint12(i1)+alpha^2*sint22(i2) + ...
                                 ctau*sint1(i1)*sint2(i2)));
          end
      end
  end
  [gmax,imax]=max(G);
  [mu(i),jmax]=max(gmax);
  theta1_opt(i) = theta1(imax(jmax));
  theta2_opt(i) = theta2(jmax);
end
[mu_min,i_min]=min(mu);
omega_opt=omega(i_min);
str_title=sprintf(['omega = %.3f, mu = %.3f, alpha = %.2f, tau = ' ...
'%.1f'], omega_opt, mu_min, alpha, tau);

figure
subplot(211)
plot(omega,mu,'LineWidth',2);
xlabel('\omega'); ylabel('\mu')
grid on
title(str_title);
subplot(212)
plot(omega, theta1_opt, omega, theta2_opt,'LineWidth',2);
legend('\theta_{1opt}', '\theta_{2opt}')
xlabel('\omega'); ylabel('optimum \theta')
grid on
%
%   \mu for Gauss-Seidel relaxation
%
Ggs = zeros(nth1,nth2);
exp1=complex(cos(theta1),sin(theta1));
exp2=complex(cos(theta2),sin(theta2));
c=2*(1+alpha^2);
ctau=complex(alpha^2, (0.5*alpha*tau).*sint1);
for i1=1:nth1
    for i2=1:nth2
        if or(abs(theta1(i1))>= pi/2, abs(theta2(i2)) >= pi/2);
            Num = exp1(i1) + ctau(i1)*exp2(i2);
            Ggs(i1,i2) = abs( Num/(c-conj(Num)) );
        end
    end
end
[gmax,imax]=max(Ggs);
[mugs,jmax]=max(gmax);
theta1_gs_opt = theta1(imax(jmax));
theta2_gs_opt = theta2(jmax);
fprintf('alpha = %.2f, tau = %.1f, Ggs = %.4f, theta1 = %.4f, theta2 = %.4f\n', alpha, tau, ...
        mugs, theta1_gs_opt, theta2_gs_opt);
subplot(211)
hold on
plot(omega,repmat(mugs,1,length(omega)),'r--','LineWidth',2)