% Computes and saves the operator H_k acting on the Fourier basis 
% for a collection of grid points in the k-plane. 
% H_k is the following integral operator:
%
%   (H_k f)(z) = \int_{unitcircle} H1(k*(z-w))*f(w) dsigma(w), 
%
% where z is a point on the unit circle, H1(z) = G1(z)-G0(z),
% and G1 is Faddeev's Green function for the Laplacian.
% The function H1 can be evaluated using routine H1.m.
% 
% Order Ntrig of trigonometric approximation is defined in this routine.
% Collection of grid points in the k-plane is created by the routine 
% tDN1_kvec_comp.m, which should be run before this file.
%
% Samuli Siltanen May 2012

% Load the k-grid
load data/kvec kvec 

% Limit for up to how large |k| we will use the boundary integral equation
% approach for computing the scattering transform. That approach does not
% work for too large |k|
BIEindex = max(find(kvec<10));
save data/BIEindex BIEindex

% Load order of trigonometric approximation
load data/Ntrig Ntrig

% Construct integration points (angles) on the circle
Ntheta = 256;
theta  = 2*pi*[0:(Ntheta-1)]/Ntheta;
theta  = theta(:);
Dtheta = theta(2)-theta(1);
save data/theta theta Ntheta Dtheta

% Loop over points in the k-grid
for kkk = 1:BIEindex
    k = kvec(kkk);
    tic
    
    % Initialization of the result
    Nvec = [-Ntrig : Ntrig];
    Hk   = zeros(length(Nvec));
    
    % Loop over Fourier basis functions and apply the operator to each of them.
    % Then compute inner products of the result with basis functions.
    % This way we build a matrix for the operator H_k.
    for nnn = 1:length(Nvec)
        % Trigonometric basis function
        n    = Nvec(nnn); 
        bfun = 1/sqrt(2*pi)*exp(i*n*theta); 
        
        % Compute result of operator H_k applied to the basis function
        Hk_bfun = zeros(Ntheta,1);
        for iii = 1:Ntheta
            HH = H1(k*(exp(i*theta(iii)) - exp(i*theta)));
            Hk_bfun(iii) = Dtheta*sum(HH .* bfun);
        end
        % Expand the result in trigonometric basis. 
        for jjj = 1:length(Nvec)
            Hk(jjj,nnn) = 1/sqrt(2*pi)*Dtheta*exp(i*Nvec(jjj)*theta)'*Hk_bfun;
        end 
        disp(['Done ', num2str(nnn), ' out of ', num2str(length(Nvec))])
    end
    
    % Save the Hk matrix to file
    eval(['save data/Hk_', num2str(kkk), ' Hk'])
    
    % Monitor the run
    disp(['Computation ', num2str(kkk), ' of ', num2str(BIEindex),...
        ' took ', num2str(toc), ' seconds'])
end % for kkk


