function [L_IW] = TSbox_interweave(LECC, LOCS, LOSC, LESS, harmo, basis)
% TSbox_interweave combines matrices LYXX into an interweaved discretization L_IW

% harmo (string, 'odd_harmo' or 'even_harmo')
% describes the parity of the harmonics in the operator. For instance:
%     L = cos theta + sin theta     has parity 1  (odd harmonics only)
%     L = cos^2 theta + d_theta     has parity 0 (even harmonics only)
% ..........................................................................
% basis (string, 'a' or 'b')
% designates which basis is used for the expansion, as defined in Miquel & Julien, submitted to JCP:
% 'a' corresponds to equation (18a): (1, sin 2t, cos 2t, sin 4t, cos 4t, ...)
% 'b' corresponds to equation (18b): (cos t, sin t, cos 3t, sin3t, cos 5t, sin 5t, ...)
% ..........................................................................
% The basis for projection is automatically deduced from the variables harmo and basis.
% Hence, matrices LXC and LXS are interwoven correctly. For (harmo, basis) set to:
%  .  (odd_harmo, a) the projection is sparse on b
%  . (even_harmo, a) the projection is sparse on a
%  .  (odd_harmo, b) the projection is sparse on a
%  . (even_harmo, b) the projection is sparse on b
% ..........................................................................

% check the size of the matrices:
ntotal = size(LECC,1);
L_IW=sparse(ntotal, ntotal);

if strcmp(harmo,'odd_harmo')
    if strcmp(basis,'a')
    % variables are expanded on (18a), equations are projected on (18b)
    % indices for the expansion:
    icos_expan = (1:2:ntotal); 
    isin_expan = (2:2:ntotal); 
    icos_proj  = (2:2:ntotal); 
    isin_proj  = (1:2:ntotal); 
    elseif strcmp(basis,'b')
    % variables are expanded on (18b), equations are projected on (18a)
    % indices for the expansion:
    icos_expan = (2:2:ntotal); 
    isin_expan = (1:2:ntotal); 
    icos_proj  = (1:2:ntotal); 
    isin_proj  = (2:2:ntotal); 
    else 
    disp('wrong ''basis'' flag. Must be ''a'' or ''b''');
    return;
    end;
elseif strcmp(harmo,'even_harmo')
    if strcmp(basis,'a')
    % variables are expanded and projected (18a)
    % indices for the expansion:
    icos_expan = (1:2:ntotal); 
    isin_expan = (2:2:ntotal); 
    icos_proj  = icos_expan;
    isin_proj  = isin_expan;
    elseif strcmp(basis,'b')
    % variables are expanded and projected (18b)
    % indices for the expansion:
    icos_expan = (2:2:ntotal); 
    isin_expan = (1:2:ntotal); 
    icos_proj  = icos_expan;
    isin_proj  = isin_expan;
    else 
    disp('wrong ''basis'' flag. Must be ''a'' or ''b''');
    return;
    end;
else
disp('wrong ''harmo'' flag. Must be ''odd_harmo'' or ''even_harmo''');
return;
end;

iodd  = (1:2:ntotal); 
ieven = (2:2:ntotal); 

L_IW (iodd,  iodd)  = LECC (icos_proj, icos_expan); 
L_IW (ieven, ieven) = LESS (isin_proj, isin_expan); 
L_IW (iodd,  ieven) = LOCS (icos_proj, isin_expan); 
L_IW (ieven, iodd)  = LOSC (isin_proj, icos_expan); 

end

