function [CM_Periodic_All,CM_Perturbed_All,RHS_All] = SystemSourceAbove(K,PML,PointSource,alpha_nodes,mesh,curve)

ng = mesh.numberpoints;
alpha_number = length(alpha_nodes);
%% Compute quadrature nodes and weights on each triangles:

[QuadTriangle] = QuadraturePointsTriangles(mesh);


%% Compute the coefficient and evaluation of the piecewise linear basis functions (PWL):

[PWL]= PWL_function(mesh,QuadTriangle);


%% Compute the PML function, i.e., S(x2):

S = PML_function(QuadTriangle,PML);


%% Construct the periodic & perturbed matrices:

[CM_Periodic_All] = CoeffMatrixPeriodic(K,S,PWL,alpha_nodes,mesh,QuadTriangle);

[CM_Perturbed_All]= B_Perturbed(K,S,PML,PWL,alpha_nodes,mesh,curve,QuadTriangle);

  
%% Define some terms 
PointSource_x = PointSource(1);
PointSource_y = PointSource(2);
Lowerbnd_indx = mesh.bnd.Lower;
Lowerbnd_xnode = mesh.points(1,mesh.bnd.Lower); %x1
Lowerbnd_ynode = mesh.points(2,mesh.bnd.Lower); %x2

Term01 = (Lowerbnd_xnode - PointSource(1)).^2;
Term02 = curve.psi(Lowerbnd_xnode,Lowerbnd_ynode);

Term1 = ( Term02 - PointSource(2) ).^2;
Term2 = ( Term02 + PointSource(2) ).^2;
Term3 = ( Lowerbnd_ynode - PointSource(2) ).^2;
Term4 = ( Lowerbnd_ynode + PointSource(2) ).^2;

main_Term = ( besselh(0,1,K * sqrt( Term01 + Term1 ))...
            - besselh(0,1,K * sqrt( Term01 + Term2 ))...
            - besselh(0,1,K * sqrt( Term01 + Term3 )) ...
            + besselh(0,1,K * sqrt( Term01 + Term4 )));

%% Construct the right-hand side for each alpha
RHS_All = cell(1,alpha_number);
parfor j = 1:alpha_number
    
    % === Compute the periodic FB transform of the Greens function for the lower boundary
    FB_Greensfunction = zeros(1,length(Lowerbnd_xnode));
    
    for l = -20:20
        
        Bc = sqrt(K^2 - (l + alpha_nodes(j)).^2);
        
        argument = angle(Bc);
        
        % Consider the branch cuts for the square root function
        if argument < -pi/4 || argument > 0.75 * pi
            Bc = - Bc;
        end
        
        FB_Greensfunction = FB_Greensfunction - (1/(2*pi))...
            .* exp( 1i * l .* Lowerbnd_xnode )...
            .* exp(-1i * (l + alpha_nodes(j)).* PointSource_x)...
            .* exp( 1i * Bc * PointSource_y)...
            .* Lowerbnd_ynode .* sinc(Bc .* Lowerbnd_ynode );
        
    end
    
     FB_Greensfunction = FB_Greensfunction ...
                       - exp( -1i * alpha_nodes(j).* Lowerbnd_xnode)...
                       * 1i .* main_Term./4;


    RHS = zeros(ng,1);
    
    RHS(Lowerbnd_indx) =  FB_Greensfunction.';

    RHS_All{j} = RHS;
    
end
