T = .3;
mu = 0.5;
NF[E_] := 1/(1+Exp[E/T]);
Q = 1.8;
alpha = .001;
Q1 = Q {Cos[alpha], Sin[alpha]};
Q2 = Q {Cos[alpha], -Sin[alpha]};


En[px_,py_]:= ((px+.000001)^2+py^2)/2 -mu;
E1[px_,py_] :=En[px,py];
E2[px_,py_] :=En[px+Q1[[1]],py+Q1[[2]]];
E3[px_,py_] :=En[px+Q1[[1]]+Q2[[1]],py+Q1[[2]]+Q2[[2]]];
E4[px_,py_] :=En[px+Q2[[1]],py+Q2[[2]]];
E5[px_,py_] :=En[px-Q2[[1]],py-Q2[[2]]];




Omega21[px_,py_] = -NF[E1[px,py]](1/(E1[px,py] - E2[px,py])^2/(E1[px,py] - E4[px,py]) +1/(E1[px,py] - E2[px,py])/(E1[px,py] - E4[px,py])^2) +
NF'[E1[px,py]]/(E1[px,py] - E2[px,py])/(E1[px,py] - E4[px,py])+
NF[E2[px,py]]/(E2[px,py] - E1[px,py])^2/(E2[px,py]-E4[px,py]) + NF[E4[px,py]]/(E4[px,py]-E1[px,py])^2/(E4[px,py]-E2[px,py]);


Omega22[px_,py_] = -NF[E1[px,py]](1/(E1[px,py] - E2[px,py])^2/(E1[px,py] - E5[px,py]) +1/(E1[px,py] - E2[px,py])/(E1[px,py] - E5[px,py])^2) +
NF'[E1[px,py]]/(E1[px,py] - E2[px,py])/(E1[px,py] - E5[px,py])+
NF[E2[px,py]]/(E2[px,py] - E1[px,py])^2/(E2[px,py]-E5[px,py]) + NF[E5[px,py]]/(E5[px,py]-E1[px,py])^2/(E5[px,py]-E2[px,py]);

Omega23[px_,py_] = NF[E1[px,py]]/(E1[px,py] - E2[px,py])/(E1[px,py] - E3[px,py])/(E1[px,py] - E4[px,py]) +  NF[E2[px,py]]/(E2[px,py] - E1[px,py])/(E2[px,py] - E3[px,py])/(E2[px,py] - E4[px,py]) +
NF[E3[px,py]]/(E3[px,py] - E1[px,py])/(E3[px,py] - E2[px,py])/(E3[px,py] - E4[px,py]) +
NF[E4[px,py]]/(E4[px,py] - E1[px,py])/(E4[px,py] - E2[px,py])/(E4[px,py] - E3[px,py]) ;

E4[0,0]

Omega23[0.1,.003]

Plot3D[NF[E4[x,y]], {x,-2,2}, {y,-2,2}, PlotRange -> All]

Omega4A[px_,py_]:= ((NF[E1[px,py]])/((E1[px,py]-E2[px,py])^2(E1[px,py]-E3[px,py])))+((NF[E3[px,py]])/((E3[px,py]-E2[px,py])^2(E3[px,py]-E1[px,py])))+((NF'[E2[px,py]])/((E2[px,py]-E1[px,py])(E2[px,py]-E3[px,py])))-
((NF[E2[px,py]])/((E2[px,py]-E1[px,py])^2(E2[px,py]-E3[px,py])))-((NF[E2[px,py]])/((E2[px,py]-E1[px,py])(E2[px,py]-E3[px,py])^2))

Plot3D[Omega21[x,y], {x, -5,5}, {y, -5, 5},  PlotRange -> All, PlotPoints-> 100]

Omega2[px_,py_]:= Re[((NF[E1[px,py]]-NF[E2[px,py]])/((E1[px,py]-E2[px,py] + I 0.0)))]

Plot3D[Omega2[x,y], {x,-3,3},{y,-3,3}, PlotRange -> All, PlotPoints-> 40]

Omega4B[px_,py_]:= -2((NF[E1[px,py]]-NF[E2[px,py]])/((E1[px,py]-E2[px,py])^3))+((NF'[E1[px,py]]+NF'[E2[px,py]])/((E1[px,py]-E2[px,py])^2))

Plot3D[Omega4A[x,y], {x, -5,5}, {y, -5, 5}, PlotRange -> All, PlotPoints-> 40]