%% Solving non-periodic scattering problems in a domain with a local perturbation:
% incident field: Green's function with the point source below the surface

%% Note that input has special format:
% 'mesh size' for example '01' for meshsize 0.01
% source = [x_coordinate, y_coordinate]

function [Errors] = MainPML(meshsize,source,K)

format long e

%% parameters:
kappa = K-floor(K);

PML.H = 2.5;
PML.Lambda = 1.5;
PML.chi = exp(1i*pi/4);

rho = 2:2:20;

number_alpha = 20;  

%% Open and Read the generated mesh by the Gmsh in the x-space:

file = strcat('Mesh_', meshsize, '_test2.msh');
mesh = Read_Mesh(file);

Index_below_H = find( mesh.points(2,:) <= PML.H );

%% Quadrature nodes and weights on alpha in [-1/2,1/2] with the special path:

alpha = Split_Integration(number_alpha,kappa);

%% Define the periodic and locally perturbed curves and their 1st order derivatives:

curve.zeta = @(x) 1.5 + sin(x)/3 - cos(2*x)/4;  % periodic lower boundary

curve.diffzeta = @(x) cos(x)/3 + sin(2*x)/2;    % local perturbation on the lower boundary

curve.zetap_zeta = @(x) exp( ( (x .* (-1 < x  & x < 1)).^2 - 1 ).^(-1) ) ...
    .* sin( pi * (x + 1) ) .* (-1 < x  & x < 1);

curve.diffzetap_zeta = @(x) (pi * cos(pi * (x + 1)) - 2 * x .* sin(pi*(x + 1)).* (x.^2 - 1).^(-2)) ...
    .* exp( ((x .* (-1 < x  & x < 1)).^2 - 1).^(-1) ).* (-1 < x  & x < 1);


%%  Transform the periodic structure to the locally perturbed periodic structure

curve.psi = @(x,y) y + ( (y - PML.H).^3 .* (curve.zeta(x) - PML.H).^(-3) ) ...
    .* curve.zetap_zeta(x) .* (y < PML.H);

%% Compute the exact scattered field & Periodic transformed field:

[ScatteredField.exact] = Greensfunction(mesh,K,source,curve);

ScatteredField.norm = norm( ScatteredField.exact(Index_below_H),2 );


for index= 1:length(rho)
    tic
    PML.rho = rho(index);

    %% Construct the Components of the linear system:
    tic
    [CM_Periodic,CM_Perturbed,F] = System(K,PML,source,alpha.nodes,mesh,curve);

    %% Use the Schur complement method to solve the linear system:

    ScatteredField.num = Solve_System_Schurcomplement(CM_Periodic,CM_Perturbed,F,mesh,alpha);
    fprintf('Computational time for solving the system with schur complement~= %.6f (s).\n', toc);


    fprintf('Computational time for one alpha ~= %.6f (s).\n', toc);

    Error.Absolute = ScatteredField.exact(Index_below_H) ...
        - ScatteredField.num(Index_below_H).' ;

    Error.Relative = Error.Absolute/ScatteredField.norm;

    RelativeError = norm(Error.Absolute,2)/ScatteredField.norm;

    Errors(index)= RelativeError


    %% uncomment this part for plotting:
    % Kstring= strcat('_k',char(string(K*10)));
    % file_solution = fopen(strcat('Solution', meshsize,'paper',Kstring),'w');

    % for j = 1:length(mesh.index_intpoints)
    %     fprintf(file_solution, '%d %e %e\n',mesh.index_intpoints(j),...
    %     real(ScatteredField.num(j)),imag(ScatteredField.num(j)));
    % end

    % file_error = fopen(strcat('Error', meshsize,'paper',Kstring),'w');

    % for j = 1:length(Index_below_H)
    %     fprintf(file_error, '%d %e %e\n',mesh.index_intpoints(Index_below_H(j)),...
    %     real(Error.Absolute(j)),imag(Error.Absolute(j)));
    % end


end

end

