function [V,b] = MHP_elmgZareni(Uzly,Trojuhelniky,kappa,Uin,phi)

k = kappa*[sin(phi);0;-cos(phi)];

m = size(Trojuhelniky,2);
V = zeros(m,m);
b = zeros(m,1);

% Gaussova kvadratura
N = 3;
w = [8;5;5]/18;
ksi = 1/2+[0;-sqrt(15);sqrt(15)]/10;
% N = 5;
% c1 = 128/225; c2 = (322+13*sqrt(70))/900; c3 = (322-13*sqrt(70))/900;
% w = [c1;c2;c2;c3;c3]/2;
% C1 = 0; C2 = sqrt(5-2*sqrt(10/7))/3; C3 = sqrt(5+2*sqrt(10/7))/3;
% ksi = 0.5+[C1;-C2;C2;-C3;C3]/2;

for l=1:m
    
    Uzlyl = Uzly(:,Trojuhelniky(:,l));
    xl = sum(Uzlyl,2)/3;
    
    for j=1:m
        
        Uzlyj = Uzly(:,Trojuhelniky(:,j));
        Rj = [ Uzlyj(:,2)-Uzlyj(:,1) , Uzlyj(:,3)-Uzlyj(:,2) ];
        Jakobian = norm( [ det(Rj([2,3],:)) ; -det(Rj([1,3],:)); det(Rj([1,2],:)) ] );
        xj1 = Uzlyj(:,1);
        
        if l==j
            
            Rj1 = [ Uzlyj(:,1)-xl , Uzlyj(:,2)-Uzlyj(:,1) ];
            Jakobian1 = norm( [ det(Rj1([2,3],:)) ; -det(Rj1([1,3],:)); det(Rj1([1,2],:)) ] );
            Rj2 = [ Uzlyj(:,2)-xl , Uzlyj(:,3)-Uzlyj(:,2) ];
            Jakobian2 = norm( [ det(Rj2([2,3],:)) ; -det(Rj2([1,3],:)); det(Rj2([1,2],:)) ] );
            Rj3 = [ Uzlyj(:,3)-xl , Uzlyj(:,1)-Uzlyj(:,3) ];
            Jakobian3 = norm( [ det(Rj3([2,3],:)) ; -det(Rj3([1,3],:)); det(Rj3([1,2],:)) ] );
            
            for alfa=1:N
                for beta=1:N
                    norma = norm(xl-Rj*[ksi(alfa);ksi(alfa)*ksi(beta)]-xj1);
                    
                    % regularni cast
                    V(j,j) = V(j,j) + w(alfa)*w(beta)*ksi(alfa) * Jakobian * (exp(i*kappa*norma)-1) / (4*pi*norma);
                    
                    % singularni cast
                    norma = norm(-Rj1*[1;ksi(beta)]);
                    V(j,j) = V(j,j) + w(alfa)*w(beta) * Jakobian1 / (4*pi*norma);
                    norma = norm(-Rj2*[1;ksi(beta)]);
                    V(j,j) = V(j,j) + w(alfa)*w(beta) * Jakobian2 / (4*pi*norma);
                    norma = norm(-Rj3*[1;ksi(beta)]);
                    V(j,j) = V(j,j) + w(alfa)*w(beta) * Jakobian3 / (4*pi*norma);
                    
                end
            end
           
        else
            
            for alfa=1:N
                for beta=1:N
                    norma = norm(xl-Rj*[ksi(alfa);ksi(alfa)*ksi(beta)]-xj1);
                    V(l,j) = V(l,j) + w(alfa)*w(beta)*ksi(alfa) * Jakobian * exp(i*kappa*norma) / (4*pi*norma);
                end
            end
            
        end
        
    end
    
    
    b(l) = -Uin*exp(i*k'*xl);
    
end

