function [Q, alpha, Val1, Val2, Val3, Val4] = boxd2d_new(g)

global T Q1 Q2 mu z G



mu = 0.5; % k_F = 1

G = 1i*g;

z = 0; 
T = .01;
        
L = 1 + 5*T;
        
L = 1.5*L;


alpha = 0.001:pi/64:pi/2 + 0.001
% alpha = .5;

Na = length(alpha);

Q = 1.5:.01:2.5;
pause(1)
%Q= 1.;

NQ = length(Q);

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

for iQ = 1:NQ 
    Q(iQ)
    %zero angle calculation
    Q1= Q(iQ)*[1,0];
    Q2= Q(iQ)*[1,0];
    Val4(iQ) = 2*integral2(@integrand2A,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5)...
        +1*integral2(@integrand2B,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5);
    for ia = 1:Na
        disp('**** new iter ****')
        Q1= Q(iQ)*[cos(alpha(ia)), sin(alpha(ia))]
        Q2= Q(iQ)*[cos(alpha(ia)), -sin(alpha(ia))]
        
%         [X,Y] = meshgrid(-L-2*Q(iQ):.2:L+Q(iQ));
%         [X,Y] = meshgrid(-2:.01:2);
%         
%         figure(1); surf(X,Y, integrand23(X,Y))
%         
%         pause
        
        Val1(iQ, ia) = integral2(@integrand21,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5);
        Val2(iQ, ia) = integral2(@integrand22,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5);
        Val3(iQ, ia) = integral2(@integrand23,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),'AbsTol', 1e-5);

        
        
        
        %        Val3(iQ, ia) = intim2x(@integrand22,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ));
        
        
        
        
    end
end

end

%% regularization relative to zero angle
%% 
function f = integrand21(x,y) % 1,-1,2,-2 

global T Q1 Q2 mu z G

x = x + 1e-7;
y = y + 1e-9;

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)).^2.*real(1./(E1-E4+G)) + real(1./(E1-E2+G)).*real(1./(E1-E4+G)).^2)...
    + nfp(E1/T)/T.*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));


end


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

global T Q1 Q2 mu z G

x = x + 1e-7;
y = y + 1e-9;

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)).^2.*real(1./(E1-E4+G)) + real(1./(E1-E2+G)).*real(1./(E1-E4+G)).^2)...
    + nfp(E1/T)/T.*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));

end



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

global T Q1 Q2 mu z G

x = x + 1e-7;
y = y + 1e-9;

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-E3+G)).*real(1./(E1-E4+G)) ...
    + 1./(1+ Exp(E2/T)).*real(1./(E2-E1+G)).*real(1./(E2-E3+G)).*real(1./(E2-E4+G))...
    + 1./(1+ Exp(E3/T)).*real(1./(E3-E1+G)).*real(1./(E3-E2+G)).*real(1./(E3-E4+G))...
    + 1./(1+ Exp(E4/T)).*real(1./(E4-E1+G)).*real(1./(E4-E2+G)).*real(1./(E4-E3+G));




end

%% zeroo-anlge functions

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

global T Q1 Q2 mu  z G

x = x + 1e-7;
y = y + 1e-9;

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 = nf(E1/T).*real(1./(E1-E2+G)).^2.*real(1./(E1-E3+G)) ...
    +nf(E3/T).*real(1./(E3-E2+G)).^2.*real(1./(E3-E1+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  z G

x = x + 1e-7;
y = y + 1e-9;

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 = (-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

%%
function nf_ = nf(E)

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

end

function nfp_ = nfp(E)

nfp_ = -Exp(E)./(1+Exp(E)).^2;

end


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

disp('****INTIM****')

N = 3000;

[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;

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;

abserr = int__ - int_
relerr = abserr/int_
end