function [Q, alpha, Val1] = triad2p1_loop(N)

global T Q1 Q2 mu counter z G

counter = 0;

mu = 0.5; % k_F = 1


G = 0.0005i;
T = .005;


L = 1 + 5*T;
L = 1.5*L; % limits of integration. Result should not depend on wheter we take L -> A*L as long as A >1.

zmi = -L; zma = L; dz = (zma-zmi)/N;

%alpha = (.05):.05: (pi/2 - .05); 
%alpha = (.2):.01: (pi/2 - .2); 
%alpha = (.2):.01: (.2); 
alpha = pi/6; %makes equilateral triangle needed for 3rd order term
Na = length(alpha);

Q = 1:.02:2.5;
%Q= sqrt(2) -.1;
%Qqc = 2*cos(.5536);
%Q = Qqc
NQ = length(Q);

Val1=zeros(NQ,Na);
Val2=zeros(NQ,Na);
Val3=zeros(NQ,Na);

for ia = 1:Na; 
    alpha(ia)
    for iQ = 1:NQ
        Q(iQ)
        
        Q1= Q(iQ)*[cos(alpha(ia)), sin(alpha(ia))];
        Q2= Q(iQ)*[cos(alpha(ia)), -sin(alpha(ia))];

%         tic
%         for z = zmi:dz:zma
            z = 0; % 2D case
            
            Val1(iQ,ia) = Val1(iQ,ia) + quad2d(@integrand23,-L-Q(iQ),L, -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-4);
%         end
%         toc
        
%         Val1(iQ,ia) = Val1(iQ,ia)*dz;
    end
end

end




%%%%%% 2D fcn %%%%%
function f = integrand23(x,y)

global T Q1 Q2 mu counter z G

counter = counter + 1;

E1 = (x.^2 + y.^2 + z^2)/2 - mu;
E2 = ((x + Q1(1)).^2 + (y+ Q1(2)).^2 + z^2)/2 - mu;
E3 = ((x + Q1(1) + Q2(1)).^2 + (y+ Q1(2) + Q2(2)).^2 + z^2)/2 - mu;
E4 = ((x + Q2(1)).^2 + (y+ Q2(2)).^2 + z^2)/2 - mu;


f =   1./(1+ exp(E1/T)).*real(1./(E1-E2+G)).*real(1./(E1-E4+G)) ...
    + 1./(1+ exp(E2/T)).*real(1./(E2-E1+G)).*real(1./(E2-E4+G))...
    + 1./(1+ exp(E4/T)).*real(1./(E4-E1+G)).*real(1./(E4-E2+G));




end


