function [Q, alpha, Val1, Val2, Val3, Val4, Val5, Val6] = boxd3d_new(g)

global T Q1 Q2 mu z G



mu = 0.5; % k_F = 1

G = 1i*g;

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


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

Na = length(alpha);
%Na = 0

Q = 1.8:.025:2.1;
%Q = 1.8
%Q = 1/tan(pi/5)
pause(1)
%Q= 1.;

NQ = length(Q);

Val1=zeros(NQ,Na);
Val2=zeros(NQ,Na);
Val3=zeros(NQ,Na);
Val4 = zeros(NQ,1); % zero angle result
Val5 = zeros(NQ,1); % zero T analytical result
Val6 = zeros(NQ,1); % q-dependent non-coplanar vertex for FCC
for iQ = 1:NQ 
    Q(iQ)
    %zero angle calculation
    Q1= Q(iQ)*[1,0];
    Q2= Q(iQ)*[1,0];
% %     Val4(iQ) = 2*integral3(@integrand3A,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),-L, L,'AbsTol', 1e-5)...
% %         +1*integral3(@integrand3B,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),-L, L, 'AbsTol', 1e-5);
    
     Val6(iQ) = integral3(@integrand4,-L-1*Q(iQ),L+1*Q(iQ), -L-2*Q(iQ), L+0*Q(iQ),-L-1*Q(iQ), L,'AbsTol', 1e-5);
%     if Q(iQ) > 2.3
% %         Val5(iQ) = 16*pi/Q(iQ)^3*(2 - Q(iQ)^2)/(Q(iQ)^2 - 4)^(1.5);
% %     else
% %         Val5(iQ) = 0; % just set to 0 to avoid plotting divergent stuff
% %     end
% %     
%     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) = integral3(@integrand31,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),-L, L,'AbsTol', 1e-5);
%         Val2(iQ, ia) = integral3(@integrand32,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),-L, L,'AbsTol', 1e-5);
%         Val3(iQ, ia) = integral3(@integrand33,-L-2*Q(iQ),L+Q(iQ), -L-Q(iQ), L+Q(iQ),-L, L,'AbsTol', 1e-5);
%         
%     end
end

end

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

global T Q1 Q2 mu 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 = integrand32(x,y,z) % 1,-1,-2,2 

global T Q1 Q2 mu 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 = integrand33(x,y,z)  % 1,2,-1,-2

global T Q1 Q2 mu 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


%% %%%% 2D fcn %%%%%
%% noncoplanar diagram for FCC
function f = integrand4(x,y,z)  % 1,2,-1,-2

global T Q1 Q2 mu G


q0 = sqrt(Q1(1)^2 + Q1(2)^2);

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

q1 = q0/sqrt(3)*[1,1,1];
q2 = q0/sqrt(3)*[-1,1,-1];
q3 = q0/sqrt(3)*[-1,-1,1];


E1 = (x.^2 + y.^2 + z.^2)/2 - mu;
E2 = ((x + q1(1)).^2 + (y+ q1(2)).^2 + (z+q1(3)).^2)/2 - mu;
E3 = ((x + q1(1)+ q2(1)).^2 + (y+ q1(2)+ q2(2)).^2 + (z+q1(3)+q2(3)).^2)/2 - mu;
E4 = ((x + q1(1)+ q2(1)+ q3(1)).^2 + (y+ q1(2)+ q2(2)+ q3(2)).^2 + (z+q1(3)+q2(3)+q3(3)).^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 = integrand3A(x,y,z) % 1,2,-2,-1 

global T Q1 Q2 mu 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 = integrand3B(x,y,z) % 1,-1,2,-2 

global T Q1 Q2 mu 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