from scipy.optimize import fsolve
from scipy.integrate import dblquad
from numpy import sqrt,cos,pi,absolute

Ueh=2320.5
Uhh=2192.5
ED=120
LCP=-59.28179796
UCP=-58.96721714


def Ih(mu,x):
    var,err= dblquad(lambda ky, kx: 1/sqrt((ED-mu-2540*(2-cos(kx)-cos(ky)))**2+x**2)/8/pi/pi, -pi, pi, lambda kx:-pi,lambda kx: pi)
    return var

def Ie(mu,x):
    var,err= dblquad(lambda ky, kx: 1/sqrt((2540*(2-cos(kx)-cos(ky))-mu)**2+x**2)/8/pi/pi, -pi, pi, lambda kx:-pi,lambda kx: pi)
    return var

def spp(mu,p):
    x, y = p
    return (absolute(x)-2*Ueh*absolute(y)*Ih(mu,y),4*Ueh**2*Ie(mu,x)*Ih(mu,y)/(1+Uhh*Ih(mu,y))-1)

def spm(mu,y):
    return Uhh*Ih(mu,y)-1

def sis(mu,p):
    x,y=p
    return(Uhh*Ih(mu,y)-1,  Ie(mu,x)*2*Ueh**2/Uhh-1)


def solution(mu):
    if mu<=LCP:
        Dele=0
        Delh=fsolve(lambda y:spm(mu,y),4)
        phase=0
    elif mu> LCP and mu<UCP:
        Dele,Delh=fsolve(lambda p: sis(mu,p),(4,4))
        phase=-Uhh/2/Ueh*Dele/Delh
    elif mu>=UCP:
        Dele,Delh=fsolve(lambda p: spp(mu,p),(9,4))
        phase=-1
    return Dele,Delh,phase

print(solution(-60))
(0, array([ 4.1349274]), 0)

print(solution(-59))
(8.2553379473241311, 4.1210736864119193, -0.94635157721569463)

print(solution(-58))
(8.9677408867771558, 4.2580061198120394, -1)