% TSbox_exemple_translated_harmonics_well.m  exemplifies how to
% use the TSbox to find the modes of a translated harmonic well
% See section  4.1, equation [40] in Miquel & Julien, J. Comput.
% Phys., 349C (2017).
% https://doi.org/10.1016/j.jcp.2017.08.034
% 
% Code by Ben Miquel, August 2017.

clear
close all

N    =64;   % set resolution.
Lmap = 3 ;      % set the mapping parameter L (see eq. [20]).
Nbuff= N+20;   % discretization of the operators of size N  are
               % obtained  as products  of matrices of  greater
               % size Nbuff, to avoid truncation errors (see Ap
               % pendix B).

% define elementary matrices
[SSC, SCS, CCC, CSS, DSC, DCS] = TSbox_SCD_matrices(Nbuff);

% #############################################################
% ... with these building blocks in hand, we proceed to construct
% the even and odd components of the M and L operators, following
% the prescription of appendix C.
% #############################################################

% discretization of the even component of the M and L operators
temp = SCS*SSC*SCS*SSC*SCS*SSC*DCS*DSC;
LECC = -temp/(Lmap^2);
temp = SCS*SSC*SCS*SSC*CCC*SCS*DSC;
LECC = LECC - 2*temp/(Lmap^2);
LECC = LECC + Lmap^2 *CCC*CCC;
%..............................................................
MECC = SCS*SSC;
%..............................................................
temp = SSC*SCS*SSC*SCS*SSC*SCS*DSC*DCS;
LESS = -temp/(Lmap^2);
temp = SSC*SCS*SSC*SCS*CSS*SSC*DCS;
LESS = LESS - 2*temp/(Lmap^2);
LESS = LESS + Lmap^2 *CSS*CSS;
%..............................................................
MESS = SSC*SCS;


% discretization of the odd  component of the M and L operators
LOSC = -6*Lmap*SSC*CCC;
%..............................................................
MOSC = sparse(Nbuff, Nbuff);  % M is even here, M_(odd) = 0
%..............................................................
LOCS = -6*Lmap*SCS*CSS;
%..............................................................
MOCS = sparse(Nbuff, Nbuff);  % M is even here, M_(odd) = 0

% We interweave the matrices using basis [18a], hence the flag 'a'.
% The mapped and regularized equation [41] contains exclusively even
% harmonics, hence the flag 'even_harmo'.
% Alternatively, we could have chosen a discretisation using basis 
% [18b], by setting the last argument to 'b'.
% See the description at the top of TSbox_interweave for more details.

L = TSbox_interweave(LECC, LOCS, LOSC, LESS, 'even_harmo', 'a');
M = TSbox_interweave(MECC, MOCS, MOSC, MESS, 'even_harmo', 'a');

% truncate
L(N+1:end,:) = [];
M(N+1:end,:) = [];
L(:,N+1:end) = [];
M(:,N+1:end) = [];
% dispose of unnecessary variables
%clear *CS *SS *CC *SC temp Nbuff

% solve partially the eigenproblem
[X,V] = eigs(full(L),full(M),3,4);  %  find the 3 eigenvalues
                                    %  closest to 4. 

% #############################################################
% plot the mode with the closest eigenvalue to our guess
spec_coefs = X(:,1);
[psi,gtheta] = TS_itransform(spec_coefs,3*N,'a');
psi = psi / max(abs(psi(:))); % normalization

figure(1)
subplot(2,2,1)
spy(L)
title('matrix L')
subplot(2,2,2)
spy(M)
title('matrix M')
subplot(2,2,3)
semilogy(abs(spec_coefs))
axis([0 N 1e-15 1])
title('spectrum')
subplot(2,2,4)
plot(Lmap*cot(gtheta),psi,'o-')
axis([-4 10 -1.2 1.2 ])
title('eigenmode')
