function [Q, ValA, ValAi, ValBi, ValB1,ValB2] = tt(g)

global T Q1 Q2 mu counter z G

counter = 0;

mu = 0.5; % k_F = 1
%
%fid = fopen('ivardataqc.txt', 'w');

G = 1i*g;
T = .1;
        
L = 1 + 5*T;
        
L = 2*L;


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

ValA=zeros(NQ,1);
ValAi=zeros(NQ,1);
ValBi=zeros(NQ,1);
ValB1=zeros(NQ,1);
ValB2=zeros(NQ,1);


for iQ = 1:NQ
    
    iQ
    
    Q1= Q(iQ)*[1,0];
    Q2= Q(iQ)*[1,0];
    
    %         tic
    %         for z = zmi:dz:zma
    z = 0; % 2D case
    
%    ValA(iQ) = ValA(iQ) + integral2(@integrand2A,-L-2*Q(iQ),L+2*Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5,'Method', 'iterated');
%    ValB(iQ) = ValB(iQ) + integral2(@integrand2B,-L-2*Q(iQ),L+2*Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5, 'Method', 'iterated');
%    ValBi(iQ) = ValBi(iQ) + intim(@integrand2B,-L-2*Q(iQ),L+2*Q(iQ), -L-Q(iQ), L+Q(iQ));
%    ValB1(iQ) = ValB1(iQ) + integral2(@integrand2B1,-L-2*Q(iQ),L+2*Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5, 'Method', 'iterated');
%    ValB2(iQ) = ValB2(iQ) + integral2(@integrand2B2,-L-2*Q(iQ),L+2*Q(iQ), -L-Q(iQ), L+Q(iQ), 'AbsTol', 1e-5, 'Method', 'iterated');

%     ValA(iQ) = ValA(iQ) + integral2(@integrand2A,-inf, inf, -inf, inf);
%     ValB1(iQ) = ValB1(iQ) + integral2(@integrand2B1,-inf, inf, -inf, inf);
%     ValB2(iQ) = ValB2(iQ) + integral2(@integrand2B2,-inf, inf, -inf, inf);

       disp('at zero angle')
%       ValA(iQ) = intim2x(@integrand2A,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ));
%        ValAi(iQ) = integral2(@integrand2A,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-4);
        ValBi(iQ) = integral2(@integrand2B,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-4);
        ValB1(iQ) = intim2x(@integrand2B,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ));
        ValB2(iQ) = intim2x(@integrand2B2,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ));

end

%fclose(fid)

end

%%%%%%% 2D fcn %%%%%
function f = integrand21(x,y) % 1,-1,2,-2 

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;



% %% w FF derivatives
% f = - 1./(1+ Exp(E1/T)).*(real(1./(E1-E2+G)).^2.*real(1./(E1-E4+G)) + real(1./(E1-E2+G)).*real(1./(E1-E4+G)).^2)...
%     - Exp(E1/T)/T./(1+ Exp(E1/T)).^2.*real(1./(E1-E2+G)).*real(1./(E1-E4+G))...
%     + 1./(1+ Exp(E2/T)).*real(1./(E2-E1+G)).^2.*real(1./(E2-E4+G))...
%     + 1./(1+ Exp(E4/T)).*real(1./(E4-E1+G)).^2.*real(1./(E4-E2+G));

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


end

%%%%%%% 2D fcn %%%%%
function f = integrand2A(x,y) % 1,2,-2,-1 

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;

%regularization 1
f = 2*nf(E1/T).*real(1./(E1-E2+G)).^2.*real(1./(E1-E3+G)) ...
    - nf(E2/T).*(real(1./(E2-E1+G)).*real(1./(E2-E3+G)).^2 + real(1./(E2-E1+G)).^2.*real(1./(E2-E3+G)))...
    - Exp(E2/T)/T.*nf(E2/T).^2.*real(1./(E2-E1+G)).*real(1./(E2-E3+G));

%regularization 2
% f = 2*nf(E1/T).*real(1./(E1-E2+G).^2).*real(1./(E1-E3+G)) ...
%     - nf(E2/T).*(real(1./(E2-E1+G)).*real(1./(E2-E3+G).^2) + real(1./(E2-E1+G).^2).*real(1./(E2-E3+G)))...
%     - Exp(E2/T)/T.*nf(E2/T).^2.*real(1./(E2-E1+G)).*real(1./(E2-E3+G));


end


%%%%%%% 2D fcn %%%%%
function f = integrand2B(x,y) % 1,-1,2,-2 

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;

% regularization 1
%f = (-4*nf(E1/T) - 2*(E1-E2).*Exp(E1/T)/T.*nf(E1/T).^2).*real(1./(E1-E2+G)).^3;

% regularization 2
f = (-2*nf(E1/T) + 2*nf(E2/T) - (E1-E2).*(Exp(E1/T)/T.*nf(E1/T).^2+ Exp(E2/T)/T.*nf(E2/T).^2)).*real(1./(E1-E2+G)).^3;


end



%%%%%%% 2D fcn %%%%%
function f = integrand2B1(x,y) % 1,-1,2,-2 

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;

% regularization 1
f = -4*nf(E1/T).*real(1./(E1-E2+G)).^3;

% regularization 2
%f = -4*nf(E1/T).*real(1./(E1-E2+G).^3);


end
function f = integrand2B2(x,y) % 1,-1,2,-2 

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;

% regularization 1
f =  - 2*Exp(E1/T)/T.*nf(E1/T).^2.*(E1-E2).*real(1./(E1-E2+G)).^3;
%f =  - 2*Exp(E1/T)/T.*nf(E1/T).^2.*real(1./(E1-E2+G)).^2;

% regularization 2
%f =  - 2*Exp(E1/T)/T.*nf(E1/T).^2.*real(1./(E1-E2+G).^2);


end

function nf_ = nf(E)

nf_ = 1./(1+Exp(E));

end

function int_ = intim(fun, xmin, xmax, ymin, ymax)

N = 2000;

[X,Y] = meshgrid(xmin:(xmax - xmin)/N: xmax, ymin: (ymax - ymin)/N: ymax);


int_= sum(sum(fun(X,Y)))*(xmax - xmin)/N*(ymax - ymin)/N;

end

