Gamma1 for a mix of ideal gas and radiation

    Hansen & Kawaler, page 177, Eqn. 3.110
    
    Parameters
    ----------
    beta : float
        Gas pressure fraction Pgas/(Pgas+Prad)
198061
def Gamma1_gasrad(beta):
    ''' 
    

    '''
    Gamma3minus1 = (old_div(2.,3.))*(4.-3.*beta)/(8.-7.*beta) 
    Gamma1 = beta + (4.-3.*beta) * Gamma3minus1
    return Gamma1
1355867
public static double gammaCdf(double x, double a, double b) {
        if(a<=0 || b<=0) {
            throw new IllegalArgumentException("All the parameters must be positive.");
        }
        
        double GammaCdf = ContinuousDistributions.gammaCdf(x/b, a);
        
        return GammaCdf;
    }
1370914
def lpgamma(p, a):
    
    res = p * (p - 1.0) / 4.0 * log(pi)
    for ii in range(1, p+1):
        res += gammaln(a + (1.0 - ii) / 2.0)

    return res
789882
def meff_e_Gamma(self, **kwargs):
        '''
        
        '''
        Eg = self.Eg_Gamma(**kwargs)
        Delta_SO = self.Delta_SO(**kwargs)
        Ep = self.Ep(**kwargs)
        F = self.F(**kwargs)
        return 1./((1.+2.*F)+(Ep*(Eg+2.*Delta_SO/3.))/(Eg*(Eg+Delta_SO)))
154349
def cdf(x, df, loc=0.0, scale=1.0, gamma = 1.0):
        
        result = np.zeros(x.shape[0])
        result[x<0] = 2.0/(np.power(gamma,2) + 1.0)*ss.t.cdf(gamma*(x[x-loc < 0]-loc[x-loc < 0])/scale, df=df)
        result[x>=0] = 1.0/(np.power(gamma,2) + 1.0) + 2.0/((1.0/np.power(gamma,2)) + 1.0)*(ss.t.cdf((x[x-loc >= 0]-loc[x-loc >= 0])/(gamma*scale), df=df)-0.5)
        return result
1337553
def beta(C, HIGHSCALE, newphys=True):
    """"""

    g = C["g"]
    gp = C["gp"]
    gs = C["gs"]
    m2 = C["m2"]
    Lambda = C["Lambda"]
    Gu = C["Gu"]
    Gd = C["Gd"]
    Ge = C["Ge"]


    Eta1 = (3*np.trace(C["uphi"] @ Gu.conj().T) \
      + 3*np.trace(C["dphi"] @ Gd.conj().T) \
      + np.trace(C["ephi"] @ Ge.conj().T) \
      + 3*np.conj(np.trace(C["uphi"] @ Gu.conj().T)) \
      + 3*np.conj(np.trace(C["dphi"] @ Gd.conj().T)) \
      + np.conj(np.trace(C["ephi"] @ Ge.conj().T)))/2
    Eta2 = -6*np.trace(C["phiq3"] @ Gu @ Gu.conj().T) \
      - 6*np.trace(C["phiq3"] @ Gd @ Gd.conj().T) \
      - 2*np.trace(C["phil3"] @ Ge @ Ge.conj().T) \
      + 3*(np.trace(C["phiud"] @ Gd.conj().T @ Gu) \
      + np.conj(np.trace(C["phiud"] @ Gd.conj().T @ Gu)))
    Eta3 = 3*np.trace(C["phiq1"] @ Gd @ Gd.conj().T) \
      - 3*np.trace(C["phiq1"] @ Gu @ Gu.conj().T) \
      + 9*np.trace(C["phiq3"] @ Gd @ Gd.conj().T) \
      + 9*np.trace(C["phiq3"] @ Gu @ Gu.conj().T) \
      + 3*np.trace(C["phiu"] @ Gu.conj().T @ Gu) \
      - 3*np.trace(C["phid"] @ Gd.conj().T @ Gd) \
      - 3*(np.trace(C["phiud"] @ Gd.conj().T @ Gu) \
      + np.conj(np.trace(C["phiud"] @ Gd.conj().T @ Gu))) \
      + np.trace(C["phil1"] @ Ge @ Ge.conj().T) \
      + 3*np.trace(C["phil3"] @ Ge @ Ge.conj().T) \
      - np.trace(C["phie"] @ Ge.conj().T @ Ge)
    Eta4 = 12*np.trace(C["phiq1"] @ Gd @ Gd.conj().T) \
      - 12*np.trace(C["phiq1"] @ Gu @ Gu.conj().T) \
      + 12*np.trace(C["phiu"] @ Gu.conj().T @ Gu) \
      - 12*np.trace(C["phid"] @ Gd.conj().T @ Gd) \
      + 6*(np.trace(C["phiud"] @ Gd.conj().T @ Gu) \
      + np.conj(np.trace(C["phiud"] @ Gd.conj().T @ Gu))) \
      + 4*np.trace(C["phil1"] @ Ge @ Ge.conj().T) \
      - 4*np.trace(C["phie"] @ Ge.conj().T @ Ge)
    Eta5 = 1j*3/2*(np.trace(Gd @ C["dphi"].conj().T) \
      - np.conj(np.trace(Gd @ C["dphi"].conj().T))) \
      - 1j*3/2*(np.trace(Gu @ C["uphi"].conj().T) \
      - np.conj(np.trace(Gu @ C["uphi"].conj().T))) \
      + 1j*1/2*(np.trace(Ge @ C["ephi"].conj().T) \
      - np.conj(np.trace(Ge @ C["ephi"].conj().T)))

    GammaH = np.trace(3*Gu @ Gu.conj().T + 3*Gd @ Gd.conj().T + Ge @ Ge.conj().T)
    Gammaq = 1/2*(Gu @ Gu.conj().T + Gd @ Gd.conj().T)
    Gammau = Gu.conj().T @ Gu
    Gammad = Gd.conj().T @ Gd
    Gammal = 1/2*Ge @ Ge.conj().T
    Gammae = Ge.conj().T @ Ge

    Beta = OrderedDict()

    Beta["g"] = -19/6*g**3 - 8*g*m2/HIGHSCALE**2*C["phiW"]

    Beta["gp"] = 41/6*gp**3  - 8*gp*m2/HIGHSCALE**2*C["phiB"]

    Beta["gs"] = -7*gs**3  - 8*gs*m2/HIGHSCALE**2*C["phiG"]

    Beta["Lambda"] = 12*Lambda**2 \
      + 3/4*gp**4 + 3/2*g**2*gp**2 + 9/4*g**4 - 3*(gp**2 + 3*g**2)*Lambda \
      + 4*Lambda*GammaH \
      - 4*(3*np.trace(Gd @ Gd.conj().T @ Gd @ Gd.conj().T) \
      + 3*np.trace(Gu @ Gu.conj().T @ Gu @ Gu.conj().T) \
      + np.trace(Ge @ Ge.conj().T @ Ge @ Ge.conj().T)) \
      + 4*m2/HIGHSCALE**2*(12*C["phi"] \
      + (-16*Lambda + 10/3*g**2)*C["phiBox"] \
      + (6*Lambda + 3/2*(gp**2 - g**2))*C["phiD"] \
      + 2*(Eta1 + Eta2) \
      + 9*g**2*C["phiW"] \
      + 3*gp**2*C["phiB"] \
      + 3*g*gp*C["phiWB"] \
      + 4/3*g**2*(np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"])))

    Beta["m2"] = m2*(6*Lambda - 9/2*g**2 - 3/2*gp**2 \
      + 2*GammaH + 4*m2/HIGHSCALE**2*(C["phiD"] \
      - 2*C["phiBox"]))

    Beta["Gu"] = 3/2*(Gu @ Gu.conj().T @ Gu - Gd @ Gd.conj().T @ Gu) \
      + (GammaH - 9/4*g**2 - 17/12*gp**2 - 8*gs**2)*Gu \
      + 2*m2/HIGHSCALE**2*(3*C["uphi"] \
      + 1/2*(C["phiD"] - 2*C["phiBox"])*Gu \
      - C["phiq1"].conj().T @ Gu \
      + 3*C["phiq3"].conj().T @ Gu \
      + Gu @ C["phiu"].conj().T \
      - Gd @ C["phiud"].conj().T \
      - 2*(my_einsum("rpts,pt", C["qu1"], Gu) \
      + 4/3*my_einsum("rpts,pt", C["qu8"], Gu)) \
      - my_einsum("ptrs,pt", C["lequ1"], np.conj(Ge)) \
      + 3*my_einsum("rspt,pt", C["quqd1"], np.conj(Gd)) \
      + 1/2*(my_einsum("psrt,pt", C["quqd1"], np.conj(Gd)) \
      + 4/3*my_einsum("psrt,pt", C["quqd8"], np.conj(Gd))))

    Beta["Gd"] = 3/2*(Gd @ Gd.conj().T @ Gd - Gu @ Gu.conj().T @ Gd) \
      + (GammaH - 9/4*g**2 - 5/12*gp**2 - 8*gs**2)*Gd \
      + 2*m2/HIGHSCALE**2*(3*C["dphi"] + 1/2*(C["phiD"] \
      - 2*C["phiBox"])*Gd \
      + C["phiq1"].conj().T @ Gd \
      + 3*C["phiq3"].conj().T @ Gd \
      - Gd @ C["phid"].conj().T \
      - Gu @ C["phiud"] \
      - 2*(my_einsum("rpts,pt", C["qd1"], Gd) \
      + 4/3*my_einsum("rpts,pt", C["qd8"], Gd)) \
      + my_einsum("ptsr,pt", np.conj(C["ledq"]), Ge) \
      + 3*my_einsum("ptrs,pt", C["quqd1"], np.conj(Gu)) \
      + 1/2*(my_einsum("rpts,tp", C["quqd1"], np.conj(Gu)) \
      + 4/3*my_einsum("rpts,tp", C["quqd8"], np.conj(Gu))))

    Beta["Ge"] = 3/2*Ge @ Ge.conj().T @ Ge + (GammaH \
      - 3/4*(3*g**2 + 5*gp**2))*Ge + 2*m2/HIGHSCALE**2*(3*C["ephi"] \
      + 1/2*(C["phiD"] - 2*C["phiBox"])*Ge \
      + C["phil1"].conj().T @ Ge \
      + 3*C["phil3"].conj().T @ Ge \
      - Ge @ C["phie"].conj().T \
      - 2*my_einsum("rpts,pt", C["le"], Ge) \
      + 3*my_einsum("rspt,tp", C["ledq"], Gd) \
      - 3*my_einsum("rspt,pt", C["lequ1"], np.conj(Gu)))

    Beta["Theta"] = -128*np.pi**2/g**2*m2/HIGHSCALE**2*C["phiWtilde"]

    Beta["Thetap"] = -128*np.pi**2/gp**2*m2/HIGHSCALE**2*C["phiBtilde"]

    Beta["Thetas"] = -128*np.pi**2/gs**2*m2/HIGHSCALE**2*C["phiGtilde"]

    if not newphys:
        # if there is no new physics, generate a dictionary with zero
        # Wilson coefficients (i.e. zero beta functions)
        BetaSM = C_array2dict(np.zeros(5000))
        BetaSM.update(Beta)
        return BetaSM


    XiB = 2/3*(C["phiBox"] + C["phiD"]) \
      + 8/3*( - np.trace(C["phil1"]) + np.trace(C["phiq1"]) \
      - np.trace(C["phie"]) \
      + 2*np.trace(C["phiu"]) - np.trace(C["phid"]))
    Xie = 2*my_einsum("prst,rs", C["le"], Ge) \
      - 3*my_einsum("ptsr,rs", C["ledq"], Gd) \
      + 3*my_einsum("ptsr,sr", C["lequ1"], np.conj(Gu))
    Xid = 2*(my_einsum("prst,rs", C["qd1"], Gd) \
      + 4/3*my_einsum("prst,rs", C["qd8"], Gd)) \
      - (3*my_einsum("srpt,sr", C["quqd1"], np.conj(Gu)) \
      + 1/2*(my_einsum("prst,sr", C["quqd1"], np.conj(Gu)) \
      + 4/3*my_einsum("prst,sr", C["quqd8"], np.conj(Gu)))) \
      - my_einsum("srtp,sr", np.conj(C["ledq"]), Ge)
    Xiu = 2*(my_einsum("prst,rs", C["qu1"], Gu) \
      + 4/3*my_einsum("prst,rs", C["qu8"], Gu)) \
      - (3*my_einsum("ptsr,sr", C["quqd1"], np.conj(Gd)) \
      + 1/2*(my_einsum("stpr,sr", C["quqd1"], np.conj(Gd)) \
      + 4/3*my_einsum("stpr,sr", C["quqd8"], np.conj(Gd)))) \
      + my_einsum("srpt,sr", C["lequ1"], np.conj(Ge))

    Beta["G"] = 15*gs**2*C["G"]

    Beta["Gtilde"] = 15*gs**2*C["Gtilde"]

    Beta["W"] = 29/2*g**2*C["W"]

    Beta["Wtilde"] = 29/2*g**2*C["Wtilde"]

    #c.c.
    Beta["phi"] = -9/2*(3*g**2 \
      + gp**2)*C["phi"] \
      + Lambda*(20/3*g**2*C["phiBox"] \
      + 3*(gp**2 \
      - g**2)*C["phiD"]) \
      - 3/4*(g**2 \
      + gp**2)**2*C["phiD"] \
      + 6*Lambda*(3*g**2*C["phiW"] \
      + gp**2*C["phiB"] \
      + g*gp*C["phiWB"]) \
      - 3*(g**2*gp**2 \
      + 3*g**4)*C["phiW"] \
      - 3*(gp**4 \
      + g**2*gp**2)*C["phiB"] \
      - 3*(g*gp**3 \
      + g**3*gp)*C["phiWB"] \
      + 8/3*Lambda*g**2*(np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"])) \
      + 54*Lambda*C["phi"] \
      - 40*Lambda**2*C["phiBox"] \
      + 12*Lambda**2*C["phiD"] \
      + 4*Lambda*(Eta1 \
      + Eta2) \
      - 4*(3*np.trace(C["uphi"] @ Gu.conj().T @ Gu @ Gu.conj().T) \
      + 3*np.trace(C["dphi"] @ Gd.conj().T @ Gd @ Gd.conj().T) \
      + np.trace(C["ephi"] @ Ge.conj().T @ Ge @ Ge.conj().T) \
      + 3*np.conj(np.trace(C["uphi"] @ Gu.conj().T @ Gu @ Gu.conj().T)) \
      + 3*np.conj(np.trace(C["dphi"] @ Gd.conj().T @ Gd @ Gd.conj().T)) \
      + np.conj(np.trace(C["ephi"] @ Ge.conj().T @ Ge @ Ge.conj().T))) \
      + 6*GammaH*C["phi"]

    Beta["phiBox"] = -(4*g**2 \
      + 4/3*gp**2)*C["phiBox"] \
      + 5/3*gp**2*C["phiD"] \
      + 2*g**2*(np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"])) \
      + 2/3*gp**2*(2*np.trace(C["phiu"]) \
      - np.trace(C["phid"]) \
      - np.trace(C["phie"]) \
      + np.trace(C["phiq1"]) \
      - np.trace(C["phil1"])) \
      + 12*Lambda*C["phiBox"] \
      - 2*Eta3 \
      + 4*GammaH*C["phiBox"]

    Beta["phiD"] = 20/3*gp**2*C["phiBox"] \
      + (9/2*g**2 \
      - 5/6*gp**2)*C["phiD"] \
      + 8/3*gp**2*(2*np.trace(C["phiu"]) \
      - np.trace(C["phid"]) \
      - np.trace(C["phie"]) \
      + np.trace(C["phiq1"]) \
      - np.trace(C["phil1"])) \
      + 6*Lambda*C["phiD"] \
      - 2*Eta4 \
      + 4*GammaH*C["phiD"]

    #c.c.
    Beta["phiG"] = (-3/2*gp**2 \
      - 9/2*g**2 \
      - 14*gs**2)*C["phiG"] \
      + 6*Lambda*C["phiG"] \
      - 2*gs*(np.trace(C["uG"] @ Gu.conj().T) \
      + np.trace(C["dG"] @ Gd.conj().T) \
      + np.conj(np.trace(C["uG"] @ Gu.conj().T)) \
      + np.conj(np.trace(C["dG"] @ Gd.conj().T))) \
      + 2*GammaH*C["phiG"]

    #c.c.
    Beta["phiB"] = (85/6*gp**2 \
      - 9/2*g**2)*C["phiB"] \
      + 3*g*gp*C["phiWB"] \
      + 6*Lambda*C["phiB"] \
      + gp*( \
      - 5*np.trace(C["uB"] @ Gu.conj().T) \
      + np.trace(C["dB"] @ Gd.conj().T) \
      + 3*np.trace(C["eB"] @ Ge.conj().T) \
      - 5*np.conj(np.trace(C["uB"] @ Gu.conj().T)) \
      + np.conj(np.trace(C["dB"] @ Gd.conj().T)) \
      + 3*np.conj(np.trace(C["eB"] @ Ge.conj().T))) \
      + 2*GammaH*C["phiB"]

    #c.c.
    Beta["phiW"] = (-3/2*gp**2 \
      - 53/6*g**2)*C["phiW"] \
      + g*gp*C["phiWB"] \
      - 15*g**3*C["W"] \
      + 6*Lambda*C["phiW"] \
      - g*(3*np.trace(C["uW"] @ Gu.conj().T) \
      + 3*np.trace(C["dW"] @ Gd.conj().T) \
      + np.trace(C["eW"] @ Ge.conj().T) \
      + 3*np.conj(np.trace(C["uW"] @ Gu.conj().T)) \
      + 3*np.conj(np.trace(C["dW"] @ Gd.conj().T)) \
      + np.conj(np.trace(C["eW"] @ Ge.conj().T))) \
      + 2*GammaH*C["phiW"]

    #c.c.
    Beta["phiWB"] = (19/3*gp**2 \
      + 4/3*g**2)*C["phiWB"] \
      + 2*g*gp*(C["phiB"] \
      + C["phiW"]) \
      + 3*g**2*gp*C["W"] \
      + 2*Lambda*C["phiWB"] \
      + g*(3*np.trace(C["uB"] @ Gu.conj().T) \
      - 3*np.trace(C["dB"] @ Gd.conj().T) \
      - np.trace(C["eB"] @ Ge.conj().T) \
      + 3*np.conj(np.trace(C["uB"] @ Gu.conj().T)) \
      - 3*np.conj(np.trace(C["dB"] @ Gd.conj().T)) \
      - np.conj(np.trace(C["eB"] @ Ge.conj().T))) \
      + gp*(5*np.trace(C["uW"] @ Gu.conj().T) \
      + np.trace(C["dW"] @ Gd.conj().T) \
      + 3*np.trace(C["eW"] @ Ge.conj().T) \
      + 5*np.conj(np.trace(C["uW"] @ Gu.conj().T)) \
      + np.conj(np.trace(C["dW"] @ Gd.conj().T)) \
      + 3*np.conj(np.trace(C["eW"] @ Ge.conj().T))) \
      + 2*GammaH*C["phiWB"]

    #problem with i as I*iCPV
    Beta["phiGtilde"] = (-3/2*gp**2 \
      - 9/2*g**2 \
      - 14*gs**2)*C["phiGtilde"] \
      + 6*Lambda*C["phiGtilde"] \
      + 2j*gs*(np.trace(C["uG"] @ Gu.conj().T) \
      + np.trace(C["dG"] @ Gd.conj().T) \
      - np.conj(np.trace(C["uG"] @ Gu.conj().T)) \
      - np.conj(np.trace(C["dG"] @ Gd.conj().T))) \
      + 2*GammaH*C["phiGtilde"]

    #i
    Beta["phiBtilde"] = (85/6*gp**2 \
      - 9/2*g**2)*C["phiBtilde"] \
      + 3*g*gp*C["phiWtildeB"] \
      + 6*Lambda*C["phiBtilde"] \
      - 1j*gp*( \
      - 5*np.trace(C["uB"] @ Gu.conj().T) \
      + np.trace(C["dB"] @ Gd.conj().T) \
      + 3*np.trace(C["eB"] @ Ge.conj().T) \
      + 5*np.conj(np.trace(C["uB"] @ Gu.conj().T)) \
      - np.conj(np.trace(C["dB"] @ Gd.conj().T)) \
      - 3*np.conj(np.trace(C["eB"] @ Ge.conj().T))) \
      + 2*GammaH*C["phiBtilde"]

    #i
    Beta["phiWtilde"] = (-3/2*gp**2 \
      - 53/6*g**2)*C["phiWtilde"] \
      + g*gp*C["phiWtildeB"] \
      - 15*g**3*C["Wtilde"] \
      + 6*Lambda*C["phiWtilde"] \
      + 1j*g*(3*np.trace(C["uW"] @ Gu.conj().T) \
      + 3*np.trace(C["dW"] @ Gd.conj().T) \
      + np.trace(C["eW"] @ Ge.conj().T) \
      - 3*np.conj(np.trace(C["uW"] @ Gu.conj().T)) \
      - 3*np.conj(np.trace(C["dW"] @ Gd.conj().T)) \
      - np.conj(np.trace(C["eW"] @ Ge.conj().T))) \
      + 2*GammaH*C["phiWtilde"]

    #i
    Beta["phiWtildeB"] = (19/3*gp**2 \
      + 4/3*g**2)*C["phiWtildeB"] \
      + 2*g*gp*(C["phiBtilde"] \
      + C["phiWtilde"]) \
      + 3*g**2*gp*C["Wtilde"] \
      + 2*Lambda*C["phiWtildeB"] \
      - 1j*g*(3*np.trace(C["uB"] @ Gu.conj().T) \
      - 3*np.trace(C["dB"] @ Gd.conj().T) \
      - np.trace(C["eB"] @ Ge.conj().T) \
      - 3*np.conj(np.trace(C["uB"] @ Gu.conj().T)) \
      + 3*np.conj(np.trace(C["dB"] @ Gd.conj().T)) \
      + np.conj(np.trace(C["eB"] @ Ge.conj().T))) \
      - 1j*gp*(5*np.trace(C["uW"] @ Gu.conj().T) \
      + np.trace(C["dW"] @ Gd.conj().T) \
      + 3*np.trace(C["eW"] @ Ge.conj().T) \
      - 5*np.conj(np.trace(C["uW"] @ Gu.conj().T)) \
      - np.conj(np.trace(C["dW"] @ Gd.conj().T)) \
      - 3*np.conj(np.trace(C["eW"] @ Ge.conj().T))) \
      + 2*GammaH*C["phiWtildeB"]

    """(3,3)"""
    #i  #the coefficients of Eta5 is not equal
    Beta["uphi"] = (10/3*g**2*C["phiBox"] \
      + 3/2*(gp**2 \
      - g**2)*C["phiD"] \
      + 32*gs**2*(C["phiG"] \
      + 1j*C["phiGtilde"]) \
      + 9*g**2*(C["phiW"] \
      + 1j*C["phiWtilde"]) \
      + 17/3*gp**2*(C["phiB"] \
      + 1j*C["phiBtilde"]) \
      - g*gp*(C["phiWB"] \
      + 1j*C["phiWtildeB"]) \
      + 4/3*g**2*(np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"])))*Gu \
      - (35/12*gp**2 \
      + 27/4*g**2 \
      + 8*gs**2)*C["uphi"] \
      - gp*(5*gp**2 \
      - 3*g**2)*C["uB"] \
      + g*(5*gp**2 \
      - 9*g**2)*C["uW"] \
      - (3*g**2 \
      - gp**2)*Gu @ C["phiu"] \
      + 3*g**2*Gd @ C["phiud"].conj().T \
      + 4*gp**2*C["phiq1"] @ Gu \
      - 4*gp**2*C["phiq3"] @ Gu \
      - 5*gp*(C["uB"] @ Gu.conj().T @ Gu \
      + Gu @ Gu.conj().T @ C["uB"]) \
      - 3*g*(C["uW"] @ Gu.conj().T @ Gu \
      - Gu @ Gu.conj().T @ C["uW"]) \
      - 16*gs*(C["uG"] @ Gu.conj().T @ Gu \
      + Gu @ Gu.conj().T @ C["uG"]) \
      - 12*g*Gd @ Gd.conj().T @ C["uW"] \
      - 6*g*C["dW"] @ Gd.conj().T @ Gu \
      + Lambda*(12*C["uphi"] \
      - 2*C["phiq1"] @ Gu \
      + 6*C["phiq3"] @ Gu \
      + 2*Gu @ C["phiu"] \
      - 2*Gd @ C["phiud"].conj().T \
      - 2*C["phiBox"]*Gu \
      + C["phiD"]*Gu \
      - 4*my_einsum("rpts,pt", C["qu1"], Gu) \
      - 16/3*my_einsum("rpts,pt", C["qu8"], Gu) \
      - 2*my_einsum("ptrs,pt", C["lequ1"], np.conj(Ge)) \
      + 6*my_einsum("rspt,pt", C["quqd1"], np.conj(Gd)) \
      + my_einsum("psrt,pt", C["quqd1"], np.conj(Gd)) \
      + 4/3*my_einsum("psrt,pt", C["quqd8"], np.conj(Gd))) \
      + 2*(Eta1 \
      + Eta2 \
      - 1j*Eta5)*Gu \
      + (C["phiD"] \
      - 6*C["phiBox"])*Gu @ Gu.conj().T @ Gu \
      - 2*C["phiq1"] @ Gu @ Gu.conj().T @ Gu \
      + 6*C["phiq3"] @ Gd @ Gd.conj().T @ Gu \
      + 2*Gu @ Gu.conj().T @ Gu @ C["phiu"] \
      - 2*Gd @ Gd.conj().T @ Gd @ C["phiud"].conj().T \
      + 8*(my_einsum("rpts,pt", C["qu1"], Gu @ Gu.conj().T @ Gu) \
      + 4/3*my_einsum("rpts,pt", C["qu8"], Gu @ Gu.conj().T @ Gu)) \
      - 2*(my_einsum("tsrp,pt", C["quqd1"], Gd.conj().T @ Gd @ Gd.conj().T) \
      + 4/3*my_einsum("tsrp,pt", C["quqd8"], Gd.conj().T @ Gd @ Gd.conj().T)) \
      - 12*my_einsum("rstp,pt", C["quqd1"], Gd.conj().T @ Gd @ Gd.conj().T) \
      + 4*my_einsum("tprs,pt", C["lequ1"], Ge.conj().T @ Ge @ Ge.conj().T) \
      + 4*C["uphi"] @ Gu.conj().T @ Gu \
      + 5*Gu @ Gu.conj().T @ C["uphi"] \
      - 2*Gd @ C["dphi"].conj().T @ Gu \
      - C["dphi"] @ Gd.conj().T @ Gu \
      - 2*Gd @ Gd.conj().T @ C["uphi"] \
      + 3*GammaH*C["uphi"] \
      + Gammaq @ C["uphi"] \
      + C["uphi"] @ Gammau

    #i  #Eta5
    Beta["dphi"] = (10/3*g**2*C["phiBox"] \
      + 3/2*(gp**2 \
      - g**2)*C["phiD"] \
      + 32*gs**2*(C["phiG"] \
      + 1j*C["phiGtilde"]) \
      + 9*g**2*(C["phiW"] \
      + 1j*C["phiWtilde"]) \
      + 5/3*gp**2*(C["phiB"] \
      + 1j*C["phiBtilde"]) \
      + g*gp*(C["phiWB"] \
      + 1j*C["phiWtildeB"]) \
      + 4/3*g**2*(np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"])))*Gd \
      - (23/12*gp**2 \
      + 27/4*g**2 \
      + 8*gs**2)*C["dphi"] \
      - gp*(3*g**2 \
      - gp**2)*C["dB"] \
      - g*(9*g**2 \
      - gp**2)*C["dW"] \
      + (3*g**2 \
      + gp**2)*Gd @ C["phid"] \
      + 3*g**2*Gu @ C["phiud"] \
      - 2*gp**2*C["phiq1"] @ Gd \
      - 2*gp**2*C["phiq3"] @ Gd \
      + gp*(C["dB"] @ Gd.conj().T @ Gd \
      + Gd @ Gd.conj().T @ C["dB"]) \
      - 3*g*(C["dW"] @ Gd.conj().T @ Gd \
      - Gd @ Gd.conj().T @ C["dW"]) \
      - 16*gs*(C["dG"] @ Gd.conj().T @ Gd \
      + Gd @ Gd.conj().T @ C["dG"]) \
      - 12*g*Gu @ Gu.conj().T @ C["dW"] \
      - 6*g*C["uW"] @ Gu.conj().T @ Gd \
      + Lambda*(12*C["dphi"] \
      + 2*C["phiq1"] @ Gd \
      + 6*C["phiq3"] @ Gd \
      - 2*Gd @ C["phid"] \
      - 2*Gu @ C["phiud"] \
      - 2*C["phiBox"]*Gd \
      + C["phiD"]*Gd \
      - 4*my_einsum("rpts,pt", C["qd1"], Gd) \
      - 16/3*my_einsum("rpts,pt", C["qd8"], Gd) \
      + 2*my_einsum("ptsr,pt", np.conj(C["ledq"]), Ge) \
      + 6*my_einsum("ptrs,pt", C["quqd1"], np.conj(Gu)) \
      + my_einsum("rtps,pt", C["quqd1"], np.conj(Gu)) \
      + 4/3*my_einsum("rtps,pt", C["quqd8"], np.conj(Gu))) \
      + 2*(Eta1 \
      + Eta2 \
      + 1j*Eta5)*Gd \
      + (C["phiD"] \
      - 6*C["phiBox"])*Gd @ Gd.conj().T @ Gd \
      + 2*C["phiq1"] @ Gd @ Gd.conj().T @ Gd \
      + 6*C["phiq3"] @ Gu @ Gu.conj().T @ Gd \
      - 2*Gd @ Gd.conj().T @ Gd @ C["phid"] \
      - 2*Gu @ Gu.conj().T @ Gu @ C["phiud"] \
      + 8*(my_einsum("rpts,pt", C["qd1"], Gd @ Gd.conj().T @ Gd) \
      + 4/3*my_einsum("rpts,pt", C["qd8"], Gd @ Gd.conj().T @ Gd)) \
      - 2*(my_einsum("rpts,pt", C["quqd1"], Gu.conj().T @ Gu @ Gu.conj().T) \
      + 4/3*my_einsum("rpts,pt", C["quqd8"], Gu.conj().T @ Gu @ Gu.conj().T)) \
      - 12*my_einsum("tprs,pt", C["quqd1"], Gu @ Gu.conj().T @ Gu) \
      - 4*my_einsum("ptsr,pt", np.conj(C["ledq"]), Ge @ Ge.conj().T @ Ge) \
      + 4*C["dphi"] @ Gd.conj().T @ Gd \
      + 5*Gd @ Gd.conj().T @ C["dphi"] \
      - 2*Gu @ C["uphi"].conj().T @ Gd \
      - C["uphi"] @ Gu.conj().T @ Gd \
      - 2*Gu @ Gu.conj().T @ C["dphi"] \
      + 3*GammaH*C["dphi"] \
      + Gammaq @ C["dphi"] \
      + C["dphi"] @ Gammad

    #i
    Beta["ephi"] = (10/3*g**2*C["phiBox"] \
      + 3/2*(gp**2 \
      - g**2)*C["phiD"] \
      + 9*g**2*(C["phiW"] \
      + 1j*C["phiWtilde"]) \
      + 15*gp**2*(C["phiB"] \
      + 1j*C["phiBtilde"]) \
      - 3*g*gp*(C["phiWB"] \
      + 1j*C["phiWtildeB"]) \
      + 4/3*g**2*(np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"])))*Ge \
      - 3/4*(7*gp**2 \
      + 9*g**2)*C["ephi"] \
      - 3*gp*(g**2 \
      - 3*gp**2)*C["eB"] \
      - 9*g*(g**2 \
      - gp**2)*C["eW"] \
      + 3*(g**2 \
      - gp**2)*Ge @ C["phie"] \
      - 6*gp**2*C["phil1"] @ Ge \
      - 6*gp**2*C["phil3"] @ Ge \
      + 9*gp*(C["eB"] @ Ge.conj().T @ Ge \
      + Ge @ Ge.conj().T @ C["eB"]) \
      - 3*g*(C["eW"] @ Ge.conj().T @ Ge \
      - Ge @ Ge.conj().T @ C["eW"]) \
      + Lambda*(12*C["ephi"] \
      + 2*C["phil1"] @ Ge \
      + 6*C["phil3"] @ Ge \
      - 2*Ge @ C["phie"] \
      - 2*C["phiBox"]*Ge \
      + C["phiD"]*Ge \
      - 4*my_einsum("rpts,pt", C["le"], Ge) \
      + 6*my_einsum("rspt,tp", C["ledq"], Gd) \
      - 6*my_einsum("rspt,pt", C["lequ1"], np.conj(Gu))) \
      + 2*(Eta1 \
      + Eta2 \
      + 1j*Eta5)*Ge \
      + (C["phiD"] \
      - 6*C["phiBox"])*Ge @ Ge.conj().T @ Ge \
      + 2*C["phil1"] @ Ge @ Ge.conj().T @ Ge \
      - 2*Ge @ Ge.conj().T @ Ge @ C["phie"] \
      + 8*my_einsum("rpts,pt", C["le"], Ge @ Ge.conj().T @ Ge) \
      - 12*my_einsum("rspt,tp", C["ledq"], Gd @ Gd.conj().T @ Gd) \
      + 12*my_einsum("rstp,pt", C["lequ1"], Gu.conj().T @ Gu @ Gu.conj().T) \
      + 4*C["ephi"] @ Ge.conj().T @ Ge \
      + 5*Ge @ Ge.conj().T @ C["ephi"] \
      + 3*GammaH*C["ephi"] \
      + Gammal @ C["ephi"] \
      + C["ephi"] @ Gammae

    #i
    Beta["eW"] = 1/12*(3*gp**2 \
      - 11*g**2)*C["eW"] \
      - 1/2*g*gp*C["eB"] \
      - (g*(C["phiW"] \
      + 1j*C["phiWtilde"]) \
      - 3/2*gp*(C["phiWB"] \
      + 1j*C["phiWtildeB"]))*Ge \
      - 6*g*my_einsum("rspt,pt", C["lequ3"], np.conj(Gu)) \
      + C["eW"] @ Ge.conj().T @ Ge \
      + GammaH*C["eW"] \
      + Gammal @ C["eW"] \
      + C["eW"] @ Gammae

    #i
    Beta["eB"] = 1/4*(151/3*gp**2 \
      - 9*g**2)*C["eB"] \
      - 3/2*g*gp*C["eW"] \
      - (3/2*g*(C["phiWB"] \
      + 1j*C["phiWtildeB"]) \
      - 3*gp*(C["phiB"] \
      + 1j*C["phiBtilde"]))*Ge \
      + 10*gp*my_einsum("rspt,pt", C["lequ3"], np.conj(Gu)) \
      + C["eB"] @ Ge.conj().T @ Ge \
      + 2*Ge @ Ge.conj().T @ C["eB"] \
      + GammaH*C["eB"] \
      + Gammal @ C["eB"] \
      + C["eB"] @ Gammae

    #i
    Beta["uG"] = -1/36*(81*g**2 \
      + 19*gp**2 \
      + 204*gs**2)*C["uG"] \
      + 6*g*gs*C["uW"] \
      + 10/3*gp*gs*C["uB"] \
      - gs*(4*(C["phiG"] \
      + 1j*C["phiGtilde"]) \
      - 9*gs*(C["G"] \
      + 1j*C["Gtilde"]))*Gu \
      - gs*(my_einsum("psrt,pt", C["quqd1"], np.conj(Gd)) \
      - 1/6*my_einsum("psrt,pt", C["quqd8"], np.conj(Gd))) \
      + 2*Gu @ Gu.conj().T @ C["uG"] \
      - 2*Gd @ Gd.conj().T @ C["uG"] \
      - C["dG"] @ Gd.conj().T @ Gu \
      + C["uG"] @ Gu.conj().T @ Gu \
      + GammaH*C["uG"] \
      + Gammaq @ C["uG"] \
      + C["uG"] @ Gammau

    #i
    Beta["uW"] = -1/36*(33*g**2 \
      + 19*gp**2 \
      - 96*gs**2)*C["uW"] \
      + 8/3*g*gs*C["uG"] \
      - 1/6*g*gp*C["uB"] \
      - (g*(C["phiW"] \
      + 1j*C["phiWtilde"]) \
      - 5/6*gp*(C["phiWB"] \
      + 1j*C["phiWtildeB"]))*Gu \
      + g/4*(my_einsum("psrt,pt", C["quqd1"], np.conj(Gd)) \
      + 4/3*my_einsum("psrt,pt", C["quqd8"], np.conj(Gd))) \
      - 2*g*my_einsum("ptrs,pt", C["lequ3"], np.conj(Ge)) \
      + 2*Gd @ Gd.conj().T @ C["uW"] \
      - C["dW"] @ Gd.conj().T @ Gu \
      + C["uW"] @ Gu.conj().T @ Gu \
      + GammaH*C["uW"] \
      + Gammaq @ C["uW"] \
      + C["uW"] @ Gammau

    #i
    Beta["uB"] = -1/36*(81*g**2 \
      - 313*gp**2 \
      - 96*gs**2)*C["uB"] \
      + 40/9*gp*gs*C["uG"] \
      - 1/2*g*gp*C["uW"] \
      - (-3/2*g*(C["phiWB"] \
      + 1j*C["phiWtildeB"]) \
      + 5/3*gp*(C["phiB"] \
      + 1j*C["phiBtilde"]))*Gu \
      + gp/12*(my_einsum("psrt,pt", C["quqd1"], np.conj(Gd)) \
      + 4/3*my_einsum("psrt,pt", C["quqd8"], np.conj(Gd))) \
      - 6*gp*my_einsum("ptrs,pt", C["lequ3"], np.conj(Ge)) \
      + 2*Gu @ Gu.conj().T @ C["uB"] \
      - 2*Gd @ Gd.conj().T @ C["uB"] \
      - C["dB"] @ Gd.conj().T @ Gu \
      + C["uB"] @ Gu.conj().T @ Gu \
      + GammaH*C["uB"] \
      + Gammaq @ C["uB"] \
      + C["uB"] @ Gammau

    #i
    Beta["dG"] = -1/36*(81*g**2 \
      + 31*gp**2 \
      + 204*gs**2)*C["dG"] \
      + 6*g*gs*C["dW"] \
      - 2/3*gp*gs*C["dB"] \
      - gs*(4*(C["phiG"] \
      + 1j*C["phiGtilde"]) \
      - 9*gs*(C["G"] \
      + 1j*C["Gtilde"]))*Gd \
      - gs*(my_einsum("rtps,pt", C["quqd1"], np.conj(Gu)) \
      - 1/6*my_einsum("rtps,pt", C["quqd8"], np.conj(Gu))) \
      - 2*Gu @ Gu.conj().T @ C["dG"] \
      + 2*Gd @ Gd.conj().T @ C["dG"] \
      - C["uG"] @ Gu.conj().T @ Gd \
      + C["dG"] @ Gd.conj().T @ Gd \
      + GammaH*C["dG"] \
      + Gammaq @ C["dG"] \
      + C["dG"] @ Gammad

    #i
    Beta["dW"] = -1/36*(33*g**2 \
      + 31*gp**2 \
      - 96*gs**2)*C["dW"] \
      + 8/3*g*gs*C["dG"] \
      + 5/6*g*gp*C["dB"] \
      - (g*(C["phiW"] \
      + 1j*C["phiWtilde"]) \
      - gp/6*(C["phiWB"] \
      + 1j*C["phiWtildeB"]))*Gd \
      + g/4*(my_einsum("rtps,pt", C["quqd1"], np.conj(Gu)) \
      + 4/3*my_einsum("rtps,pt", C["quqd8"], np.conj(Gu))) \
      + 2*Gu @ Gu.conj().T @ C["dW"] \
      - C["uW"] @ Gu.conj().T @ Gd \
      + C["dW"] @ Gd.conj().T @ Gd \
      + GammaH*C["dW"] \
      + Gammaq @ C["dW"] \
      + C["dW"] @ Gammad

    #i
    Beta["dB"] = -1/36*(81*g**2 \
      - 253*gp**2 \
      - 96*gs**2)*C["dB"] \
      - 8/9*gp*gs*C["dG"] \
      + 5/2*g*gp*C["dW"] \
      - (3/2*g*(C["phiWB"] \
      + 1j*C["phiWtildeB"]) \
      - gp/3*(C["phiB"] \
      + 1j*C["phiBtilde"]))*Gd \
      - 5/12*gp*(my_einsum("rtps,pt", C["quqd1"], np.conj(Gu)) \
      + 4/3*my_einsum("rtps,pt", C["quqd8"], np.conj(Gu))) \
      - 2*Gu @ Gu.conj().T @ C["dB"] \
      + 2*Gd @ Gd.conj().T @ C["dB"] \
      - C["uB"] @ Gu.conj().T @ Gd \
      + C["dB"] @ Gd.conj().T @ Gd \
      + GammaH*C["dB"] \
      + Gammaq @ C["dB"] \
      + C["dB"] @ Gammad

    #I3 #coefficient not equal with manual!!!!!!
    Beta["phil1"] = -1/4*XiB*gp**2*I3 \
      + 1/3*gp**2*C["phil1"] \
      - 2/3*gp**2*(my_einsum("rstt", C["ld"]) \
      + my_einsum("rstt", C["le"]) \
      + 2*my_einsum("rstt", C["ll"]) \
      + my_einsum("rtts", C["ll"]) \
      - my_einsum("rstt", C["lq1"]) \
      - 2*my_einsum("rstt", C["lu"])) \
      - 1/2*(C["phiBox"] \
      + C["phiD"])*Ge @ Ge.conj().T \
      - Ge @ C["phie"] @ Ge.conj().T \
      + 3/2*(Ge @ Ge.conj().T @ C["phil1"] \
      + C["phil1"] @ Ge @ Ge.conj().T \
      + 3*Ge @ Ge.conj().T @ C["phil3"] \
      + 3*C["phil3"] @ Ge @ Ge.conj().T) \
      + 2*my_einsum("rspt,tp", C["le"], Ge.conj().T @ Ge) \
      - 2*(2*my_einsum("rspt,tp", C["ll"], Ge @ Ge.conj().T) \
      + my_einsum("rtps,tp", C["ll"], Ge @ Ge.conj().T)) \
      - 6*my_einsum("rspt,tp", C["lq1"], Gd @ Gd.conj().T) \
      + 6*my_einsum("rspt,tp", C["lq1"], Gu @ Gu.conj().T) \
      - 6*my_einsum("rspt,tp", C["lu"], Gu.conj().T @ Gu) \
      + 6*my_einsum("rspt,tp", C["ld"], Gd.conj().T @ Gd) \
      + 2*GammaH*C["phil1"] \
      + Gammal @ C["phil1"] \
      + C["phil1"] @ Gammal

    #I3 #coefficient
    Beta["phil3"] = 2/3*g**2*(1/4*C["phiBox"] \
      + np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"]))*I3 \
      - 17/3*g**2*C["phil3"] \
      + 2/3*g**2*my_einsum("rtts", C["ll"]) \
      + 2*g**2*my_einsum("rstt", C["lq3"]) \
      - 1/2*C["phiBox"]*Ge @ Ge.conj().T \
      + 1/2*(3*Ge @ Ge.conj().T @ C["phil1"] \
      + 3*C["phil1"] @ Ge @ Ge.conj().T \
      + Ge @ Ge.conj().T @ C["phil3"] \
      + C["phil3"] @ Ge @ Ge.conj().T) \
      - 2*(my_einsum("rtps,tp", C["ll"], Ge @ Ge.conj().T)) \
      - 6*my_einsum("rspt,tp", C["lq3"], Gd @ Gd.conj().T) \
      - 6*my_einsum("rspt,tp", C["lq3"], Gu @ Gu.conj().T) \
      + 2*GammaH*C["phil3"] \
      + Gammal @ C["phil3"] \
      + C["phil3"] @ Gammal

    #I3  #coefficient even terms not equal...
    Beta["phie"] = -1/2*XiB*gp**2*I3 \
      + 1/3*gp**2*C["phie"] \
      - 2/3*gp**2*(my_einsum("rstt", C["ed"]) \
      + 4*my_einsum("rstt", C["ee"]) \
      - 2*my_einsum("rstt", C["eu"]) \
      + my_einsum("ttrs", C["le"]) \
      - my_einsum("ttrs", C["qe"])) \
      + (C["phiBox"] \
      + C["phiD"])*Ge.conj().T @ Ge \
      - 2*Ge.conj().T @ C["phil1"] @ Ge \
      + 3*(Ge.conj().T @ Ge @ C["phie"] \
      + C["phie"] @ Ge.conj().T @ Ge) \
      - 2*my_einsum("ptrs,tp", C["le"], Ge @ Ge.conj().T) \
      + 8*my_einsum("rspt,tp", C["ee"], Ge.conj().T @ Ge) \
      - 6*my_einsum("rspt,tp", C["eu"], Gu.conj().T @ Gu) \
      + 6*my_einsum("rspt,tp", C["ed"], Gd.conj().T @ Gd) \
      - 6*my_einsum("ptrs,tp", C["qe"], Gd @ Gd.conj().T) \
      + 6*my_einsum("ptrs,tp", C["qe"], Gu @ Gu.conj().T) \
      + 2*GammaH*C["phie"] \
      + Gammae @ C["phie"] \
      + C["phie"] @ Gammae

    #I3  #coefficient???
    Beta["phiq1"] = 1/12*XiB*gp**2*I3 \
      + 1/3*gp**2*C["phiq1"] \
      - 2/3*gp**2*(my_einsum("ttrs", C["lq1"]) \
      + my_einsum("rstt", C["qd1"]) \
      - 2*my_einsum("rstt", C["qu1"]) \
      + my_einsum("rstt", C["qe"]) \
      - 2*my_einsum("rstt", C["qq1"]) \
      - 1/3*my_einsum("rtts", C["qq1"]) \
      - my_einsum("rtts", C["qq3"])) \
      + 1/2*(C["phiBox"] \
      + C["phiD"])*(Gu @ Gu.conj().T \
      - Gd @ Gd.conj().T) \
      - Gu @ C["phiu"] @ Gu.conj().T \
      - Gd @ C["phid"] @ Gd.conj().T \
      + 2*my_einsum("rspt,tp", C["qe"], Ge.conj().T @ Ge) \
      - 2*my_einsum("ptrs,tp", C["lq1"], Ge @ Ge.conj().T) \
      + 3/2*(Gd @ Gd.conj().T @ C["phiq1"] \
      + Gu @ Gu.conj().T @ C["phiq1"] \
      + C["phiq1"] @ Gd @ Gd.conj().T \
      + C["phiq1"] @ Gu @ Gu.conj().T \
      + 3*Gd @ Gd.conj().T @ C["phiq3"] \
      - 3*Gu @ Gu.conj().T @ C["phiq3"] \
      + 3*C["phiq3"] @ Gd @ Gd.conj().T \
      - 3*C["phiq3"] @ Gu @ Gu.conj().T) \
      - 2*(6*my_einsum("ptrs,tp", C["qq1"], Gd @ Gd.conj().T) \
      + my_einsum("psrt,tp", C["qq1"], Gd @ Gd.conj().T) \
      + 3*my_einsum("psrt,tp", C["qq3"], Gd @ Gd.conj().T) \
      - 6*my_einsum("ptrs,tp", C["qq1"], Gu @ Gu.conj().T) \
      - my_einsum("psrt,tp", C["qq1"], Gu @ Gu.conj().T) \
      - 3*my_einsum("psrt,tp", C["qq3"], Gu @ Gu.conj().T)) \
      - 6*my_einsum("rspt,tp", C["qu1"], Gu.conj().T @ Gu) \
      + 6*my_einsum("rspt,tp", C["qd1"], Gd.conj().T @ Gd) \
      + 2*GammaH*C["phiq1"] \
      + Gammaq @ C["phiq1"] \
      + C["phiq1"] @ Gammaq

    #I3 #co
    Beta["phiq3"] = 2/3*g**2*(1/4*C["phiBox"] \
      + np.trace(C["phil3"]) \
      + 3*np.trace(C["phiq3"]))*I3 \
      - 17/3*g**2*C["phiq3"] \
      + 2/3*g**2*(my_einsum("ttrs", C["lq3"]) \
      + my_einsum("rtts", C["qq1"]) \
      + 6*my_einsum("rstt", C["qq3"]) \
      - my_einsum("rtts", C["qq3"])) \
      - 1/2*C["phiBox"]*(Gu @ Gu.conj().T \
      + Gd @ Gd.conj().T) \
      + 1/2*(3*Gd @ Gd.conj().T @ C["phiq1"] \
      - 3*Gu @ Gu.conj().T @ C["phiq1"] \
      + 3*C["phiq1"] @ Gd @ Gd.conj().T \
      - 3*C["phiq1"] @ Gu @ Gu.conj().T \
      + Gd @ Gd.conj().T @ C["phiq3"] \
      + Gu @ Gu.conj().T @ C["phiq3"] \
      + C["phiq3"] @ Gd @ Gd.conj().T \
      + C["phiq3"] @ Gu @ Gu.conj().T) \
      - 2*(6*my_einsum("rspt,tp", C["qq3"], Gd @ Gd.conj().T) \
      + my_einsum("rtps,tp", C["qq1"], Gd @ Gd.conj().T) \
      - my_einsum("rtps,tp", C["qq3"], Gd @ Gd.conj().T) \
      + 6*my_einsum("rspt,tp", C["qq3"], Gu @ Gu.conj().T) \
      + my_einsum("rtps,tp", C["qq1"], Gu @ Gu.conj().T) \
      - my_einsum("rtps,tp", C["qq3"], Gu @ Gu.conj().T)) \
      - 2*my_einsum("ptrs,tp", C["lq3"], Ge @ Ge.conj().T) \
      + 2*GammaH*C["phiq3"] \
      + Gammaq @ C["phiq3"] \
      + C["phiq3"] @ Gammaq

    #I3 #co
    Beta["phiu"] = 1/3*XiB*gp**2*I3 \
      + 1/3*gp**2*C["phiu"] \
      - 2/3*gp**2*(my_einsum("ttrs", C["eu"]) \
      + my_einsum("ttrs", C["lu"]) \
      - my_einsum("ttrs", C["qu1"]) \
      + my_einsum("rstt", C["ud1"]) \
      - 4*my_einsum("rstt", C["uu"]) \
      - 4/3*my_einsum("rtts", C["uu"])) \
      - (C["phiBox"] \
      + C["phiD"])*Gu.conj().T @ Gu \
      - 2*Gu.conj().T @ C["phiq1"] @ Gu \
      + 3*(Gu.conj().T @ Gu @ C["phiu"] \
      + C["phiu"] @ Gu.conj().T @ Gu) \
      + Gu.conj().T @ Gd @ C["phiud"].conj().T \
      + C["phiud"] @ Gd.conj().T @ Gu \
      - 4*(3*my_einsum("rspt,tp", C["uu"], Gu.conj().T @ Gu) \
      + my_einsum("rtps,tp", C["uu"], Gu.conj().T @ Gu)) \
      + 2*my_einsum("ptrs,tp", C["eu"], Ge.conj().T @ Ge) \
      - 2*my_einsum("ptrs,tp", C["lu"], Ge @ Ge.conj().T) \
      + 6*my_einsum("rspt,tp", C["ud1"], Gd.conj().T @ Gd) \
      - 6*my_einsum("ptrs,tp", C["qu1"], Gd @ Gd.conj().T) \
      + 6*my_einsum("ptrs,tp", C["qu1"], Gu @ Gu.conj().T) \
      + 2*GammaH*C["phiu"] \
      + Gammau @ C["phiu"] \
      + C["phiu"] @ Gammau

    #I3 #co
    Beta["phid"] = -1/6*XiB*gp**2*I3 \
      + 1/3*gp**2*C["phid"] \
      - 2/3*gp**2*(2*my_einsum("rstt", C["dd"]) \
      + 2/3*my_einsum("rtts", C["dd"]) \
      + my_einsum("ttrs", C["ed"]) \
      + my_einsum("ttrs", C["ld"]) \
      - my_einsum("ttrs", C["qd1"]) \
      - 2*my_einsum("ttrs", C["ud1"])) \
      + (C["phiBox"] \
      + C["phiD"])*Gd.conj().T @ Gd \
      - 2*Gd.conj().T @ C["phiq1"] @ Gd \
      + 3*(Gd.conj().T @ Gd @ C["phid"] \
      + C["phid"] @ Gd.conj().T @ Gd) \
      - Gd.conj().T @ Gu @ C["phiud"] \
      - C["phiud"].conj().T @ Gu.conj().T @ Gd \
      + 4*(3*my_einsum("rspt,tp", C["dd"], Gd.conj().T @ Gd) \
      + my_einsum("rtps,tp", C["dd"], Gd.conj().T @ Gd)) \
      + 2*my_einsum("ptrs,tp", C["ed"], Ge.conj().T @ Ge) \
      - 2*my_einsum("ptrs,tp", C["ld"], Ge @ Ge.conj().T) \
      - 6*my_einsum("ptrs,tp", C["ud1"], Gu.conj().T @ Gu) \
      - 6*my_einsum("ptrs,tp", C["qd1"], Gd @ Gd.conj().T) \
      + 6*my_einsum("ptrs,tp", C["qd1"], Gu @ Gu.conj().T) \
      + 2*GammaH*C["phid"] \
      + Gammad @ C["phid"] \
      + C["phid"] @ Gammad

        #co
    Beta["phiud"] = -3*gp**2*C["phiud"] \
      + (2*C["phiBox"] \
      - C["phiD"])*Gu.conj().T @ Gd \
      - 2*Gu.conj().T @ Gd @ C["phid"] \
      + 2*C["phiu"] @ Gu.conj().T @ Gd \
      + 4*(my_einsum("rtps,tp", C["ud1"], Gu.conj().T @ Gd) \
      + 4/3*my_einsum("rtps,tp", C["ud8"], Gu.conj().T @ Gd)) \
      + 2*Gu.conj().T @ Gu @ C["phiud"] \
      + 2*C["phiud"] @ Gd.conj().T @ Gd \
      + 2*GammaH*C["phiud"] \
      + Gammau @ C["phiud"] \
      + C["phiud"] @ Gammad

    """Dimension-5"""
    Beta["llphiphi"] = (2*Lambda \
      - 3*g**2 \
      + 2*GammaH)*C["llphiphi"]-3/2*(C["llphiphi"] @ Ge @ Ge.conj().T \
      + Ge.conj() @ Ge.T @ C["llphiphi"])

    """(3,3,3,3)"""
    # the einsum function is strong
    Beta["ll"] = -1/6*gp**2*my_einsum("st,pr", C["phil1"], I3) \
      - 1/6*g**2*(my_einsum("st,pr", C["phil3"], I3) \
      - 2*my_einsum("sr,pt", C["phil3"], I3)) \
      + 1/3*gp**2*(2*my_einsum("prww,st", C["ll"], I3) \
      + my_einsum("pwwr,st", C["ll"], I3)) \
      - 1/3*g**2*my_einsum("pwwr,st", C["ll"], I3) \
      + 2/3*g**2*my_einsum("swwr,pt", C["ll"], I3) \
      - 1/3*gp**2*my_einsum("prww,st", C["lq1"], I3) \
      - g**2*my_einsum("prww,st", C["lq3"], I3) \
      + 2*g**2*my_einsum("ptww,rs", C["lq3"], I3) \
      + 1/3*gp**2*( \
      - 2*my_einsum("prww,st", C["lu"], I3) \
      + my_einsum("prww,st", C["ld"], I3) \
      + my_einsum("prww,st", C["le"], I3)) \
      - 1/2*(my_einsum("pr,st", Ge @ Ge.conj().T, C["phil1"]) \
      - my_einsum("pr,st", Ge @ Ge.conj().T, C["phil3"])) \
      - my_einsum("pt,sr", Ge @ Ge.conj().T, C["phil3"]) \
      - 1/2*my_einsum("sv,tw,prvw", Ge, np.conj(Ge), C["le"]) \
      + my_einsum("pv,vrst", Gammal, C["ll"]) \
      + my_einsum("pvst,vr", C["ll"], Gammal) \
      - 1/6*gp**2*my_einsum("pr,st", C["phil1"], I3) \
      - 1/6*g**2*(my_einsum("pr,st", C["phil3"], I3) \
      - 2*my_einsum("pt,sr", C["phil3"], I3)) \
      + 1/3*gp**2*(2*my_einsum("stww,pr", C["ll"], I3) \
      + my_einsum("swwt,pr", C["ll"], I3)) \
      - 1/3*g**2*my_einsum("swwt,pr", C["ll"], I3) \
      + 2/3*g**2*my_einsum("pwwt,sr", C["ll"], I3) \
      - 1/3*gp**2*my_einsum("stww,pr", C["lq1"], I3) \
      - g**2*my_einsum("stww,pr", C["lq3"], I3) \
      + 2*g**2*my_einsum("srww,tp", C["lq3"], I3) \
      + 1/3*gp**2*( \
      - 2*my_einsum("stww,pr", C["lu"], I3) \
      + my_einsum("stww,pr", C["ld"], I3) \
      + my_einsum("stww,pr", C["le"], I3)) \
      - 1/2*(my_einsum("st,pr", Ge @ Ge.conj().T, C["phil1"]) \
      - my_einsum("st,pr", Ge @ Ge.conj().T, C["phil3"])) \
      - my_einsum("sr,pt", Ge @ Ge.conj().T, C["phil3"]) \
      - 1/2*my_einsum("pv,rw,stvw", Ge, np.conj(Ge), C["le"]) \
      + my_einsum("sv,vtpr", Gammal, C["ll"]) \
      + my_einsum("svpr,vt", C["ll"], Gammal) \
      + 6*g**2*my_einsum("ptsr", C["ll"]) \
      + 3*(gp**2 \
      - g**2)*my_einsum("prst", C["ll"])

    Beta["qq1"] = 1/18*gp**2*my_einsum("st,pr", C["phiq1"], I3) \
      - 1/9*gp**2*my_einsum("wwst,pr", C["lq1"], I3) \
      + 1/9*gp**2*(2*my_einsum("prww,st", C["qq1"], I3) \
      + 1/3*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3))) \
      + 1/3*gs**2*(my_einsum("swwr,pt", C["qq1"], I3) \
      + 3*my_einsum("swwr,pt", C["qq3"], I3)) \
      - 2/9*gs**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3)) \
      + 2/9*gp**2*my_einsum("prww,st", C["qu1"], I3) \
      - 1/9*gp**2*my_einsum("prww,st", C["qd1"], I3) \
      + 1/12*gs**2*(my_einsum("srww,pt", C["qu8"], I3) \
      + my_einsum("srww,pt", C["qd8"], I3)) \
      - 1/18*gs**2*(my_einsum("prww,st", C["qu8"], I3) \
      + my_einsum("prww,st", C["qd8"], I3)) \
      - 1/9*gp**2*my_einsum("prww,st", C["qe"], I3) \
      + 1/2*(my_einsum("pr,st", Gu @ Gu.conj().T, C["phiq1"]) \
      - my_einsum("pr,st", Gd @ Gd.conj().T, C["phiq1"])) \
      - 1/2*(my_einsum("pv,rw,stvw", Gu, np.conj(Gu), C["qu1"]) \
      - 1/6*my_einsum("pv,rw,stvw", Gu, np.conj(Gu), C["qu8"])) \
      - 1/2*(my_einsum("pv,rw,stvw", Gd, np.conj(Gd), C["qd1"]) \
      - 1/6*my_einsum("pv,rw,stvw", Gd, np.conj(Gd), C["qd8"])) \
      - 1/8*(my_einsum("pv,tw,srvw", Gu, np.conj(Gu), C["qu8"]) \
      + my_einsum("pv,tw,srvw", Gd, np.conj(Gd), C["qd8"])) \
      - 1/8*(my_einsum("tw,rv,pvsw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - 1/6*my_einsum("tw,rv,pvsw", np.conj(Gd), np.conj(Gu), C["quqd8"])) \
      - 1/8*(my_einsum("sw,pv,rvtw", Gd, Gu, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("sw,pv,rvtw", Gd, Gu, np.conj(C["quqd8"]))) \
      + 1/16*(my_einsum("tw,rv,svpw", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("sw,pv,tvrw", Gd, Gu, np.conj(C["quqd8"]))) \
      + my_einsum("pv,vrst", Gammaq, C["qq1"]) \
      + my_einsum("pvst,vr", C["qq1"], Gammaq) \
      + 1/18*gp**2*my_einsum("pr,st", C["phiq1"], I3) \
      - 1/9*gp**2*my_einsum("wwpr,st", C["lq1"], I3) \
      + 1/9*gp**2*(2*my_einsum("stww,pr", C["qq1"], I3) \
      + 1/3*(my_einsum("swwt,pr", C["qq1"], I3) \
      + 3*my_einsum("swwt,pr", C["qq3"], I3))) \
      + 1/3*gs**2*(my_einsum("pwwt,sr", C["qq1"], I3) \
      + 3*my_einsum("pwwt,sr", C["qq3"], I3)) \
      - 2/9*gs**2*(my_einsum("swwt,pr", C["qq1"], I3) \
      + 3*my_einsum("swwt,pr", C["qq3"], I3)) \
      + 2/9*gp**2*my_einsum("stww,pr", C["qu1"], I3) \
      - 1/9*gp**2*my_einsum("stww,pr", C["qd1"], I3) \
      + 1/12*gs**2*(my_einsum("ptww,sr", C["qu8"], I3) \
      + my_einsum("ptww,sr", C["qd8"], I3)) \
      - 1/18*gs**2*(my_einsum("stww,pr", C["qu8"], I3) \
      + my_einsum("stww,pr", C["qd8"], I3)) \
      - 1/9*gp**2*my_einsum("stww,pr", C["qe"], I3) \
      + 1/2*(my_einsum("st,pr", Gu @ Gu.conj().T, C["phiq1"]) \
      - my_einsum("st,pr", Gd @ Gd.conj().T, C["phiq1"])) \
      - 1/2*(my_einsum("sv,tw,prvw", Gu, np.conj(Gu), C["qu1"]) \
      - 1/6*my_einsum("sv,tw,prvw", Gu, np.conj(Gu), C["qu8"])) \
      - 1/2*(my_einsum("sv,tw,prvw", Gd, np.conj(Gd), C["qd1"]) \
      - 1/6*my_einsum("sv,tw,prvw", Gd, np.conj(Gd), C["qd8"])) \
      - 1/8*(my_einsum("sv,rw,ptvw", Gu, np.conj(Gu), C["qu8"]) \
      + my_einsum("sv,rw,ptvw", Gd, np.conj(Gd), C["qd8"])) \
      - 1/8*(my_einsum("rw,tv,svpw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - 1/6*my_einsum("rw,tv,svpw", np.conj(Gd), np.conj(Gu), C["quqd8"])) \
      - 1/8*(my_einsum("pw,sv,tvrw", Gd, Gu, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("pw,sv,tvrw", Gd, Gu, np.conj(C["quqd8"]))) \
      + 1/16*(my_einsum("rw,tv,pvsw", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("pw,sv,rvtw", Gd, Gu, np.conj(C["quqd8"]))) \
      + my_einsum("sv,vtpr", Gammaq, C["qq1"]) \
      + my_einsum("svpr,vt", C["qq1"], Gammaq) \
      + 9*g**2*my_einsum("prst", C["qq3"]) \
      - 2*(gs**2 \
      - 1/6*gp**2)*my_einsum("prst", C["qq1"]) \
      + 3*gs**2*(my_einsum("ptsr", C["qq1"]) \
      + 3*my_einsum("ptsr", C["qq3"]))

    Beta["qq3"] = 1/6*g**2*my_einsum("st,pr", C["phiq3"], I3) \
      + 1/3*g**2*my_einsum("wwst,pr", C["lq3"], I3) \
      + 1/3*g**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      - my_einsum("pwwr,st", C["qq3"], I3)) \
      + 2*g**2*my_einsum("prww,st", C["qq3"], I3) \
      + 1/3*gs**2*(my_einsum("swwr,pt", C["qq1"], I3) \
      + 3*my_einsum("swwr,pt", C["qq3"], I3)) \
      + 1/12*gs**2*(my_einsum("srww,pt", C["qu8"], I3) \
      + my_einsum("srww,pt", C["qd8"], I3)) \
      - 1/2*(my_einsum("pr,st", Gu @ Gu.conj().T, C["phiq3"]) \
      + my_einsum("pr,st", Gd @ Gd.conj().T, C["phiq3"])) \
      - 1/8*(my_einsum("pv,tw,srvw", Gu, np.conj(Gu), C["qu8"]) \
      + my_einsum("pv,tw,srvw", Gd, np.conj(Gd), C["qd8"])) \
      + 1/8*(my_einsum("tw,rv,pvsw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - 1/6*my_einsum("tw,rv,pvsw", np.conj(Gd), np.conj(Gu), C["quqd8"])) \
      + 1/8*(my_einsum("sw,pv,rvtw", Gd, Gu, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("sw,pv,rvtw", Gd, Gu, np.conj(C["quqd8"]))) \
      - 1/16*(my_einsum("tw,rv,svpw", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("sw,pv,tvrw", Gd, Gu, np.conj(C["quqd8"]))) \
      + my_einsum("pv,vrst", Gammaq, C["qq3"]) \
      + my_einsum("pvst,vr", C["qq3"], Gammaq) \
      + 1/6*g**2*my_einsum("pr,st", C["phiq3"], I3) \
      + 1/3*g**2*my_einsum("wwpr,st", C["lq3"], I3) \
      + 1/3*g**2*(my_einsum("swwt,pr", C["qq1"], I3) \
      - my_einsum("swwt,pr", C["qq3"], I3)) \
      + 2*g**2*my_einsum("stww,pr", C["qq3"], I3) \
      + 1/3*gs**2*(my_einsum("pwwt,sr", C["qq1"], I3) \
      + 3*my_einsum("pwwt,sr", C["qq3"], I3)) \
      + 1/12*gs**2*(my_einsum("ptww,sr", C["qu8"], I3) \
      + my_einsum("ptww,sr", C["qd8"], I3)) \
      - 1/2*(my_einsum("st,pr", Gu @ Gu.conj().T, C["phiq3"]) \
      + my_einsum("st,pr", Gd @ Gd.conj().T, C["phiq3"])) \
      - 1/8*(my_einsum("sv,rw,ptvw", Gu, np.conj(Gu), C["qu8"]) \
      + my_einsum("sv,rw,ptvw", Gd, np.conj(Gd), C["qd8"])) \
      + 1/8*(my_einsum("rw,tv,svpw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - 1/6*my_einsum("rw,tv,svpw", np.conj(Gd), np.conj(Gu), C["quqd8"])) \
      + 1/8*(my_einsum("pw,sv,tvrw", Gd, Gu, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("pw,sv,tvrw", Gd, Gu, np.conj(C["quqd8"]))) \
      - 1/16*(my_einsum("rw,tv,pvsw", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("pw,sv,rvtw", Gd, Gu, np.conj(C["quqd8"]))) \
      + my_einsum("sv,vtpr", Gammaq, C["qq3"]) \
      + my_einsum("svpr,vt", C["qq3"], Gammaq) \
      + 3*gs**2*(my_einsum("ptsr", C["qq1"]) \
      - my_einsum("ptsr", C["qq3"])) \
      - 2*(gs**2 \
      + 3*g**2 \
      - 1/6*gp**2)*my_einsum("prst", C["qq3"]) \
      + 3*g**2*my_einsum("prst", C["qq1"])

    #the terms are equal, but the order is not. No wonder if you check some differences inside
    Beta["lq1"] = -1/3*gp**2*my_einsum("st,pr", C["phiq1"], I3) \
      + 1/9*gp**2*my_einsum("pr,st", C["phil1"], I3) \
      - 2/9*gp**2*(2*my_einsum("prww,st", C["ll"], I3) \
      + my_einsum("pwwr,st", C["ll"], I3)) \
      + 2/9*gp**2*my_einsum("prww,st", C["lq1"], I3) \
      + 2/3*gp**2*my_einsum("wwst,pr", C["lq1"], I3) \
      - 2/9*gp**2*(6*my_einsum("stww,pr", C["qq1"], I3) \
      + my_einsum("swwt,pr", C["qq1"], I3) \
      + 3*my_einsum("swwt,pr", C["qq3"], I3)) \
      - 2/3*gp**2*(2*my_einsum("stww,pr", C["qu1"], I3) \
      - my_einsum("stww,pr", C["qd1"], I3) \
      - my_einsum("stww,pr", C["qe"], I3)) \
      + 2/9*gp**2*(2*my_einsum("prww,st", C["lu"], I3) \
      - my_einsum("prww,st", C["ld"], I3) \
      - my_einsum("prww,st", C["le"], I3)) \
      - gp**2*my_einsum("prst", C["lq1"]) \
      + 9*g**2*my_einsum("prst", C["lq3"]) \
      - my_einsum("pr,st", Ge @ Ge.conj().T, C["phiq1"]) \
      + my_einsum("st,pr", Gu @ Gu.conj().T, C["phil1"]) \
      - my_einsum("st,pr", Gd @ Gd.conj().T, C["phil1"]) \
      + 1/4*(my_einsum("tw,rv,pvsw", np.conj(Gu), np.conj(Ge), C["lequ1"]) \
      - 12*my_einsum("tw,rv,pvsw", np.conj(Gu), np.conj(Ge), C["lequ3"]) \
      + my_einsum("sw,pv,rvtw", Gu, Ge, np.conj(C["lequ1"])) \
      - 12*my_einsum("sw,pv,rvtw", Gu, Ge, np.conj(C["lequ3"]))) \
      - my_einsum("sv,tw,prvw", Gu, np.conj(Gu), C["lu"]) \
      - my_einsum("sv,tw,prvw", Gd, np.conj(Gd), C["ld"]) \
      - my_einsum("pv,rw,stvw", Ge, np.conj(Ge), C["qe"]) \
      + 1/4*(my_einsum("sw,rv,pvwt", Gd, np.conj(Ge), C["ledq"]) \
      + my_einsum("pv,tw,rvws", Ge, np.conj(Gd), np.conj(C["ledq"]))) \
      + my_einsum("pv,vrst", Gammal, C["lq1"]) \
      + my_einsum("sv,prvt", Gammaq, C["lq1"]) \
      + my_einsum("pvst,vr", C["lq1"], Gammal) \
      + my_einsum("prsv,vt", C["lq1"], Gammaq)

    Beta["lq3"] = 1/3*g**2*(my_einsum("st,pr", C["phiq3"], I3) \
      + my_einsum("pr,st", C["phil3"], I3)) \
      + 2/3*g**2*(3*my_einsum("prww,st", C["lq3"], I3) \
      + my_einsum("wwst,pr", C["lq3"], I3)) \
      + 2/3*g**2*(6*my_einsum("stww,pr", C["qq3"], I3) \
      + my_einsum("swwt,pr", C["qq1"], I3) \
      - my_einsum("swwt,pr", C["qq3"], I3)) \
      + 2/3*g**2*my_einsum("pwwr,st", C["ll"], I3) \
      + 3*g**2*my_einsum("prst", C["lq1"]) \
      - (6*g**2 \
      + gp**2)*my_einsum("prst", C["lq3"]) \
      - my_einsum("pr,st", Ge @ Ge.conj().T, C["phiq3"]) \
      - my_einsum("st,pr", Gu @ Gu.conj().T, C["phil3"]) \
      - my_einsum("st,pr", Gd @ Gd.conj().T, C["phil3"]) \
      - 1/4*(my_einsum("tw,rv,pvsw", np.conj(Gu), np.conj(Ge), C["lequ1"]) \
      - 12*my_einsum("tw,rv,pvsw", np.conj(Gu), np.conj(Ge), C["lequ3"]) \
      + my_einsum("sw,pv,rvtw", Gu, Ge, np.conj(C["lequ1"])) \
      - 12*my_einsum("sw,pv,rvtw", Gu, Ge, np.conj(C["lequ3"]))) \
      + 1/4*(my_einsum("sw,rv,pvwt", Gd, np.conj(Ge), C["ledq"]) \
      + my_einsum("pv,tw,rvws", Ge, np.conj(Gd), np.conj(C["ledq"]))) \
      + my_einsum("pv,vrst", Gammal, C["lq3"]) \
      + my_einsum("sv,prvt", Gammaq, C["lq3"]) \
      + my_einsum("pvst,vr", C["lq3"], Gammal) \
      + my_einsum("prsv,vt", C["lq3"], Gammaq)

    #order
    Beta["ee"] = -1/3*gp**2*my_einsum("st,pr", C["phie"], I3) \
      + 2/3*gp**2*(my_einsum("wwpr,st", C["le"], I3) \
      - my_einsum("wwpr,st", C["qe"], I3) \
      - 2*my_einsum("prww,st", C["eu"], I3) \
      + my_einsum("prww,st", C["ed"], I3) \
      + 4*my_einsum("prww,st", C["ee"], I3)) \
      + my_einsum("pr,st", Ge.conj().T @ Ge, C["phie"]) \
      - my_einsum("wr,vp,vwst", Ge, np.conj(Ge), C["le"]) \
      + my_einsum("pv,vrst", Gammae, C["ee"]) \
      + my_einsum("pvst,vr", C["ee"], Gammae) \
      - 1/3*gp**2*my_einsum("pr,st", C["phie"], I3) \
      + 2/3*gp**2*(my_einsum("wwst,pr", C["le"], I3) \
      - my_einsum("wwst,pr", C["qe"], I3) \
      - 2*my_einsum("stww,pr", C["eu"], I3) \
      + my_einsum("stww,pr", C["ed"], I3) \
      + 4*my_einsum("wwst,pr", C["ee"], I3)) \
      + my_einsum("st,pr", Ge.conj().T @ Ge, C["phie"]) \
      - my_einsum("wt,vs,vwpr", Ge, np.conj(Ge), C["le"]) \
      + my_einsum("sv,vtpr", Gammae, C["ee"]) \
      + my_einsum("svpr,vt", C["ee"], Gammae) \
      + 12*gp**2*my_einsum("prst", C["ee"])

    #order
    Beta["uu"] = 2/9*gp**2*my_einsum("st,pr", C["phiu"], I3) \
      - 4/9*gp**2*(my_einsum("wwst,pr", C["eu"], I3) \
      + my_einsum("wwst,pr", C["lu"], I3) \
      - my_einsum("wwst,pr", C["qu1"], I3) \
      - 4*my_einsum("wwst,pr", C["uu"], I3) \
      - 4/3*my_einsum("swwt,pr", C["uu"], I3)) \
      - 1/9*gs**2*(my_einsum("wwst,pr", C["qu8"], I3) \
      - 3*my_einsum("wwsr,pt", C["qu8"], I3)) \
      + 2/3*gs**2*my_einsum("pwwt,rs", C["uu"], I3) \
      - 2/9*gs**2*my_einsum("swwt,pr", C["uu"], I3) \
      - 4/9*gp**2*my_einsum("stww,pr", C["ud1"], I3) \
      - 1/18*gs**2*(my_einsum("stww,pr", C["ud8"], I3) \
      - 3*my_einsum("srww,pt", C["ud8"], I3)) \
      - my_einsum("pr,st", Gu.conj().T @ Gu, C["phiu"]) \
      - (my_einsum("wr,vp,vwst", Gu, np.conj(Gu), C["qu1"]) \
      - 1/6*my_einsum("wr,vp,vwst", Gu, np.conj(Gu), C["qu8"])) \
      - 1/2*my_einsum("wr,vs,vwpt", Gu, np.conj(Gu), C["qu8"]) \
      + my_einsum("pv,vrst", Gammau, C["uu"]) \
      + my_einsum("pvst,vr", C["uu"], Gammau) \
      + 2/9*gp**2*my_einsum("pr,st", C["phiu"], I3) \
      - 4/9*gp**2*(my_einsum("wwpr,st", C["eu"], I3) \
      + my_einsum("wwpr,st", C["lu"], I3) \
      - my_einsum("wwpr,st", C["qu1"], I3) \
      - 4*my_einsum("wwpr,st", C["uu"], I3) \
      - 4/3*my_einsum("pwwr,st", C["uu"], I3)) \
      - 1/9*gs**2*(my_einsum("wwpr,st", C["qu8"], I3) \
      - 3*my_einsum("wwpt,sr", C["qu8"], I3)) \
      + 2/3*gs**2*my_einsum("swwr,tp", C["uu"], I3) \
      - 2/9*gs**2*my_einsum("pwwr,st", C["uu"], I3) \
      - 4/9*gp**2*my_einsum("prww,st", C["ud1"], I3) \
      - 1/18*gs**2*(my_einsum("prww,st", C["ud8"], I3) \
      - 3*my_einsum("ptww,sr", C["ud8"], I3)) \
      - my_einsum("st,pr", Gu.conj().T @ Gu, C["phiu"]) \
      - (my_einsum("wt,vs,vwpr", Gu, np.conj(Gu), C["qu1"]) \
      - 1/6*my_einsum("wt,vs,vwpr", Gu, np.conj(Gu), C["qu8"])) \
      - 1/2*my_einsum("wt,vp,vwsr", Gu, np.conj(Gu), C["qu8"]) \
      + my_einsum("sv,vtpr", Gammau, C["uu"]) \
      + my_einsum("svpr,vt", C["uu"], Gammau) \
      + 2*(8/3*gp**2 \
      - gs**2)*my_einsum("prst", C["uu"]) \
      + 6*gs**2*my_einsum("ptsr", C["uu"])

    #order
    Beta["dd"] = -1/9*gp**2*my_einsum("st,pr", C["phid"], I3) \
      + 2/9*gp**2*(my_einsum("wwst,pr", C["ed"], I3) \
      + my_einsum("wwst,pr", C["ld"], I3) \
      - my_einsum("wwst,pr", C["qd1"], I3) \
      + 2*my_einsum("wwst,pr", C["dd"], I3) \
      + 2/3*my_einsum("swwt,pr", C["dd"], I3)) \
      - 1/9*gs**2*(my_einsum("wwst,pr", C["qd8"], I3) \
      - 3*my_einsum("wwsr,pt", C["qd8"], I3)) \
      + 2/3*gs**2*my_einsum("pwwt,rs", C["dd"], I3) \
      - 2/9*gs**2*my_einsum("swwt,pr", C["dd"], I3) \
      - 4/9*gp**2*my_einsum("wwst,pr", C["ud1"], I3) \
      - 1/18*gs**2*(my_einsum("wwst,pr", C["ud8"], I3) \
      - 3*my_einsum("wwsr,pt", C["ud8"], I3)) \
      + my_einsum("pr,st", Gd.conj().T @ Gd, C["phid"]) \
      - (my_einsum("wr,vp,vwst", Gd, np.conj(Gd), C["qd1"]) \
      - 1/6*my_einsum("wr,vp,vwst", Gd, np.conj(Gd), C["qd8"])) \
      - 1/2*my_einsum("wr,vs,vwpt", Gd, np.conj(Gd), C["qd8"]) \
      + my_einsum("pv,vrst", Gammad, C["dd"]) \
      + my_einsum("pvst,vr", C["dd"], Gammad) \
      - 1/9*gp**2*my_einsum("pr,st", C["phid"], I3) \
      + 2/9*gp**2*(my_einsum("wwpr,st", C["ed"], I3) \
      + my_einsum("wwpr,st", C["ld"], I3) \
      - my_einsum("wwpr,st", C["qd1"], I3) \
      + 2*my_einsum("wwpr,st", C["dd"], I3) \
      + 2/3*my_einsum("pwwr,st", C["dd"], I3)) \
      - 1/9*gs**2*(my_einsum("wwpr,st", C["qd8"], I3) \
      - 3*my_einsum("wwpt,sr", C["qd8"], I3)) \
      + 2/3*gs**2*my_einsum("swwr,tp", C["dd"], I3) \
      - 2/9*gs**2*my_einsum("pwwr,st", C["dd"], I3) \
      - 4/9*gp**2*my_einsum("wwpr,st", C["ud1"], I3) \
      - 1/18*gs**2*(my_einsum("wwpr,st", C["ud8"], I3) \
      - 3*my_einsum("wwpt,sr", C["ud8"], I3)) \
      + my_einsum("st,pr", Gd.conj().T @ Gd, C["phid"]) \
      - (my_einsum("wt,vs,vwpr", Gd, np.conj(Gd), C["qd1"]) \
      - 1/6*my_einsum("wt,vs,vwpr", Gd, np.conj(Gd), C["qd8"])) \
      - 1/2*my_einsum("wt,vp,vwsr", Gd, np.conj(Gd), C["qd8"]) \
      + my_einsum("sv,vtpr", Gammad, C["dd"]) \
      + my_einsum("svpr,vt", C["dd"], Gammad) \
      + 2*(2/3*gp**2 \
      - gs**2)*my_einsum("prst", C["dd"]) \
      + 6*gs**2*my_einsum("ptsr", C["dd"])

    Beta["eu"] = -2/3*gp**2*(my_einsum("st,pr", C["phiu"], I3) \
      + 2*(my_einsum("wwst,pr", C["qu1"], I3) \
      - my_einsum("wwst,pr", C["lu"], I3) \
      + 4*my_einsum("wwst,pr", C["uu"], I3) \
      - my_einsum("wwst,pr", C["eu"], I3) \
      - my_einsum("stww,pr", C["ud1"], I3)) \
      + 8/3*my_einsum("swwt,pr", C["uu"], I3)) \
      + 4/9*gp**2*(my_einsum("pr,st", C["phie"], I3) \
      + 2*(my_einsum("wwpr,st", C["qe"], I3) \
      - my_einsum("wwpr,st", C["le"], I3) \
      - 4*my_einsum("prww,st", C["ee"], I3) \
      + 2*my_einsum("prww,st", C["eu"], I3) \
      - my_einsum("prww,st", C["ed"], I3))) \
      - 8*gp**2*my_einsum("prst", C["eu"]) \
      + 2*my_einsum("pr,st", Ge.conj().T @ Ge, C["phiu"]) \
      - 2*my_einsum("st,pr", Gu.conj().T @ Gu, C["phie"]) \
      + my_einsum("vp,ws,vrwt", np.conj(Ge), np.conj(Gu), C["lequ1"]) \
      - 12*my_einsum("vp,ws,vrwt", np.conj(Ge), np.conj(Gu), C["lequ3"]) \
      + my_einsum("vr,wt,vpws", Ge, Gu, np.conj(C["lequ1"])) \
      - 12*my_einsum("vr,wt,vpws", Ge, Gu, np.conj(C["lequ3"])) \
      - 2*my_einsum("vp,wr,vwst", np.conj(Ge), Ge, C["lu"]) \
      - 2*my_einsum("vs,wt,vwpr", np.conj(Gu), Gu, C["qe"]) \
      + my_einsum("pv,vrst", Gammae, C["eu"]) \
      + my_einsum("sv,prvt", Gammau, C["eu"]) \
      + my_einsum("pvst,vr", C["eu"], Gammae) \
      + my_einsum("prsv,vt", C["eu"], Gammau)

    Beta["ed"] = -2/3*gp**2*(my_einsum("st,pr", C["phid"], I3) \
      + 2*(my_einsum("wwst,pr", C["qd1"], I3) \
      - my_einsum("wwst,pr", C["ld"], I3) \
      - 2*my_einsum("wwst,pr", C["dd"], I3) \
      - my_einsum("wwst,pr", C["ed"], I3) \
      + 2*my_einsum("wwst,pr", C["ud1"], I3)) \
      - 4/3*my_einsum("swwt,pr", C["dd"], I3)) \
      - 2/9*gp**2*(my_einsum("pr,st", C["phie"], I3) \
      + 2*(my_einsum("wwpr,st", C["qe"], I3) \
      - my_einsum("wwpr,st", C["le"], I3) \
      - 4*my_einsum("prww,st", C["ee"], I3) \
      - my_einsum("prww,st", C["ed"], I3) \
      + 2*my_einsum("prww,st", C["eu"], I3))) \
      + 4*gp**2*my_einsum("prst", C["ed"]) \
      + 2*my_einsum("pr,st", Ge.conj().T @ Ge, C["phid"]) \
      + 2*my_einsum("st,pr", Gd.conj().T @ Gd, C["phie"]) \
      - 2*my_einsum("vp,wr,vwst", np.conj(Ge), Ge, C["ld"]) \
      - 2*my_einsum("vs,wt,vwpr", np.conj(Gd), Gd, C["qe"]) \
      + my_einsum("vp,wt,vrsw", np.conj(Ge), Gd, C["ledq"]) \
      + my_einsum("vr,ws,vptw", Ge, np.conj(Gd), np.conj(C["ledq"])) \
      + my_einsum("pv,vrst", Gammae, C["ed"]) \
      + my_einsum("sv,prvt", Gammad, C["ed"]) \
      + my_einsum("pvst,vr", C["ed"], Gammae) \
      + my_einsum("prsv,vt", C["ed"], Gammad)

    #order
    Beta["ud1"] = 4/9*gp**2*(my_einsum("st,pr", C["phid"], I3) \
      + 2*(my_einsum("wwst,pr", C["qd1"], I3) \
      - my_einsum("wwst,pr", C["ld"], I3) \
      - 2*my_einsum("wwst,pr", C["dd"], I3) \
      + 2*my_einsum("wwst,pr", C["ud1"], I3) \
      - my_einsum("wwst,pr", C["ed"], I3)) \
      - 4/3*my_einsum("swwt,pr", C["dd"], I3)) \
      - 2/9*gp**2*(my_einsum("pr,st", C["phiu"], I3) \
      + 2*(my_einsum("wwpr,st", C["qu1"], I3) \
      - my_einsum("wwpr,st", C["lu"], I3) \
      + 4*my_einsum("wwpr,st", C["uu"], I3) \
      - my_einsum("prww,st", C["ud1"], I3) \
      - my_einsum("wwpr,st", C["eu"], I3)) \
      + 8/3*my_einsum("pwwr,st", C["uu"], I3)) \
      - 8/3*(gp**2*my_einsum("prst", C["ud1"]) \
      - gs**2*my_einsum("prst", C["ud8"])) \
      - 2*my_einsum("pr,st", Gu.conj().T @ Gu, C["phid"]) \
      + 2*my_einsum("st,pr", Gd.conj().T @ Gd, C["phiu"]) \
      + 2/3*my_einsum("sr,pt", Gd.conj().T @ Gu, C["phiud"]) \
      + 2/3*my_einsum("pt,rs", Gu.conj().T @ Gd, np.conj(C["phiud"])) \
      + 1/3*(my_einsum("vs,wp,vrwt", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      + 4/3*my_einsum("vs,wp,vrwt", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("vt,wr,vpws", Gd, Gu, np.conj(C["quqd1"])) \
      + 4/3*my_einsum("vt,wr,vpws", Gd, Gu, np.conj(C["quqd8"]))) \
      - my_einsum("ws,vp,vrwt", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - my_einsum("wt,vr,vpws", Gd, Gu, np.conj(C["quqd1"])) \
      - 2*my_einsum("vp,wr,vwst", np.conj(Gu), Gu, C["qd1"]) \
      - 2*my_einsum("vs,wt,vwpr", np.conj(Gd), Gd, C["qu1"]) \
      + my_einsum("pv,vrst", Gammau, C["ud1"]) \
      + my_einsum("sv,prvt", Gammad, C["ud1"]) \
      + my_einsum("pvst,vr", C["ud1"], Gammau) \
      + my_einsum("prsv,vt", C["ud1"], Gammad)

    #order
    Beta["ud8"] = 8/3*gs**2*my_einsum("pwwr,st", C["uu"], I3) \
      + 8/3*gs**2*my_einsum("swwt,pr", C["dd"], I3) \
      + 4/3*gs**2*my_einsum("wwpr,st", C["qu8"], I3) \
      + 4/3*gs**2*my_einsum("wwst,pr", C["qd8"], I3) \
      + 2/3*gs**2*my_einsum("prww,st", C["ud8"], I3) \
      + 2/3*gs**2*my_einsum("wwst,pr", C["ud8"], I3) \
      - 4*(2/3*gp**2 \
      + gs**2)*my_einsum("prst", C["ud8"]) \
      + 12*gs**2*my_einsum("prst", C["ud1"]) \
      + 4*my_einsum("sr,pt", Gd.conj().T @ Gu, C["phiud"]) \
      + 4*my_einsum("pt,rs", Gu.conj().T @ Gd, np.conj(C["phiud"])) \
      + 2*(my_einsum("vs,wp,vrwt", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - 1/6*my_einsum("vs,wp,vrwt", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("vt,wr,vpws", Gd, Gu, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("vt,wr,vpws", Gd, Gu, np.conj(C["quqd8"]))) \
      - 2*my_einsum("vp,wr,vwst", np.conj(Gu), Gu, C["qd8"]) \
      - 2*my_einsum("vs,wt,vwpr", np.conj(Gd), Gd, C["qu8"]) \
      - (my_einsum("ws,vp,vrwt", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + my_einsum("wt,vr,vpws", Gd, Gu, np.conj(C["quqd8"]))) \
      + my_einsum("pv,vrst", Gammau, C["ud8"]) \
      + my_einsum("sv,prvt", Gammad, C["ud8"]) \
      + my_einsum("pvst,vr", C["ud8"], Gammau) \
      + my_einsum("prsv,vt", C["ud8"], Gammad)

    Beta["le"] = -1/3*gp**2*my_einsum("st,pr", C["phie"], I3) \
      - 2/3*gp**2*my_einsum("pr,st", C["phil1"], I3) \
      + 8/3*gp**2*my_einsum("prww,st", C["ll"], I3) \
      + 4/3*gp**2*my_einsum("pwwr,st", C["ll"], I3) \
      - 4/3*gp**2*my_einsum("prww,st", C["lq1"], I3) \
      - 2/3*gp**2*my_einsum("wwst,pr", C["qe"], I3) \
      + 4/3*gp**2*my_einsum("prww,st", C["le"], I3) \
      + 2/3*gp**2*my_einsum("wwst,pr", C["le"], I3) \
      - 8/3*gp**2*my_einsum("prww,st", C["lu"], I3) \
      + 4/3*gp**2*my_einsum("prww,st", C["ld"], I3) \
      - 4/3*gp**2*my_einsum("stww,pr", C["eu"], I3) \
      + 2/3*gp**2*my_einsum("stww,pr", C["ed"], I3) \
      + 8/3*gp**2*my_einsum("wwst,pr", C["ee"], I3) \
      - 6*gp**2*my_einsum("prst", C["le"]) \
      + my_einsum("rs,pt", np.conj(Ge), Xie) \
      + my_einsum("pt,rs", Ge, np.conj(Xie)) \
      - my_einsum("pr,st", Ge @ Ge.conj().T, C["phie"]) \
      + 2*my_einsum("st,pr", Ge.conj().T @ Ge, C["phil1"]) \
      - 4*my_einsum("pv,rw,vtsw", Ge, np.conj(Ge), C["ee"]) \
      + my_einsum("pw,vs,vrwt", Ge, np.conj(Ge), C["le"]) \
      - 2*my_einsum("wt,vs,pwvr", Ge, np.conj(Ge), C["ll"]) \
      - 4*my_einsum("wt,vs,prvw", Ge, np.conj(Ge), C["ll"]) \
      + my_einsum("vt,rw,pvsw", Ge, np.conj(Ge), C["le"]) \
      + my_einsum("pv,vrst", Gammal, C["le"]) \
      + my_einsum("sv,prvt", Gammae, C["le"]) \
      + my_einsum("pvst,vr", C["le"], Gammal) \
      + my_einsum("prsv,vt", C["le"], Gammae)

    #order
    Beta["lu"] = -1/3*gp**2*my_einsum("st,pr", C["phiu"], I3) \
      + 4/9*gp**2*my_einsum("pr,st", C["phil1"], I3) \
      - 16/9*gp**2*my_einsum("prww,st", C["ll"], I3) \
      - 8/9*gp**2*my_einsum("pwwr,st", C["ll"], I3) \
      + 8/9*gp**2*my_einsum("prww,st", C["lq1"], I3) \
      - 2/3*gp**2*my_einsum("wwst,pr", C["qu1"], I3) \
      + 16/9*gp**2*my_einsum("prww,st", C["lu"], I3) \
      + 2/3*gp**2*my_einsum("wwst,pr", C["lu"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["ld"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["le"], I3) \
      + 2/3*gp**2*my_einsum("stww,pr", C["ud1"], I3) \
      + 2/3*gp**2*my_einsum("wwst,pr", C["eu"], I3) \
      - 8/3*gp**2*my_einsum("stww,pr", C["uu"], I3) \
      - 8/9*gp**2*my_einsum("swwt,pr", C["uu"], I3) \
      + 4*gp**2*my_einsum("prst", C["lu"]) \
      - my_einsum("pr,st", Ge @ Ge.conj().T, C["phiu"]) \
      - 2*my_einsum("st,pr", Gu.conj().T @ Gu, C["phil1"]) \
      - 1/2*(my_einsum("rv,ws,pvwt", np.conj(Ge), np.conj(Gu), C["lequ1"]) \
      + 12*my_einsum("rv,ws,pvwt", np.conj(Ge), np.conj(Gu), C["lequ3"])) \
      - 1/2*(my_einsum("pv,wt,rvws", Ge, Gu, np.conj(C["lequ1"])) \
      + 12*my_einsum("pv,wt,rvws", Ge, Gu, np.conj(C["lequ3"]))) \
      - 2*my_einsum("vs,wt,prvw", np.conj(Gu), Gu, C["lq1"]) \
      - my_einsum("rw,pv,vwst", np.conj(Ge), Ge, C["eu"]) \
      + my_einsum("pv,vrst", Gammal, C["lu"]) \
      + my_einsum("sv,prvt", Gammau, C["lu"]) \
      + my_einsum("pvst,vr", C["lu"], Gammal) \
      + my_einsum("prsv,vt", C["lu"], Gammau)

    Beta["ld"] = -1/3*gp**2*my_einsum("st,pr", C["phid"], I3) \
      - 2/9*gp**2*my_einsum("pr,st", C["phil1"], I3) \
      + 8/9*gp**2*my_einsum("prww,st", C["ll"], I3) \
      + 4/9*gp**2*my_einsum("pwwr,st", C["ll"], I3) \
      - 4/9*gp**2*my_einsum("prww,st", C["lq1"], I3) \
      - 2/3*gp**2*my_einsum("wwst,pr", C["qd1"], I3) \
      + 4/9*gp**2*my_einsum("prww,st", C["ld"], I3) \
      + 2/3*gp**2*my_einsum("wwst,pr", C["ld"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["lu"], I3) \
      + 4/9*gp**2*my_einsum("prww,st", C["le"], I3) \
      - 4/3*gp**2*my_einsum("wwst,pr", C["ud1"], I3) \
      + 2/3*gp**2*my_einsum("wwst,pr", C["ed"], I3) \
      + 4/3*gp**2*my_einsum("stww,pr", C["dd"], I3) \
      + 4/9*gp**2*my_einsum("swwt,pr", C["dd"], I3) \
      - 2*gp**2*my_einsum("prst", C["ld"]) \
      - my_einsum("pr,st", Ge @ Ge.conj().T, C["phid"]) \
      + 2*my_einsum("st,pr", Gd.conj().T @ Gd, C["phil1"]) \
      - 1/2*my_einsum("rv,wt,pvsw", np.conj(Ge), Gd, C["ledq"]) \
      - 1/2*my_einsum("pv,ws,rvtw", Ge, np.conj(Gd), np.conj(C["ledq"])) \
      - 2*my_einsum("vs,wt,prvw", np.conj(Gd), Gd, C["lq1"]) \
      - my_einsum("rw,pv,vwst", np.conj(Ge), Ge, C["ed"]) \
      + my_einsum("pv,vrst", Gammal, C["ld"]) \
      + my_einsum("sv,prvt", Gammad, C["ld"]) \
      + my_einsum("pvst,vr", C["ld"], Gammal) \
      + my_einsum("prsv,vt", C["ld"], Gammad)

    Beta["qe"] = 1/9*gp**2*my_einsum("st,pr", C["phie"], I3) \
      - 2/3*gp**2*my_einsum("pr,st", C["phiq1"], I3) \
      - 8/3*gp**2*my_einsum("prww,st", C["qq1"], I3) \
      - 4/9*gp**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3)) \
      + 4/3*gp**2*my_einsum("wwpr,st", C["lq1"], I3) \
      - 2/9*gp**2*my_einsum("wwst,pr", C["le"], I3) \
      + 4/3*gp**2*my_einsum("prww,st", C["qe"], I3) \
      + 2/9*gp**2*my_einsum("wwst,pr", C["qe"], I3) \
      - 8/3*gp**2*my_einsum("prww,st", C["qu1"], I3) \
      + 4/3*gp**2*my_einsum("prww,st", C["qd1"], I3) \
      + 4/9*gp**2*my_einsum("stww,pr", C["eu"], I3) \
      - 2/9*gp**2*my_einsum("stww,pr", C["ed"], I3) \
      - 8/9*gp**2*my_einsum("wwst,pr", C["ee"], I3) \
      + 2*gp**2*my_einsum("prst", C["qe"]) \
      + my_einsum("pr,st", Gu @ Gu.conj().T, C["phie"]) \
      - my_einsum("pr,st", Gd @ Gd.conj().T, C["phie"]) \
      + 2*my_einsum("st,pr", Ge.conj().T @ Ge, C["phiq1"]) \
      - 1/2*my_einsum("pw,vs,vtwr", Gd, np.conj(Ge), C["ledq"]) \
      - 1/2*my_einsum("vt,rw,vswp", Ge, np.conj(Gd), np.conj(C["ledq"])) \
      - 2*my_einsum("vs,wt,vwpr", np.conj(Ge), Ge, C["lq1"]) \
      - 1/2*(my_einsum("rw,vs,vtpw", np.conj(Gu), np.conj(Ge), C["lequ1"]) \
      + 12*my_einsum("rw,vs,vtpw", np.conj(Gu), np.conj(Ge), C["lequ3"])) \
      - 1/2*(my_einsum("pw,vt,vsrw", Gu, Ge, np.conj(C["lequ1"])) \
      + 12*my_einsum("pw,vt,vsrw", Gu, Ge, np.conj(C["lequ3"]))) \
      - my_einsum("rw,pv,stvw", np.conj(Gd), Gd, C["ed"]) \
      - my_einsum("rw,pv,stvw", np.conj(Gu), Gu, C["eu"]) \
      + my_einsum("pv,vrst", Gammaq, C["qe"]) \
      + my_einsum("sv,prvt", Gammae, C["qe"]) \
      + my_einsum("pvst,vr", C["qe"], Gammaq) \
      + my_einsum("prsv,vt", C["qe"], Gammae)

    Beta["qu1"] = 1/9*gp**2*my_einsum("st,pr", C["phiu"], I3) \
      + 4/9*gp**2*my_einsum("pr,st", C["phiq1"], I3) \
      + 16/9*gp**2*my_einsum("prww,st", C["qq1"], I3) \
      + 8/27*gp**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3)) \
      - 8/9*gp**2*my_einsum("wwpr,st", C["lq1"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["qe"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["qd1"], I3) \
      + 16/9*gp**2*my_einsum("prww,st", C["qu1"], I3) \
      + 2/9*gp**2*my_einsum("wwst,pr", C["qu1"], I3) \
      - 2/9*gp**2*my_einsum("wwst,pr", C["lu"], I3) \
      - 2/9*gp**2*my_einsum("wwst,pr", C["eu"], I3) \
      - 2/9*gp**2*my_einsum("stww,pr", C["ud1"], I3) \
      + 8/9*gp**2*my_einsum("stww,pr", C["uu"], I3) \
      + 8/27*gp**2*my_einsum("swwt,pr", C["uu"], I3) \
      - 4/3*gp**2*my_einsum("prst", C["qu1"]) \
      - 8/3*gs**2*my_einsum("prst", C["qu8"]) \
      + 1/3*my_einsum("rs,pt", np.conj(Gu), Xiu) \
      + 1/3*my_einsum("pt,rs", Gu, np.conj(Xiu)) \
      + my_einsum("pr,st", Gu @ Gu.conj().T, C["phiu"]) \
      - my_einsum("pr,st", Gd @ Gd.conj().T, C["phiu"]) \
      - 2*my_einsum("st,pr", Gu.conj().T @ Gu, C["phiq1"]) \
      + 1/3*(my_einsum("pw,vs,vrwt", Gu, np.conj(Gu), C["qu1"]) \
      + 4/3*my_einsum("pw,vs,vrwt", Gu, np.conj(Gu), C["qu8"])) \
      + 1/3*(my_einsum("vt,rw,pvsw", Gu, np.conj(Gu), C["qu1"]) \
      + 4/3*my_einsum("vt,rw,pvsw", Gu, np.conj(Gu), C["qu8"])) \
      + 1/3*(my_einsum("rw,vs,ptvw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      + 4/3*my_einsum("rw,vs,ptvw", np.conj(Gd), np.conj(Gu), C["quqd8"])) \
      + 1/3*(my_einsum("pw,vt,rsvw", Gd, Gu, np.conj(C["quqd1"])) \
      + 4/3*my_einsum("pw,vt,rsvw", Gd, Gu, np.conj(C["quqd8"]))) \
      + 1/2*my_einsum("rw,vs,vtpw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      + 1/2*my_einsum("pw,vt,vsrw", Gd, Gu, np.conj(C["quqd1"])) \
      - 2/3*(my_einsum("vt,ws,pvwr", Gu, np.conj(Gu), C["qq1"]) \
      + 3*my_einsum("vt,ws,pvwr", Gu, np.conj(Gu), C["qq3"])) \
      - 4*my_einsum("wt,vs,prvw", Gu, np.conj(Gu), C["qq1"]) \
      - 2/3*my_einsum("pv,rw,vtsw", Gu, np.conj(Gu), C["uu"]) \
      - 2*my_einsum("pv,rw,vwst", Gu, np.conj(Gu), C["uu"]) \
      - my_einsum("pv,rw,stvw", Gd, np.conj(Gd), C["ud1"]) \
      + my_einsum("pv,vrst", Gammaq, C["qu1"]) \
      + my_einsum("sv,prvt", Gammau, C["qu1"]) \
      + my_einsum("pvst,vr", C["qu1"], Gammaq) \
      + my_einsum("prsv,vt", C["qu1"], Gammau)

    Beta["qd1"] = 1/9*gp**2*my_einsum("st,pr", C["phid"], I3) \
      - 2/9*gp**2*my_einsum("pr,st", C["phiq1"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["qq1"], I3) \
      - 4/27*gp**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3)) \
      + 4/9*gp**2*my_einsum("wwpr,st", C["lq1"], I3) \
      + 4/9*gp**2*my_einsum("prww,st", C["qe"], I3) \
      - 8/9*gp**2*my_einsum("prww,st", C["qu1"], I3) \
      + 4/9*gp**2*my_einsum("prww,st", C["qd1"], I3) \
      + 2/9*gp**2*my_einsum("wwst,pr", C["qd1"], I3) \
      - 2/9*gp**2*my_einsum("wwst,pr", C["ld"], I3) \
      - 2/9*gp**2*my_einsum("wwst,pr", C["ed"], I3) \
      + 4/9*gp**2*my_einsum("wwst,pr", C["ud1"], I3) \
      - 4/9*gp**2*my_einsum("stww,pr", C["dd"], I3) \
      - 4/27*gp**2*my_einsum("swwt,pr", C["dd"], I3) \
      + 2/3*gp**2*my_einsum("prst", C["qd1"]) \
      - 8/3*gs**2*my_einsum("prst", C["qd8"]) \
      + 1/3*my_einsum("rs,pt", np.conj(Gd), Xid) \
      + 1/3*my_einsum("pt,rs", Gd, np.conj(Xid)) \
      + my_einsum("pr,st", Gu @ Gu.conj().T, C["phid"]) \
      - my_einsum("pr,st", Gd @ Gd.conj().T, C["phid"]) \
      + 2*my_einsum("st,pr", Gd.conj().T @ Gd, C["phiq1"]) \
      + 1/3*(my_einsum("pw,vs,vrwt", Gd, np.conj(Gd), C["qd1"]) \
      + 4/3*my_einsum("pw,vs,vrwt", Gd, np.conj(Gd), C["qd8"])) \
      + 1/3*(my_einsum("vt,rw,pvsw", Gd, np.conj(Gd), C["qd1"]) \
      + 4/3*my_einsum("vt,rw,pvsw", Gd, np.conj(Gd), C["qd8"])) \
      + 1/3*(my_einsum("rw,vs,vwpt", np.conj(Gu), np.conj(Gd), C["quqd1"]) \
      + 4/3*my_einsum("rw,vs,vwpt", np.conj(Gu), np.conj(Gd), C["quqd8"])) \
      + 1/3*(my_einsum("pw,vt,vwrs", Gu, Gd, np.conj(C["quqd1"])) \
      + 4/3*my_einsum("pw,vt,vwrs", Gu, Gd, np.conj(C["quqd8"]))) \
      + 1/2*my_einsum("ws,rv,pvwt", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      + 1/2*my_einsum("pv,wt,rvws", Gu, Gd, np.conj(C["quqd1"])) \
      - 2/3*(my_einsum("vt,ws,pvwr", Gd, np.conj(Gd), C["qq1"]) \
      + 3*my_einsum("vt,ws,pvwr", Gd, np.conj(Gd), C["qq3"])) \
      - 4*my_einsum("wt,vs,prvw", Gd, np.conj(Gd), C["qq1"]) \
      - 2/3*my_einsum("pv,rw,vtsw", Gd, np.conj(Gd), C["dd"]) \
      - 2*my_einsum("pv,rw,vwst", Gd, np.conj(Gd), C["dd"]) \
      - my_einsum("pv,rw,vwst", Gu, np.conj(Gu), C["ud1"]) \
      + my_einsum("pv,vrst", Gammaq, C["qd1"]) \
      + my_einsum("sv,prvt", Gammad, C["qd1"]) \
      + my_einsum("pvst,vr", C["qd1"], Gammaq) \
      + my_einsum("prsv,vt", C["qd1"], Gammad)

    Beta["qu8"] = 8/3*gs**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3)) \
      + 2/3*gs**2*my_einsum("prww,st", C["qu8"], I3) \
      + 2/3*gs**2*my_einsum("prww,st", C["qd8"], I3) \
      + 4/3*gs**2*my_einsum("wwst,pr", C["qu8"], I3) \
      + 2/3*gs**2*my_einsum("stww,pr", C["ud8"], I3) \
      + 8/3*gs**2*my_einsum("swwt,pr", C["uu"], I3) \
      - (4/3*gp**2 \
      + 14*gs**2)*my_einsum("prst", C["qu8"]) \
      - 12*gs**2*my_einsum("prst", C["qu1"]) \
      + 2*my_einsum("rs,pt", np.conj(Gu), Xiu) \
      + 2*my_einsum("pt,rs", Gu, np.conj(Xiu)) \
      + 2*(my_einsum("pw,vs,vrwt", Gu, np.conj(Gu), C["qu1"]) \
      - 1/6*my_einsum("pw,vs,vrwt", Gu, np.conj(Gu), C["qu8"])) \
      + 2*(my_einsum("vt,rw,pvsw", Gu, np.conj(Gu), C["qu1"]) \
      - 1/6*my_einsum("vt,rw,pvsw", Gu, np.conj(Gu), C["qu8"])) \
      + 2*(my_einsum("rw,vs,ptvw", np.conj(Gd), np.conj(Gu), C["quqd1"]) \
      - 1/6*my_einsum("rw,vs,ptvw", np.conj(Gd), np.conj(Gu), C["quqd8"])) \
      + 2*(my_einsum("pw,vt,rsvw", Gd, Gu, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("pw,vt,rsvw", Gd, Gu, np.conj(C["quqd8"]))) \
      + 1/2*my_einsum("vs,rw,vtpw", np.conj(Gu), np.conj(Gd), C["quqd8"]) \
      + 1/2*my_einsum("vt,pw,vsrw", Gu, Gd, np.conj(C["quqd8"])) \
      - 4*(my_einsum("vt,ws,pvwr", Gu, np.conj(Gu), C["qq1"]) \
      + 3*my_einsum("vt,ws,pvwr", Gu, np.conj(Gu), C["qq3"])) \
      - 4*my_einsum("pv,rw,vtsw", Gu, np.conj(Gu), C["uu"]) \
      - my_einsum("pv,rw,stvw", Gd, np.conj(Gd), C["ud8"]) \
      + my_einsum("pv,vrst", Gammaq, C["qu8"]) \
      + my_einsum("sv,prvt", Gammau, C["qu8"]) \
      + my_einsum("pvst,vr", C["qu8"], Gammaq) \
      + my_einsum("prsv,vt", C["qu8"], Gammau)

    Beta["qd8"] = 8/3*gs**2*(my_einsum("pwwr,st", C["qq1"], I3) \
      + 3*my_einsum("pwwr,st", C["qq3"], I3)) \
      + 2/3*gs**2*my_einsum("prww,st", C["qu8"], I3) \
      + 2/3*gs**2*my_einsum("prww,st", C["qd8"], I3) \
      + 4/3*gs**2*my_einsum("wwst,pr", C["qd8"], I3) \
      + 2/3*gs**2*my_einsum("wwst,pr", C["ud8"], I3) \
      + 8/3*gs**2*my_einsum("swwt,pr", C["dd"], I3) \
      - (-2/3*gp**2 \
      + 14*gs**2)*my_einsum("prst", C["qd8"]) \
      - 12*gs**2*my_einsum("prst", C["qd1"]) \
      + 2*my_einsum("rs,pt", np.conj(Gd), Xid) \
      + 2*my_einsum("pt,rs", Gd, np.conj(Xid)) \
      + 2*(my_einsum("pw,vs,vrwt", Gd, np.conj(Gd), C["qd1"]) \
      - 1/6*my_einsum("pw,vs,vrwt", Gd, np.conj(Gd), C["qd8"])) \
      + 2*(my_einsum("vt,rw,pvsw", Gd, np.conj(Gd), C["qd1"]) \
      - 1/6*my_einsum("vt,rw,pvsw", Gd, np.conj(Gd), C["qd8"])) \
      + 2*(my_einsum("rw,vs,vwpt", np.conj(Gu), np.conj(Gd), C["quqd1"]) \
      - 1/6*my_einsum("rw,vs,vwpt", np.conj(Gu), np.conj(Gd), C["quqd8"])) \
      + 2*(my_einsum("pw,vt,vwrs", Gu, Gd, np.conj(C["quqd1"])) \
      - 1/6*my_einsum("pw,vt,vwrs", Gu, Gd, np.conj(C["quqd8"]))) \
      + 1/2*my_einsum("vs,rw,pwvt", np.conj(Gd), np.conj(Gu), C["quqd8"]) \
      + 1/2*my_einsum("vt,pw,rwvs", Gd, Gu, np.conj(C["quqd8"])) \
      - 4*(my_einsum("vt,ws,pvwr", Gd, np.conj(Gd), C["qq1"]) \
      + 3*my_einsum("vt,ws,pvwr", Gd, np.conj(Gd), C["qq3"])) \
      - 4*my_einsum("pv,rw,vtsw", Gd, np.conj(Gd), C["dd"]) \
      - my_einsum("pv,rw,vwst", Gu, np.conj(Gu), C["ud8"]) \
      + my_einsum("pv,vrst", Gammaq, C["qd8"]) \
      + my_einsum("sv,prvt", Gammad, C["qd8"]) \
      + my_einsum("pvst,vr", C["qd8"], Gammaq) \
      + my_einsum("prsv,vt", C["qd8"], Gammad)

    Beta["ledq"] = -(8/3*gp**2 \
      + 8*gs**2)*my_einsum("prst", C["ledq"]) \
      - 2*my_einsum("ts,pr", np.conj(Gd), Xie) \
      - 2*my_einsum("pr,ts", Ge, np.conj(Xid)) \
      + 2*my_einsum("pv,tw,vrsw", Ge, np.conj(Gd), C["ed"]) \
      - 2*my_einsum("vr,tw,pvsw", Ge, np.conj(Gd), C["ld"]) \
      + 2*my_einsum("vr,ws,pvwt", Ge, np.conj(Gd), C["lq1"]) \
      + 6*my_einsum("vr,ws,pvwt", Ge, np.conj(Gd), C["lq3"]) \
      - 2*my_einsum("pw,vs,vtwr", Ge, np.conj(Gd), C["qe"]) \
      + 2*my_einsum("vs,tw,prvw", np.conj(Gd), np.conj(Gu), C["lequ1"]) \
      + my_einsum("pv,vrst", Gammal, C["ledq"]) \
      + my_einsum("sv,prvt", Gammad, C["ledq"]) \
      + my_einsum("pvst,vr", C["ledq"], Gammae) \
      + my_einsum("prsv,vt", C["ledq"], Gammaq)

    Beta["quqd1"] = 10/3*gp*my_einsum("st,pr", C["dB"], Gu) \
      - 6*g*my_einsum("st,pr", C["dW"], Gu) \
      - 20/9*gp*my_einsum("pt,sr", C["dB"], Gu) \
      + 4*g*my_einsum("pt,sr", C["dW"], Gu) \
      - 64/9*gs*my_einsum("pt,sr", C["dG"], Gu) \
      - 2/3*gp*my_einsum("pr,st", C["uB"], Gd) \
      - 6*g*my_einsum("pr,st", C["uW"], Gd) \
      + 4/9*gp*my_einsum("sr,pt", C["uB"], Gd) \
      + 4*g*my_einsum("sr,pt", C["uW"], Gd) \
      - 64/9*gs*my_einsum("sr,pt", C["uG"], Gd) \
      - 1/2*(11/9*gp**2 + 3*g**2 + 32*gs**2)*my_einsum("prst", C["quqd1"]) \
      - 1/3*( - 5/9*gp**2 - 3*g**2 + 64/3*gs**2)*my_einsum("srpt", C["quqd1"]) \
      - 4/9*( - 5/9*gp**2 - 3*g**2 + 28/3*gs**2)*my_einsum("srpt", C["quqd8"]) \
      + 16/9*gs**2*my_einsum("prst", C["quqd8"]) \
      - 2*my_einsum("pr,st", Gu, Xid) \
      - 2*my_einsum("st,pr", Gd, Xiu) \
      + 4/3*(my_einsum("vr,pw,svwt", Gu, Gd, C["qd1"]) \
      + 4/3*my_einsum("vr,pw,svwt", Gu, Gd, C["qd8"]) \
      + my_einsum("vt,sw,pvwr", Gd, Gu, C["qu1"]) \
      + 4/3*my_einsum("vt,sw,pvwr", Gd, Gu, C["qu8"]) \
      + my_einsum("pw,sv,vrwt", Gd, Gu, C["ud1"]) \
      + 4/3*my_einsum("pw,sv,vrwt", Gd, Gu, C["ud8"])) \
      + 8/3*(my_einsum("wt,vr,svpw", Gd, Gu, C["qq1"]) \
      - 3*my_einsum("wt,vr,svpw", Gd, Gu, C["qq3"]) \
      - 3*my_einsum("wt,vr,swpv", Gd, Gu, C["qq1"]) \
      + 9*my_einsum("wt,vr,swpv", Gd, Gu, C["qq3"])) \
      - 4*my_einsum("sw,pv,vrwt", Gd, Gu, C["ud1"]) \
      + my_einsum("pv,vrst", Gammaq, C["quqd1"]) \
      + my_einsum("sv,prvt", Gammaq, C["quqd1"]) \
      + my_einsum("pvst,vr", C["quqd1"], Gammau) \
      + my_einsum("prsv,vt", C["quqd1"], Gammad)

    Beta["quqd8"] = 8*gs*my_einsum("st,pr", C["dG"], Gu) \
      - 40/3*gp*my_einsum("pt,sr", C["dB"], Gu) \
      + 24*g*my_einsum("pt,sr", C["dW"], Gu) \
      + 16/3*gs*my_einsum("pt,sr", C["dG"], Gu) \
      + 8*gs*my_einsum("pr,st", C["uG"], Gd) \
      + 8/3*gp*my_einsum("sr,pt", C["uB"], Gd) \
      + 24*g*my_einsum("sr,pt", C["uW"], Gd) \
      + 16/3*gs*my_einsum("sr,pt", C["uG"], Gd) \
      + 8*gs**2*my_einsum("prst", C["quqd1"]) \
      + (10/9*gp**2 + 6*g**2 + 16/3*gs**2)*my_einsum("srpt", C["quqd1"]) \
      + (-11/18*gp**2 - 3/2*g**2 + 16/3*gs**2)*my_einsum("prst", C["quqd8"]) \
      - 1/3*(5/9*gp**2 + 3*g**2 \
      + 44/3*gs**2)*my_einsum("srpt", C["quqd8"]) \
      + 8*(my_einsum("vr,pw,svwt", Gu, Gd, C["qd1"]) \
      - 1/6*my_einsum("vr,pw,svwt", Gu, Gd, C["qd8"]) \
      + my_einsum("vt,sw,pvwr", Gd, Gu, C["qu1"]) \
      - 1/6*my_einsum("vt,sw,pvwr", Gd, Gu, C["qu8"]) \
      + my_einsum("pw,sv,vrwt", Gd, Gu, C["ud1"]) \
      - 1/6*my_einsum("pw,sv,vrwt", Gd, Gu, C["ud8"])) \
      + 16*(my_einsum("wt,vr,svpw", Gd, Gu, C["qq1"]) \
      - 3*my_einsum("wt,vr,svpw", Gd, Gu, C["qq3"])) \
      - 4*my_einsum("sw,pv,vrwt", Gd, Gu, C["ud8"]) \
      + my_einsum("pv,vrst", Gammaq, C["quqd8"]) \
      + my_einsum("sv,prvt", Gammaq, C["quqd8"]) \
      + my_einsum("pvst,vr", C["quqd8"], Gammau) \
      + my_einsum("prsv,vt", C["quqd8"], Gammad)

    Beta["lequ1"] = -(11/3*gp**2 + 8*gs**2)*my_einsum("prst", C["lequ1"]) \
      + (30*gp**2 + 18*g**2)*my_einsum("prst", C["lequ3"]) \
      + 2*my_einsum("st,pr", Gu, Xie) \
      + 2*my_einsum("pr,st", Ge, Xiu) \
      + 2*my_einsum("sv,wt,prvw", Gd, Gu, C["ledq"]) \
      + 2*my_einsum("pv,sw,vrwt", Ge, Gu, C["eu"]) \
      + 2*my_einsum("vr,wt,pvsw", Ge, Gu, C["lq1"]) \
      - 6*my_einsum("vr,wt,pvsw", Ge, Gu, C["lq3"]) \
      - 2*my_einsum("vr,sw,pvwt", Ge, Gu, C["lu"]) \
      - 2*my_einsum("pw,vt,svwr", Ge, Gu, C["qe"]) \
      + my_einsum("pv,vrst", Gammal, C["lequ1"]) \
      + my_einsum("sv,prvt", Gammaq, C["lequ1"]) \
      + my_einsum("pvst,vr", C["lequ1"], Gammae) \
      + my_einsum("prsv,vt", C["lequ1"], Gammau)

    Beta["lequ3"] = 5/6*gp*my_einsum("pr,st", C["eB"], Gu) \
      - 3/2*g*my_einsum("st,pr", C["uW"], Ge) \
      - 3/2*gp*my_einsum("st,pr", C["uB"], Ge) \
      - 3/2*g*my_einsum("pr,st", C["eW"], Gu) \
      + (2/9*gp**2 - 3*g**2 + 8/3*gs**2)*my_einsum("prst", C["lequ3"]) \
      + 1/8*(5*gp**2 + 3*g**2)*my_einsum("prst", C["lequ1"]) \
      - 1/2*my_einsum("sw,pv,vrwt", Gu, Ge, C["eu"]) \
      - 1/2*my_einsum("vr,wt,pvsw", Ge, Gu, C["lq1"]) \
      + 3/2*my_einsum("vr,wt,pvsw", Ge, Gu, C["lq3"]) \
      - 1/2*my_einsum("vr,sw,pvwt", Ge, Gu, C["lu"]) \
      - 1/2*my_einsum("pw,vt,svwr", Ge, Gu, C["qe"]) \
      + my_einsum("pv,vrst", Gammal, C["lequ3"]) \
      + my_einsum("sv,prvt", Gammaq, C["lequ3"]) \
      + my_einsum("pvst,vr", C["lequ3"], Gammae) \
      + my_einsum("prsv,vt", C["lequ3"], Gammau)

    Beta["duql"] = -(9/2*g**2 \
      + 11/6*gp**2 \
      + 4*gs**2)*my_einsum("prst", C["duql"]) \
      - my_einsum("sv,wp,vrwt", np.conj(Gd), Gd, C["duql"]) \
      - my_einsum("sv,wr,pvwt", np.conj(Gu), Gu, C["duql"]) \
      + 2*my_einsum("tv,sw,prwv", np.conj(Ge), np.conj(Gu), C["duue"]) \
      + my_einsum("tv,sw,pwrv", np.conj(Ge), np.conj(Gu), C["duue"]) \
      + 4*my_einsum("vp,wr,vwst", Gd, Gu, C["qqql"]) \
      + 4*my_einsum("vp,wr,wvst", Gd, Gu, C["qqql"]) \
      - my_einsum("vp,wr,vswt", Gd, Gu, C["qqql"]) \
      - my_einsum("vp,wr,wsvt", Gd, Gu, C["qqql"]) \
      + 2*my_einsum("wp,tv,wsrv", Gd, np.conj(Ge), C["qque"]) \
      + my_einsum("vp,vrst", Gd.conj().T @ Gd, C["duql"]) \
      + my_einsum("vr,pvst", Gu.conj().T @ Gu, C["duql"]) \
      + 1/2*(my_einsum("vs,prvt", Gu @ Gu.conj().T, C["duql"]) \
      + my_einsum("vs,prvt", Gd @ Gd.conj().T, C["duql"])) \
      + 1/2*my_einsum("vt,prsv", Ge @ Ge.conj().T, C["duql"])

    Beta["qque"] = -(9/2*g**2 \
      + 23/6*gp**2 + 4*gs**2)*my_einsum("prst", C["qque"]) \
      - my_einsum("rv,ws,pwvt", np.conj(Gu), Gu, C["qque"]) \
      + 1/2*my_einsum("wt,rv,vspw", Ge, np.conj(Gd), C["duql"]) \
      - 1/2*(2*my_einsum("pv,rw,vwst", np.conj(Gd), np.conj(Gu), C["duue"]) \
      + my_einsum("pv,rw,vswt", np.conj(Gd), np.conj(Gu), C["duue"])) \
      + 1/2*( \
      - 2*my_einsum("ws,vt,prwv", Gu, Ge, C["qqql"]) \
      + my_einsum("ws,vt,pwrv", Gu, Ge, C["qqql"]) \
      - 2*my_einsum("ws,vt,wprv", Gu, Ge, C["qqql"])) \
      + 1/2*(my_einsum("vp,vrst", Gu @ Gu.conj().T, C["qque"]) \
      + my_einsum("vp,vrst", Gd @ Gd.conj().T, C["qque"])) \
      - my_einsum("pv,ws,rwvt", np.conj(Gu), Gu, C["qque"]) \
      + 1/2*my_einsum("wt,pv,vsrw", Ge, np.conj(Gd), C["duql"]) \
      - 1/2*(2*my_einsum("rv,pw,vwst", np.conj(Gd), np.conj(Gu), C["duue"]) \
      + my_einsum("rv,pw,vswt", np.conj(Gd), np.conj(Gu), C["duue"])) \
      + 1/2*( \
      - 2*my_einsum("ws,vt,rpwv", Gu, Ge, C["qqql"]) \
      + my_einsum("ws,vt,rwpv", Gu, Ge, C["qqql"]) \
      - 2*my_einsum("ws,vt,wrpv", Gu, Ge, C["qqql"])) \
      + 1/2*(my_einsum("vr,vpst", Gu @ Gu.conj().T, C["qque"]) \
      + my_einsum("vr,vpst", Gd @ Gd.conj().T, C["qque"])) \
      + my_einsum("vs,prvt", Gu.conj().T @ Gu, C["qque"]) \
      + my_einsum("vt,prsv", Ge.conj().T @ Ge, C["qque"])

    Beta["qqql"] = -(3*g**2 \
      + 1/3*gp**2 + 4*gs**2)*my_einsum("prst", C["qqql"]) \
      - 4*g**2*(my_einsum("rpst", C["qqql"]) \
      + my_einsum("srpt", C["qqql"]) \
      + my_einsum("psrt", C["qqql"])) \
      - 4*my_einsum("tv,sw,prwv", np.conj(Ge), np.conj(Gu), C["qque"]) \
      + 2*(my_einsum("pv,rw,vwst", np.conj(Gd), np.conj(Gu), C["duql"]) \
      + my_einsum("rv,pw,vwst", np.conj(Gd), np.conj(Gu), C["duql"])) \
      + 1/2*(my_einsum("vp,vrst", Gu @ Gu.conj().T, C["qqql"]) \
      + my_einsum("vp,vrst", Gd @ Gd.conj().T, C["qqql"])) \
      + 1/2*(my_einsum("vr,pvst", Gu @ Gu.conj().T, C["qqql"]) \
      + my_einsum("vr,pvst", Gd @ Gd.conj().T, C["qqql"])) \
      + 1/2*(my_einsum("vs,prvt", Gu @ Gu.conj().T, C["qqql"]) \
      + my_einsum("vs,prvt", Gd @ Gd.conj().T, C["qqql"])) \
      + 1/2*my_einsum("vt,prsv", Ge @ Ge.conj().T, C["qqql"])

    Beta["duue"] = -(2*gp**2 + 4*gs**2)*my_einsum("prst", C["duue"]) \
      - 20/3*gp**2*my_einsum("psrt", C["duue"]) \
      + 4*my_einsum("ws,vt,prwv", Gu, Ge, C["duql"]) \
      - 8*my_einsum("vp,wr,vwst", Gd, Gu, C["qque"]) \
      + my_einsum("vp,vrst", Gd.conj().T @ Gd, C["duue"]) \
      + my_einsum("vr,pvst", Gu.conj().T @ Gu, C["duue"]) \
      + my_einsum("vs,prvt", Gu.conj().T @ Gu, C["duue"]) \
      + my_einsum("vt,prsv", Ge.conj().T @ Ge, C["duue"])

    Beta["llphiphi"] = (2*Lambda \
      - 3*g**2 \
      + 2*GammaH)*C["llphiphi"]-3/2*(C["llphiphi"] @ Ge @ Ge.conj().T \
      + Ge.conj() @ Ge.T @ C["llphiphi"])

    return Beta
852224
def rho2theta(self, rho0, gamma):
        
        fac = np.sqrt(np.pi) * special.gamma(1. / 2 * (-1 + gamma)) / special.gamma(gamma / 2.) * 2 / (3 - gamma) * rho0

        #fac = theta_E**(gamma - 1)
        theta_E = fac**(1. / (gamma - 1))
        return theta_E
1113853
def _c0(self):
        ""
        h_df = self.df / 2
        p, S = self._p, self.S

        return h_df * (logdet(S) + p * logtwo) + lpgamma(p, h_df)
1206040
def G12(x, a):
    
    alpha, a, beta, b = a
    pi26 = np.pi ** 2 / 6.0
    G = (pi26 + x) * np.exp(-x)
    tmp = 1 + b * x ** beta
    g = 1.0 / (a * x ** alpha / tmp + 1.0)
    return G * g
631837
def hugoniot_p(rho, rho0, c0, s):
    
    eta = 1. - (rho0 / rho)
    Ph = rho0 * c0 * c0 * eta / np.power((1. - s * eta), 2.)
    return Ph
367475
def _water(cls, T, P):
        """"""
        water = IAPWS95(P=P, T=T)
        prop = {}
        prop["g"] = water.h-T*water.s
        prop["gt"] = -water.s
        prop["gp"] = 1./water.rho
        prop["gtt"] = -water.cp/T
        prop["gtp"] = water.betas*water.cp/T
        prop["gpp"] = -1e6/(water.rho*water.w)**2-water.betas**2*1e3*water.cp/T
        prop["gs"] = 0
        prop["gsp"] = 0
        prop["thcond"] = water.k
        return prop
1327807
private double[][] assignPiFactors(IAtomContainerSet setAc) {
        //a,b,c,denom,chi,q
        double[][] gasteigerFactors = new double[setAc.getAtomContainerCount()][(setAc.getAtomContainer(0)
                .getAtomCount() * (STEP_SIZE + 1))];
        String AtomSymbol = "";
        double[] factors = new double[]{0.0, 0.0, 0.0};
        for (int k = 1; k < setAc.getAtomContainerCount(); k++) {
            IAtomContainer ac = setAc.getAtomContainer(k);
            for (int i = 0; i < ac.getAtomCount(); i++) {
                factors[0] = 0.0;
                factors[1] = 0.0;
                factors[2] = 0.0;
                AtomSymbol = ac.getAtom(i).getSymbol();
                if (AtomSymbol.equals("H")) {
                    factors[0] = 0.0;
                    factors[1] = 0.0;
                    factors[2] = 0.0;
                } else if (AtomSymbol.equals("C")) {/*
                                                     * if(ac.getAtom(i).getFlag(
                                                     * ISCHANGEDFC))
                                                     */
                    {
                        factors[0] = 5.60;
                        factors[1] = 8.93;
                        factors[2] = 2.94;
                    }
                } else if (AtomSymbol.equals("O")) {
                    if (ac.getMaximumBondOrder(ac.getAtom(i)) == IBond.Order.SINGLE) {
                        factors[0] = 10.0;
                        factors[1] = 13.86;
                        factors[2] = 9.68;
                    } else {
                        factors[0] = 7.91;
                        factors[1] = 14.76;
                        factors[2] = 6.85;
                    }
                } else if (AtomSymbol.equals("N")) {
                    if (ac.getMaximumBondOrder(ac.getAtom(i)) != IBond.Order.SINGLE) {
                        factors[0] = 7.95;/* 7.95 */
                        factors[1] = 9.73;/* 9.73 */
                        factors[2] = 2.67;/* 2.67 */
                    } else {
                        factors[0] = 4.54;/* 4.54 *//* 5.5 */
                        factors[1] = 11.86;/* 11.86 *//* 10.86 */
                        factors[2] = 7.32;/* 7.32 *//* 7.99 */
                    }
                } else if (AtomSymbol.equals("S")) {
                    if (ac.getMaximumBondOrder(ac.getAtom(i)) == IBond.Order.SINGLE) {
                        factors[0] = 7.73;
                        factors[1] = 8.16;
                        factors[2] = 1.81;
                    } else {
                        factors[0] = 6.60;
                        factors[1] = 10.32;
                        factors[2] = 3.72;
                    }
                } else if (AtomSymbol.equals("F")) {
                    factors[0] = 7.34;
                    factors[1] = 13.86;
                    factors[2] = 9.68;
                } else if (AtomSymbol.equals("Cl")) {
                    factors[0] = 6.50;
                    factors[1] = 11.02;
                    factors[2] = 4.52;
                } else if (AtomSymbol.equals("Br")) {
                    factors[0] = 5.20;
                    factors[1] = 9.68;
                    factors[2] = 4.48;
                } else if (AtomSymbol.equals("I")) {
                    factors[0] = 4.95;
                    factors[1] = 8.81;
                    factors[2] = 3.86;
                }

                gasteigerFactors[k][STEP_SIZE * i + i] = factors[0];
                gasteigerFactors[k][STEP_SIZE * i + i + 1] = factors[1];
                gasteigerFactors[k][STEP_SIZE * i + i + 2] = factors[2];
                gasteigerFactors[k][STEP_SIZE * i + i + 5] = ac.getAtom(i).getCharge();

                if (factors[0] == 0 && factors[1] == 0 && factors[2] == 0) {
                    gasteigerFactors[k][STEP_SIZE * i + i + 3] = 1;
                } else {
                    gasteigerFactors[k][STEP_SIZE * i + i + 3] = factors[0] + factors[1] + factors[2];
                }
            }
        }

        return gasteigerFactors;
    }
1183589
def sigma_r2(self, r, a, gamma, rho0_r0_gamma, r_ani):
        
        # first term
        prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
        prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
        hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
        hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
        fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
        return prefac1 * prefac2 * fac * (self._cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000) ** 2
442663
def _g(self, h, xp, s):
        """"""

        nphi = sum(self.phi)

        return (nphi / 2.0) * log(2 * pi) + nphi * \
            log(s) + 0.5 * sum((h - xp) ** 2) / (s ** 2)
394547
private static function hsl2rgb($h, $s, $l, &$r, &$g, &$b)
    {
        if ($s == 0) {
            $r = round($l * 255);
            $g = $r;
            $b = $r;

            return;
        }

        $temp1;
        if ($l < 0.5) {
            $temp1 = $l * (1.0 + $s);
        } else {
            $temp1 = $l + $s - $l * $s;
        }

        $temp2 = 2 * $l - $temp1;

        $h = $h / 360;

        $tempR = $h + (1 / 3);
        $tempG = $h;
        $tempB = $h - (1 / 3);

        if ($tempR < 0) {
            $tempR++;
        } elseif ($tempR > 1) {
            $tempR--;
        }

        if ($tempG < 0) {
            $tempG++;
        } elseif ($tempG > 1) {
            $tempG--;
        }

        if ($tempB < 0) {
            $tempB++;
        } elseif ($tempB > 1) {
            $tempB--;
        }

        if (6 * $tempR < 1) {
            $r = $temp2 + ($temp1 - $temp2) * 6 * $tempR;
        } elseif (2 * $tempR < 1) {
            $r = $temp1;
        } elseif (3 * $tempR < 2) {
            $r = $temp2 + ($temp1 - $temp2) * ((2 / 3) - $tempR) * 6;
        } else {
            $r = $temp2;
        }

        if (6 * $tempG < 1) {
            $g = $temp2 + ($temp1 - $temp2) * 6 * $tempG;
        } elseif (2 * $tempG < 1) {
            $g = $temp1;
        } elseif (3 * $tempG < 2) {
            $g = $temp2 + ($temp1 - $temp2) * ((2 / 3) - $tempG) * 6;
        } else {
            $g = $temp2;
        }

        if (6 * $tempB < 1) {
            $b = $temp2 + ($temp1 - $temp2) * 6 * $tempB;
        } elseif (2 * $tempB < 1) {
            $b = $temp1;
        } elseif (3 * $tempB < 2) {
            $b = $temp2 + ($temp1 - $temp2) * ((2 / 3) - $tempB) * 6;
        } else {
            $b = $temp2;
        }

        $r = round($r * 255);
        $g = round($g * 255);
        $b = round($b * 255);
    }
993727
private static double gSer(double x, double A) {
        // Good for X<A+1.
        double T9=1/A;
        double G=T9;
        double I=1;
        while (T9>G*0.00001) {
            T9=T9*x/(A+I);
            G=G+T9;
            ++I;
        }
        G=G*Math.exp(A*Math.log(x)-x-logGamma(A));

        return G;
    }
1092879
def wavelengthToRGB(wavelength):
	gamma = 0.80;
	intensityMax = 255;

	

	factor = None
	r = None
	g = None
	b = None

	if((wavelength >= 380) and (wavelength<440)):
		r = -(wavelength - 440) / (440.0 - 380.0)
		g = 0.0
		b = 1.0
	elif((wavelength >= 440) and (wavelength<490)):
		r = 0.0
		g = (wavelength - 440) / (490.0 - 440.0)
		b = 1.0
	elif((wavelength >= 490) and (wavelength<510)):
		r = 0.0
		g = 1.0
		b = -(wavelength - 510) / (510.0 - 490.0)
	elif((wavelength >= 510) and (wavelength<580)):
		r = (wavelength - 510) / (580.0 - 510.0)
		g = 1.0
		b = 0.0
	elif((wavelength >= 580) and (wavelength<645)):
		r = 1.0
		g = -(wavelength - 645) / (645.0 - 580.0)
		b = 0.0
	elif((wavelength >= 645) and (wavelength<781)):
		r = 1.0
		g = 0.0
		b = 0.0
	else:
		r = 0.0
		g = 0.0
		b = 0.0

	# Let the intensity fall off near the vision limits
	if((wavelength >= 380) and (wavelength<420)):
	    factor = 0.3 + 0.7*(wavelength - 380) / (420.0 - 380.0)
	elif((wavelength >= 420) and (wavelength<701)):
	    factor = 1.0
	elif((wavelength >= 701) and (wavelength<781)):
	    factor = 0.3 + 0.7*(780 - wavelength) / (780.0 - 700.0)
	else:
	    factor = 0.0

	# Don't want 0^x = 1 for x != 0

	r = int(round(intensityMax * pow(r * factor, gamma)))
	g = int(round(intensityMax * pow(g * factor, gamma)))
	b = int(round(intensityMax * pow(b * factor, gamma)))

	return (r, g, b)