import math
from scipy.optimize import fsolve
from sympy import nsolve

eRs=1.00918653134839
eDs2=0.00166534843904
eDs3=-4.03945520938235*10**-5
eDs4=1.17223093696419*10**-5
eRb=1.00576330135629
eDb2=0.000208948562128
eDb3=-2.36794450340036*10**-6
eDb4=3.25268397646453*10**-7

CVsb=0.000175588301585
CVs2b=-6.26860807093143*10**-6
CVsb2=-4.17120782810223*10**-6
CVs2b2=6.65898977814507*10**-7
CVs3b=1.89847158043424*10**-6
CVsb3=4.29542210089401*10**-7




Rf=1
def dv1(wS,wB,gamma,alpha,Rf):
    x=Rf+wS*eRs+wB*eRb+gamma/alpha
    x*=100
    u=eRs*(x)
    u+=eRs*gamma*(1+gamma)*(x**(-gamma-2))*(1/2*wS**2*eDs2 + wS*wB*CVsb+1/2*wB**2*eDb2)
    u+=(wS*eDs2+wB*CVsb)*(-gamma)*(x**(-gamma-1))
    u+=-eRs*gamma*(1+gamma)*(2+gamma)*(x**(-gamma-3))*(1/6*wS**3 +1/2*wS**2*wB*CVs2b + 1/2*wS*wB**2*CVsb2+1/6*wB**3*eDb3)
    u+=gamma*(1+gamma)*(x**(-gamma-2))*(1/2*wS**2*eDs3+wS*wB*CVs2b+1/2*wB**2*CVsb2)
    u+=eRs*gamma*(1+gamma)*(2+gamma)*(3+gamma)*(x**(-gamma-4))*(1/24*wS**4*eDs4+1/6*wS**3*wB+1/4*wS**2*wB**2*CVs2b2 + 1/6*wS*wB**2*CVsb3 + 1/24*wB**4*eDb4)
    u+=-gamma*(1+gamma)*(2+gamma)*(x**(-gamma-3))*(1/6*wS**3*eDs4 + 1/2*wS**2*wB*CVs3b+1/2*wS*wB**2*eDs2*eDb2 +CVs2b2+1/6*wB**3*CVsb3)
    return u*gamma

def dv2(wS,wB,gamma,alpha,Rf):
    x=Rf+wS*eRs+wB*eRb+gamma/alpha
    x*=100
    u=eRb*(x)
    u+=eRb*gamma*(1+gamma)*(x**(-gamma-2))*(1/2*wB**2*eDb2 + wS*wB*CVsb+1/2*wS**2*eDs2)
    u+=(wB*eDb2+wS*CVsb)*(-gamma)*(x**(-gamma-1))
    u+=-eRb*gamma*(1+gamma)*(2+gamma)*(x**(-gamma-3))*(1/6*wB**3 +1/2*wB**2*wS*CVsb2 + 1/2*wB*wS**2*CVs2b+1/6*wS**3*eDs3)
    u+=gamma*(1+gamma)*(x**(-gamma-2))*(1/2*wB**2*eDb3+wS*wB*CVsb2+1/2*wS**2*CVs2b)
    u+=eRb*gamma*(1+gamma)*(2+gamma)*(3+gamma)*(x**(-gamma-4))*(1/24*wB**4*eDb4+1/6*wS**3*wB+1/4*wS**2*wB**2*CVs2b2 + 1/6*wB*wS**2*CVs3b + 1/24*wS**4*eDs4)
    u+=-gamma*(1+gamma)*(2+gamma)*(x**(-gamma-3))*(1/6*wB**3*eDb4 + 1/2*wB**2*wS*CVsb3+1/2*wB*wS**2*eDs2*eDb2 +CVs2b2+1/6*wS**3*CVs3b)

    return u*gamma


def equations(x,y):
    wS, wB = x,y
    wS/=100
    wB/=100
    # print(wS,wB)
    return (dv1(wS,wB,gamma,alpha,Rf), 1-wS-wB,dv2(wS,wB,gamma,alpha,Rf))

nsolve(equations(x,y),[x,y],[50,50])