% PML

sigma_x_PML = zeros(N_p,K);
sigma_y_PML = zeros(N_p,K);

D_PML = 10;
sigma_m = 2;
beta = 2;

indexes_R = find(x > x_u + NODETOL);
indexes_L = find(x < x_l - NODETOL);
sigma_x_PML(indexes_R) = sigma_m * (abs(x(indexes_R) - x_u) / D_PML) .^ beta;
sigma_x_PML(indexes_L) = sigma_m * (abs(x(indexes_L) - x_l) / D_PML) .^ beta;

indexes_T = find(y > y_u + NODETOL);
indexes_B = find(y < y_l - NODETOL);
sigma_y_PML(indexes_T) = sigma_m * (abs(y(indexes_T) - y_u) / D_PML) .^ beta;
sigma_y_PML(indexes_B) = sigma_m * (abs(y(indexes_B) - y_l) / D_PML) .^ beta;

sigma_x_PML_plus_sigma_y_PML = sigma_x_PML + sigma_y_PML;
sigma_x_PML_times_sigma_y_PML = sigma_x_PML .* sigma_y_PML;

sigma_y_PML_times_n_x = sigma_y_PML(Fmask(:),:) .* n_x;
sigma_x_PML_times_n_y = sigma_x_PML(Fmask(:),:) .* n_y;


% PML access information
[rows1,cols1] = find(x > x_u + NODETOL);
[rows2,cols2] = find(x < x_l - NODETOL);
[rows3,cols3] = find(y > y_u + NODETOL);
[rows4,cols4] = find(y < y_l - NODETOL);
elements_PML = unique([cols1; cols2; cols3; cols4]);

elements_Euler = setdiff((1:K)',elements_PML);

Faces_PML = zeros(3,K);
Faces_PML(:,elements_PML) = ones(3,length(elements_PML));
Faces_PML = Faces_PML(:);
Nodes_PML = ones(N_fp,1) * Faces_PML';
Nodes_PML = Nodes_PML(:);
mapPML = find(Nodes_PML == 1); % faces nalezici elementum v PML
mapPML_M = unique([mapPML;mapP(mapPML)]); % body na faces na vsech elementech lezicich v PML nebo sousedicich s elementem v PML
mapPML_P = mapP(mapPML_M); % prilehajici body na body z mapPML_P
vmapPML = vmapM(mapPML);
vmapPML_M = vmapM(mapPML_M);
vmapPML_P = vmapM(mapPML_P);



% sponge zone

sigma_sponge = zeros(N_p,K);

D_sponge = 100;
sigma_m = 0.3;
beta = 2;

indexes_R = find((x > x_u) & (abs(y) < y_u + NODETOL));
sigma_sponge(indexes_R) = sigma_m * (abs(x(indexes_R) - x_u) / D_sponge) .^ beta;
indexes_L = find((x < x_l) & (abs(y) < y_u + NODETOL));
sigma_sponge(indexes_L) = sigma_m * (abs(x(indexes_L) - x_l) / D_sponge) .^ beta;
indexes_T = find((y > y_u) & (abs(x) < x_u + NODETOL));
sigma_sponge(indexes_T) = sigma_m * (abs(y(indexes_T) - y_u) / D_sponge) .^ beta;
indexes_B = find((y < y_l) & (abs(x) < x_u + NODETOL));
sigma_sponge(indexes_B) = sigma_m * (abs(y(indexes_B) - y_l) / D_sponge) .^ beta;
indexes_RT = find((x > x_u) & (y > y_u));
vectors = [x(indexes_RT) - x_u,y(indexes_RT) - y_u];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
sigma_sponge(indexes_RT) = sigma_m * (distances / D_sponge) .^ beta;
indexes_LT = find((x < x_l) & (y > y_u + NODETOL));
vectors = [x(indexes_LT) - x_l,y(indexes_LT) - y_u];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
sigma_sponge(indexes_LT) = sigma_m * (distances / D_sponge) .^ beta;
indexes_RB = find((x > x_u) & (y < y_l + NODETOL));
vectors = [x(indexes_RB) - x_u,y(indexes_RB) - y_l];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
sigma_sponge(indexes_RB) = sigma_m * (distances / D_sponge) .^ beta;
indexes_LB = find((x < x_l) & (y < y_l + NODETOL));
vectors = [x(indexes_LB) - x_l,y(indexes_LB) - y_l];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
sigma_sponge(indexes_LB) = sigma_m * (distances / D_sponge) .^ beta;

elements_sponge = elements_PML;