clc
clear

K = sqrt(10);
kappa = K-floor(K);

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


source = [0,0.5];

number_alpha = 20;    

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

meshsize = '01';
file = strcat('Mesh_', meshsize, '_test2.msh');
mesh = Read_Mesh(file);
index_below_H = find( mesh.points(2,:) <= PML.H );

Kstring= strcat('_k',char(string(K^2)));

%%
solutionFile = strcat('paperResult/Solution', meshsize,'paper',Kstring);
u0 = Read_Solution(solutionFile);

%% 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);

[CM_Periodic,CM_Perturbed,F] = System(K,PML,source,alpha.nodes,mesh,curve);

%% Extend the numerical solution to the previous and next cells:

cm1=183309;
c1=512567;

extendedIndices = [mesh.index_intpoints(index_below_H)+cm1,mesh.index_intpoints(index_below_H),mesh.index_intpoints(index_below_H)+c1];

[u1,um1,uu]= transformedFieldScatteredField(u0,CM_Periodic,CM_Perturbed,F,mesh.points(1,index_below_H),alpha,index_below_H);

extendedNumericalSolution=[um1.',u0(index_below_H).',u1.'];


file_SolutionExtended = fopen(strcat('paperResult/SolutionExtended', meshsize,'paper', Kstring),"w");
for j = 1:length(extendedIndices)
    fprintf(file_SolutionExtended, '%d %e %e\n',extendedIndices(j),...
        real(extendedNumericalSolution(j)),imag(extendedNumericalSolution(j)));
end


%% Extend the exact solution to the previous and next cells
extendedExactSolution = GreensfunctionExtended(mesh,K,source,curve,index_below_H);
l = length(index_below_H);

for idx = 0:2
    norm(extendedExactSolution(l*idx+1:l*(idx+1)),inf)
end

%% Compute errors
extendedError=extendedExactSolution - extendedNumericalSolution;

file_ErrorExtended = fopen(strcat('paperResult/ErrorExtended', meshsize,'paper',Kstring),"w");
for j = 1:length(extendedError)
    fprintf(file_ErrorExtended, '%d %e %e\n',extendedIndices(j),...
        real(extendedError(j)),imag(extendedError(j)));
end

