(*----------------------------------------------------------------------*)
(*----------------------------------------------------------------------

This mathematica package contains the two-loop decoupling coefficients for alpha_1/2/3
for matching the Georgi-Glashow SU(5) model to the Standard Model.
The main functions are:
zetaAlU1, izetaAlU1, zetaAlSU2, izetaAlSU2, zetaAlSU3, izetaAlSU3.

The arguments are:
         alpha4pi = alpha(mu)/(4*pi) the unique gauge coupling in the GUT in MSbar at the scale mu
         alpha14pi = alpha_1(mu)/(4*pi) the U(1) gauge coupling in the SM in MSbar at the scale mu
         alpha24pi = alpha_2(mu)/(4*pi) the SU(2) gauge coupling in the SM in MSbar at the scale mu
         alpha34pi = alpha_3(mu)/(4*pi) the SU(3) gauge coupling in the SM in MSbar at the scale mu
         yt = the top Yukawa coupling in the GUT (for zetaAl...) or the SM (for izetaAl...) in MSbar
         yb = the (unified) bottom-tau Yukawa coupling in the GUT (for zetaAl...) or the SM (for izetaAl...) in MSbar
         Mx = the on-shell mass of the heavy gauge boson living in (3,2bar,-5/6) + (3bar,2,5/6) of 24
         MHc = the on-shell mass of the colored Higgs triplet (3,1,-1/3) in 5_H 
         Msig = the on-shell mass of the color oktet Higgs (8,1,0) in 24_H
         M24 = the on-shell mass of the singlet Higgs (1,1,0) in 24_H
         mu = the decoupling scale
         Loops = the number of loops of the RGE analysis, 
                 i.e. Loops=3 gives the two-loop decoupling coefficient,
                 Loops=2 gives the one-loop decoupling coefficient

The description of their arguments can also be obtained by e.g. the command
?zetaAlU1

Important Note: We perform a one-step decoupling here. Hence, be careful when choosing
                the mass spectrum. Too far separated masses can lead to power-enhanced
                contributions and spoil the perturbative expansion (check whether
                the two-loop contribution is smaller than the one-loop contribution).
                Usually up to two orders of magnitude separation is still safe.



For more details, please have a look at the paper:

  "Towards a Two-Loop Matching of Gauge Couplings in Grand Unified Theories "

                                 W. Martens

Or contact me:
martens@particle.uni-karlsruhe.de

  ----------------------------------------------------------------------*)
(*----------------------------------------------------------------------*)



BeginPackage["zetaalgut`"];

usagestring[zeta_,alpha_,group_,inverse_,i_,GUTSM_]:= zeta<>"["<>alpha<>"4pi, yt, yb, Mx, MHc, Msig, M24, mu, Loops] gives the 
"<>inverse<>" (up to two-loop) decoupling coefficient for alpha"<>i<>" 
for matching the Georgi-Glashow SU(5) to the Standard Model.
The arguments are:
         "<>alpha<>"4pi = "<>alpha<>"(mu)/(4*pi) the "<>group<>" gauge coupling in the "<>GUTSM<>" in MSbar at the scale mu
         yt = the top Yukawa coupling in the "<>GUTSM<>" in MSbar
         yb = the (unified) bottom-tau Yukawa coupling in the "<>GUTSM<>" in MSbar
         Mx = the on-shell mass of the heavy gauge boson living in (3,2bar,-5/6) + (3bar,2,5/6) of 24
         MHc = the on-shell mass of the colored Higgs triplet (3,1,-1/3) in 5_H 
         Msig = the on-shell mass of the color oktet Higgs (8,1,0) in 24_H
         M24 = the on-shell mass of the singlet Higgs (1,1,0) in 24_H
         mu = the decoupling scale
         Loops = the number of loops of the RGE analysis, 
         i.e. Loops=3 gives the two-loop decoupling coefficient,
         Loops=2 gives the one-loop decoupling coefficient ";


zetaAlU1::usage=usagestring["zetaAlU1","alpha","unique","","1","GUT"]
izetaAlU1::usage=usagestring["izetaAlU1","alpha1","U(1)","inverse","1","SM"]
zetaAlSU2::usage=usagestring["zetaAlSU2","alpha","unique","","2","GUT"]
izetaAlSU2::usage=usagestring["izetaAlSU2","alpha2","SU(2)","inverse","2","SM"]
zetaAlSU3::usage=usagestring["zetaAlSU3","alpha","unique","","3","GUT"]
izetaAlSU3::usage=usagestring["izetaAlSU3","alpha3","SU(3)","inverse","3","SM"]



Begin["`Private`"];



(*----------------------------------------------------------------------*)
(*--------------- decoupling coefficient for alpha1 --------------------*)
(*----------------------------------------------------------------------*)

zetaAlU1[alpha4pi_, yt_, yb_, Mx_, MHc_, Msig_, M24_, mu_, Loops_] :=  
    1 + alpha4pi^2*((1545962 + 375*Msig*(16/(Msig - 2*Mx) + 9/(Msig - Mx)) + 
         (1875*M24)/(M24 - 2*Mx) + (1750*(5*M24^4 + 6*MHc^4 + 160*Msig^4))/
          Mx^4 + (29160*M24^2*MHc^6*(M24^2 - 4*MHc^2) + 
           2*MHc^2*(-57500*M24^6 + 173391*M24^4*MHc^2 + 239124*M24^2*MHc^4 - 
             46656*MHc^6)*Msig^2 + 2*(14375*M24^6 - 643539*M24^4*MHc^2 + 
             2340984*M24^2*MHc^4 + 11664*MHc^6)*Msig^4 + 
           299000*M24^2*(M24^2 - 4*MHc^2)*Msig^6)/(M24^2*(M24^2 - 4*MHc^2)*
           Msig^2*(4*MHc^2 - Msig^2)*Mx^2) - (240*Mx^2)/MHc^2 + 
         (3375*Msig)/(Msig + Mx) + (1875*M24)/(M24 + 2*Mx) + 
         (6000*Msig)/(Msig + 2*Mx))/3600 + 
       (32*MHc^2*B0fin[MHc^2, MHc, M24, mu])/(225*Mx^2) + 
       (128*MHc^2*B0fin[MHc^2, MHc, Msig, mu])/(45*Mx^2) + 
       ((MHc^2 - Mx^2)^2*B0fin[MHc^2, Mx, 0, mu])/(15*MHc^2*Mx^2) - 
       (35*(MHc^2 - Mx^2)^2*B0fin[Mx^2, MHc, 0, mu])/(12*Mx^4) - 
       (175*(M24^4 - 4*M24^2*Mx^2 + 12*Mx^4)*B0fin[Mx^2, Mx, M24, mu])/
        (72*Mx^4) - (70*(Msig^4 - 4*Msig^2*Mx^2 + 12*Mx^4)*
         B0fin[Mx^2, Mx, Msig, mu])/(9*Mx^4) + 
       (35*(-3 - (4*Msig^4)/Mx^4 + (4*Msig^2)/Mx^2)*B0fin[Mx^2, Mx, 2*Msig, 
          mu])/2 + (64*MHc^8*CapitalPhiLambdaM[M24, MHc, MHc])/
        (75*M24^2*(M24^2 - 4*MHc^2)^3*Mx^2) + 
       (25*Mx^2*(3*M24^4 - 12*M24^2*Mx^2 + 4*Mx^4)*CapitalPhiLambdaM[M24, Mx, 
          Mx])/(4*M24^2*(M24^2 - 4*Mx^2)^3) + 
       (256*MHc^8*CapitalPhiLambdaM[Msig, MHc, MHc])/
        (15*Msig^2*(-4*MHc^2 + Msig^2)^3*Mx^2) + 
       (20*Mx^2*(3*Msig^4 - 12*Msig^2*Mx^2 + 4*Mx^4)*CapitalPhiLambdaM[Msig, 
          Mx, Mx])/(Msig^2*(Msig^2 - 4*Mx^2)^3) + 
       (45*Mx^2*(12*Msig^4 - 12*Msig^2*Mx^2 + Mx^4)*CapitalPhiLambdaM[2*Msig, 
          Mx, Mx])/(256*Msig^2*(Msig^2 - Mx^2)^3) + 
       ((-32*MHc^2*(M24^4 - 7*M24^2*MHc^2 + 6*MHc^4))/
          (225*(M24^2 - 4*MHc^2)^2*Mx^2) - 
         (25*(-7*M24^8 + 77*M24^6*Mx^2 - 322*M24^4*Mx^4 + 468*M24^2*Mx^6 + 
            72*Mx^8))/(72*Mx^4*(M24^2 - 4*Mx^2)^2))*Log[mu^2/M24^2] + 
       Log[mu^2/MHc^2]^2/225 - (5*(-224*Msig^8 + 616*Msig^6*Mx^2 - 
          644*Msig^4*Mx^4 + 234*Msig^2*Mx^6 + 9*Mx^8)*Log[mu^2/(4*Msig^2)])/
        (16*Mx^4*(Msig^2 - Mx^2)^2) + 
       (2*((-64*MHc^2*(6*MHc^4 - 7*MHc^2*Msig^2 + Msig^4)*Mx^2)/
           (-4*MHc^2 + Msig^2)^2 + (25*(7*Msig^8 - 77*Msig^6*Mx^2 + 
             322*Msig^4*Mx^4 - 468*Msig^2*Mx^6 - 72*Mx^8))/
           (Msig^2 - 4*Mx^2)^2)*Log[mu^2/Msig^2])/(45*Mx^4) + 
       Log[mu^2/MHc^2]*((-2416 + (5250*MHc^4)/Mx^4 + 
           (4*(-80*MHc^6*(11184*MHc^4 - 4568*MHc^2*Msig^2 + 539*Msig^4) - 
              5*M24^4*(11056*MHc^6 - 4504*MHc^4*Msig^2 + 531*MHc^2*Msig^4) + 
              8*M24^2*(55408*MHc^8 - 22584*MHc^6*Msig^2 + 2663*MHc^4*
                 Msig^4)))/((M24^2 - 4*MHc^2)^2*(-4*MHc^2 + Msig^2)^2*Mx^2) + 
           (2835*MHc)/(-MHc + Mx) - (2835*MHc)/(MHc + Mx))/1800 - 
         (7*Log[mu^2/Mx^2])/3) + ((982920 - (192*Mx^2)/MHc^2 + 
          (9072*MHc^2)/(MHc^2 - Mx^2) + (500*M24^2*(-14*M24^4 + 
             19*M24^2*Mx^2 + 112*Mx^4))/(Mx^2*(M24^2 - 4*Mx^2)^2) - 
          (100*Msig^2*(728*Msig^8 - 4955*Msig^6*Mx^2 + 8220*Msig^4*Mx^4 + 
             2560*Msig^2*Mx^6 - 5824*Mx^8))/(Mx^2*(Msig^4 - 5*Msig^2*Mx^2 + 
              4*Mx^4)^2))*Log[mu^2/Mx^2])/2880 + (1225*Log[mu^2/Mx^2]^2)/4)*
      UnitStep[-3 + Loops] + alpha4pi*
      (yt^2*(1/(16*Pi^2) + (7*Log[mu^2/MHc^2])/(40*Pi^2))*
        UnitStep[-3 + Loops] + yb^2*(67/(320*Pi^2) + (3*Log[mu^2/MHc^2])/
          (16*Pi^2))*UnitStep[-3 + Loops] + (5/3 - Log[mu^2/MHc^2]/15 + 
		      (35*Log[mu^2/Mx^2])/2)*UnitStep[-2 + Loops]);


(* inverse decoupling coefficient for alpha1 *)
izetaAlU1[alpha14pi_, yt_, yb_, Mx_, MHc_, Msig_, M24_, mu_, Loops_] := Module[{A,B,izeta,alpha4pi},

     A = zetaAlU1[alpha14pi,yt,yb,Mx,MHc,Msig,M24,mu,2] - 1;
     alpha4pi = (1-UnitStep[-3 + Loops]*A)*alpha14pi;

     A = zetaAlU1[alpha4pi,yt,yb,Mx,MHc,Msig,M24,mu,2] - 1;
     B = zetaAlU1[alpha4pi,yt,yb,Mx,MHc,Msig,M24,mu,3] - A - 1;


     izeta = UnitStep[-3 + Loops]*(A^2 - B) - UnitStep[-2 + Loops]*A + 1;
     izeta
];




(*----------------------------------------------------------------------*)
(*--------------- decoupling coefficient for alpha2 --------------------*)
(*----------------------------------------------------------------------*)

zetaAlSU2[alpha4pi_, yt_, yb_, Mx_, MHc_, Msig_, M24_, mu_, Loops_] := Re[
    1 + alpha4pi^2*((25*M24^10*(Msig - 2*Mx)*(Msig - Mx)*(Msig + Mx)*
          (Msig + 2*Mx)*(42*Msig^2 - Mx^2) - 50*M24^8*(Msig - 2*Mx)*
          (Msig - Mx)*(Msig + Mx)*(Msig + 2*Mx)*(336*Msig^4 + 
           151*Msig^2*Mx^2 - 2*Mx^4) + 62208*MHc^4*Msig^4*Mx^4*
          (Msig^4 - 5*Msig^2*Mx^2 + 4*Mx^4) + 16*M24^2*Msig^2*Mx^2*
          (9*MHc^4*(Msig - 2*Mx)*(Msig - Mx)*(Msig + Mx)*(Msig + 2*Mx)*
            (452*Msig^2 - 47*Mx^2) - 5400*MHc^2*Msig^2*Mx^2*
            (Msig^4 - 5*Msig^2*Mx^2 + 4*Mx^4) + 40*(3360*Msig^10 - 
             20876*Msig^8*Mx^2 + 48588*Msig^6*Mx^4 - 89703*Msig^4*Mx^6 + 
             59198*Msig^2*Mx^8 - 648*Mx^10)) + 
         5*M24^6*(9*MHc^4*(Msig - 2*Mx)*(Msig - Mx)*(Msig + Mx)*(Msig + 2*Mx)*
            (28*Msig^2 - Mx^2) - 270*MHc^2*Msig^2*Mx^2*(Msig^4 - 
             5*Msig^2*Mx^2 + 4*Mx^4) + 2*(3360*Msig^10 - 7796*Msig^8*Mx^2 - 
             15427*Msig^6*Mx^4 - 44308*Msig^4*Mx^6 + 64738*Msig^2*Mx^8 - 
             648*Mx^10)) + 4*M24^4*(1350*MHc^2*Msig^2*(Msig - 2*Mx)*
            (Msig - Mx)*Mx^2*(Msig + Mx)*(Msig + 2*Mx)*(4*Msig^2 + Mx^2) - 
           9*MHc^4*(Msig - 2*Mx)*(Msig - Mx)*(Msig + Mx)*(Msig + 2*Mx)*
            (560*Msig^4 + 93*Msig^2*Mx^2 - 5*Mx^4) - 
           10*(13440*Msig^12 - 80144*Msig^10*Mx^2 + 180016*Msig^8*Mx^4 - 
             342924*Msig^6*Mx^6 + 173249*Msig^4*Mx^8 + 56606*Msig^2*Mx^10 - 
             648*Mx^12)))/(720*M24^2*Msig^2*(M24^2 - 16*Msig^2)*Mx^4*
         (M24^2 - 4*Mx^2)*(Msig^4 - 5*Msig^2*Mx^2 + 4*Mx^4)) + 
       (5*(M24^2 + 8*Msig^2)^2*B0fin[4*Msig^2, M24, 2*Msig, mu])/
        (144*Msig^2*Mx^2) + (-2 + (2*Msig^2)/Mx^2 + (3*Mx^2)/(2*Msig^2))*
        B0fin[4*Msig^2, Mx, Mx, mu] - (7*(MHc^2 - Mx^2)^2*
         B0fin[Mx^2, MHc, 0, mu])/(4*Mx^4) - 
       (35*(M24^4 - 4*M24^2*Mx^2 + 12*Mx^4)*B0fin[Mx^2, Mx, M24, mu])/
        (24*Mx^4) + (14*(-12 - Msig^4/Mx^4 + (4*Msig^2)/Mx^2)*
         B0fin[Mx^2, Mx, Msig, mu])/3 + 
       (-63/2 + (42*Msig^2*(-Msig^2 + Mx^2))/Mx^4)*B0fin[Mx^2, Mx, 2*Msig, 
         mu] + (15*Mx^2*(3*M24^4 - 12*M24^2*Mx^2 + 4*Mx^4)*
         CapitalPhiLambdaM[M24, Mx, Mx])/(4*M24^2*(M24^2 - 4*Mx^2)^3) + 
       (12*Mx^2*(3*Msig^4 - 12*Msig^2*Mx^2 + 4*Mx^4)*CapitalPhiLambdaM[Msig, 
          Mx, Mx])/(Msig^2*(Msig^2 - 4*Mx^2)^3) + 
       (40*Msig^4*(M24^2 + 8*Msig^2)^2*CapitalPhiLambdaM[2*Msig, M24, 
          2*Msig])/(3*M24^2*(M24^2 - 16*Msig^2)^3*Mx^2) - 
       (9*Mx^2*(Msig^2 - 2*Mx^2)*(12*Msig^4 - 12*Msig^2*Mx^2 + Mx^4)*
         CapitalPhiLambdaM[2*Msig, Mx, Mx])/(256*Msig^4*(Msig^2 - Mx^2)^3) + 
       (5*(-108 + (27*M24^2)/(M24 - 2*Mx)^2 + (531*M24)/(M24 - 2*Mx) + 
          (168*M24^4)/Mx^4 - (4*(M24^8 + 114*M24^6*Msig^2 - 
             4320*M24^4*Msig^4 + 32000*M24^2*Msig^6 + 6144*Msig^8))/
           (Msig^2*(M24^2 - 16*Msig^2)^2*Mx^2) + (27*M24^2)/(M24 + 2*Mx)^2 + 
          (531*M24)/(M24 + 2*Mx))*Log[mu^2/M24^2])/576 + 
       (7*(MHc^6 - 3*MHc^4*Mx^2 - MHc^2*Mx^4)*Log[mu^2/MHc^2])/
        (4*Mx^4*(MHc^2 - Mx^2)) + Log[mu^2/(4*Msig^2)]^2/9 - 
       (2*(-7*Msig^8 + 77*Msig^6*Mx^2 - 322*Msig^4*Mx^4 + 468*Msig^2*Mx^6 + 
          72*Mx^8)*Log[mu^2/Msig^2])/(3*Mx^4*(Msig^2 - 4*Mx^2)^2) + 
       Log[mu^2/(4*Msig^2)]*((612 + (81*Msig^2)/(Msig - Mx)^2 + 
           (1269*Msig)/(Msig - Mx) + (24192*Msig^4)/Mx^4 + 
           (16*(5*M24^6 - 1326*M24^4*Msig^2 + 35712*M24^2*Msig^4 - 
              321536*Msig^6))/((M24^2 - 16*Msig^2)^2*Mx^2) - 
           (648*Mx^2)/Msig^2 + (81*Msig^2)/(Msig + Mx)^2 + 
           (1269*Msig)/(Msig + Mx))/576 - 7*Log[mu^2/Mx^2]) - 
       ((70*M24^6*(MHc - Mx)*(MHc + Mx)*(Msig^5 - 5*Msig^3*Mx^2 + 
             4*Msig*Mx^4)^2 + M24^4*(-728*Msig^12*Mx^2 + 16684*Msig^10*Mx^4 - 
            124862*Msig^8*Mx^6 + 378737*Msig^6*Mx^8 - 448378*Msig^4*Mx^10 + 
            178592*Msig^2*Mx^12 - 288*Mx^14 + MHc^2*(728*Msig^12 - 
              16936*Msig^10*Mx^2 + 127382*Msig^8*Mx^4 - 387053*Msig^6*Mx^6 + 
              458458*Msig^4*Mx^8 - 182624*Msig^2*Mx^10 + 288*Mx^12)) + 
          16*Mx^4*(-728*Msig^12*Mx^2 + 16589*Msig^10*Mx^4 - 
            123912*Msig^8*Mx^6 + 375602*Msig^6*Mx^8 - 444578*Msig^4*Mx^10 + 
            177072*Msig^2*Mx^12 - 288*Mx^14 + MHc^2*(728*Msig^12 - 
              16841*Msig^10*Mx^2 + 126432*Msig^8*Mx^4 - 383918*Msig^6*Mx^6 + 
              454658*Msig^4*Mx^8 - 181104*Msig^2*Mx^10 + 288*Mx^12)) + 
          8*M24^2*Mx^2*(MHc^2*(-728*Msig^12 + 16771*Msig^10*Mx^2 - 
              125732*Msig^8*Mx^4 + 381608*Msig^6*Mx^6 - 451858*Msig^4*Mx^8 + 
              179984*Msig^2*Mx^10 - 288*Mx^12) + Mx^2*(728*Msig^12 - 
              16519*Msig^10*Mx^2 + 123212*Msig^8*Mx^4 - 373292*Msig^6*Mx^6 + 
              441778*Msig^4*Mx^8 - 175952*Msig^2*Mx^10 + 288*Mx^12)))*
         Log[mu^2/Mx^2])/(48*Msig^2*(MHc - Mx)*Mx^2*(MHc + Mx)*
         (M24^2 - 4*Mx^2)^2*(Msig^4 - 5*Msig^2*Mx^2 + 4*Mx^4)^2) + 
       (441*Log[mu^2/Mx^2]^2)/4)*UnitStep[-3 + Loops] + 
     alpha4pi*(yt^2*(3/(64*Pi^2) + (3*Log[mu^2/MHc^2])/(16*Pi^2))*
        UnitStep[-3 + Loops] + yb^2*(11/(64*Pi^2) + (3*Log[mu^2/MHc^2])/
          (16*Pi^2))*UnitStep[-3 + Loops] + (1 - Log[mu^2/(4*Msig^2)]/3 + 
	  (21*Log[mu^2/Mx^2])/2)*UnitStep[-2 + Loops])];


(* inverse decoupling coefficient for alpha2 *)
izetaAlSU2[alpha24pi_, yt_, yb_, Mx_, MHc_, Msig_, M24_, mu_, Loops_] := Module[{A,B,izeta,alpha4pi},

     A = zetaAlSU2[alpha24pi,yt,yb,Mx,MHc,Msig,M24,mu,2] - 1;
     alpha4pi = (1-UnitStep[-3 + Loops]*A)*alpha24pi;

     A = zetaAlSU2[alpha4pi,yt,yb,Mx,MHc,Msig,M24,mu,2] - 1;
     B = zetaAlSU2[alpha4pi,yt,yb,Mx,MHc,Msig,M24,mu,3] - A - 1;


     izeta = UnitStep[-3 + Loops]*(A^2 - B) - UnitStep[-2 + Loops]*A + 1;
     izeta
];




(*----------------------------------------------------------------------*)
(*--------------- decoupling coefficient for alpha3 --------------------*)
(*----------------------------------------------------------------------*)

zetaAlSU3[alpha4pi_, yt_, yb_, Mx_, MHc_, Msig_, M24_, mu_, Loops_] := Re[
    1 + alpha4pi^2*(-((-20160*MHc^6*Mx^2)/(-4*MHc^2 + Msig^2)^2 + 
          (34560*MHc^8*Mx^2)/(-4*MHc^2*Msig + Msig^3)^2 + 
          4*MHc^2*(165 + (16*M24^2)/(M24^2 - 4*MHc^2) + (320*Msig^2)/
             (-4*MHc^2 + Msig^2))*Mx^2 + (60*Mx^6)/MHc^2 + 
          20*MHc^4*(-21 + 16*((-M24^2 + 4*MHc^2)^(-1) + 20/(4*MHc^2 - 
                Msig^2) + (9*Msig^2)/(-4*MHc^2 + Msig^2)^2)*Mx^2) - 
          (-25*M24^6*(Msig - 2*Mx)*(Msig - Mx)*(Msig + Mx)*(Msig + 2*Mx)*
             (56*Msig^4 + 99*Msig^2*Mx^2 - 12*Mx^4) + 
            25*M24^8*(14*Msig^6 - 73*Msig^4*Mx^2 + 71*Msig^2*Mx^4 - 
              12*Mx^6) + 4*M24^2*(-11200*Msig^12 + 53895*Msig^10*Mx^2 - 
              72763*Msig^8*Mx^4 + 144892*Msig^6*Mx^6 + 40188*Msig^4*Mx^8 - 
              166352*Msig^2*Mx^10 + 12960*Mx^12 + 1800*Sqrt[3]*Msig^4*Mx^2*(
                Msig^6 - 4*Msig^4*Mx^2 - Msig^2*Mx^4 + 4*Mx^6)*Pi) + 
            2*M24^4*(5600*Msig^10 + 98556*Msig^2*Mx^8 - 6480*Mx^10 - 
              20*Msig^8*Mx^2*(1333 + 45*Sqrt[3]*Pi) - 25*Msig^4*Mx^6*(4399 + 
                144*Sqrt[3]*Pi) + 2*Msig^6*Mx^4*(19277 + 2250*Sqrt[3]*Pi)) + 
            16*Msig^2*Mx^2*(11200*Msig^10 + 179312*Msig^2*Mx^8 - 
              12960*Mx^10 - 7200*Msig^4*Mx^6*(34 + Sqrt[3]*Pi) - 
              5*Msig^8*Mx^2*(13019 + 360*Sqrt[3]*Pi) + Msig^6*Mx^4*(131533 + 
                9000*Sqrt[3]*Pi)))/(Msig^2*(M24^2 - 4*Msig^2)*
            (M24^2 - 4*Mx^2)*(Msig^4 - 5*Msig^2*Mx^2 + 4*Mx^4)))/(360*Mx^4) + 
       (90*Msig^2*S2)/Mx^2 + (8*MHc^2*B0fin[MHc^2, MHc, M24, mu])/(45*Mx^2) + 
       (32*MHc^2*B0fin[MHc^2, MHc, Msig, mu])/(9*Mx^2) + 
       ((MHc^2 - Mx^2)^2*B0fin[MHc^2, Mx, 0, mu])/(6*MHc^2*Mx^2) + 
       (5*(M24^2 + 2*Msig^2)^2*B0fin[Msig^2, M24, Msig, mu])/
        (24*Msig^2*Mx^2) + (4*MHc^4*B0fin[Msig^2, MHc, MHc, mu])/
        (Msig^2*Mx^2) + (-2 + Msig^2/(2*Mx^2) + (6*Mx^2)/Msig^2)*
        B0fin[Msig^2, Mx, Mx, mu] - (7*(MHc^2 - Mx^2)^2*B0fin[Mx^2, MHc, 0, 
          mu])/(6*Mx^4) - (35*(M24^4 - 4*M24^2*Mx^2 + 12*Mx^4)*
         B0fin[Mx^2, Mx, M24, mu])/(36*Mx^4) - 
       (28*(Msig^4 - 4*Msig^2*Mx^2 + 12*Mx^4)*B0fin[Mx^2, Mx, Msig, mu])/
        (9*Mx^4) + 7*(-3 - (4*Msig^4)/Mx^4 + (4*Msig^2)/Mx^2)*
        B0fin[Mx^2, Mx, 2*Msig, mu] + (16*MHc^8*CapitalPhiLambdaM[M24, MHc, 
          MHc])/(15*M24^2*(M24^2 - 4*MHc^2)^3*Mx^2) + 
       (5*Mx^2*(3*M24^4 - 12*M24^2*Mx^2 + 4*Mx^4)*CapitalPhiLambdaM[M24, Mx, 
          Mx])/(2*M24^2*(M24^2 - 4*Mx^2)^3) + 
       (5*Msig^4*(M24^2 + 2*Msig^2)^2*CapitalPhiLambdaM[Msig, M24, Msig])/
        (4*M24^2*(M24^2 - 4*Msig^2)^3*Mx^2) + 
       (8*MHc^8*(18*MHc^2 - Msig^2)*CapitalPhiLambdaM[Msig, MHc, MHc])/
        (3*Msig^4*(-4*MHc^2 + Msig^2)^3*Mx^2) + 
       ((-3*Msig^6*Mx^2 + 66*Msig^4*Mx^4 - 220*Msig^2*Mx^6 + 72*Mx^8)*
         CapitalPhiLambdaM[Msig, Mx, Mx])/(Msig^4*(Msig^2 - 4*Mx^2)^3) + 
       (9*Mx^2*(12*Msig^4 - 12*Msig^2*Mx^2 + Mx^4)*CapitalPhiLambdaM[2*Msig, 
          Mx, Mx])/(128*Msig^2*(Msig^2 - Mx^2)^3) + 
       (((-75*M24^10)/(Msig^2*(M24^2 - 4*Msig^2)^3*Mx^2) + 
          (50*M24^6*((9*Msig^2)/(M24^2 - 4*Msig^2)^3 - 
             77/(M24^2 - 4*Mx^2)^2))/Mx^2 + 
          (175*M24^8*((3*Mx^2)/(M24^2 - 4*Msig^2)^3 + 2/(M24^2 - 4*Mx^2)^2))/
           Mx^4 + M24^4*((4*((-16*MHc^2)/(M24^2 - 4*MHc^2)^2 - (1275*Msig^4)/
                (M24^2 - 4*Msig^2)^3))/Mx^2 + 16100/(M24^2 - 4*Mx^2)^2) + 
          (48*((-8*MHc^6)/(M24^2 - 4*MHc^2)^2 + (150*Msig^8)/
              (M24^2 - 4*Msig^2)^3 - (75*Mx^6)/(M24^2 - 4*Mx^2)^2))/Mx^2 + 
          (8*M24^2*((56*MHc^4)/(M24^2 - 4*MHc^2)^2 + 
             75*((-5*Msig^6)/(M24^2 - 4*Msig^2)^3 - (39*Mx^4)/
                (M24^2 - 4*Mx^2)^2)))/Mx^2)*Log[mu^2/M24^2])/360 + 
       Log[mu^2/MHc^2]^2/36 + ((224*Msig^8 - 616*Msig^6*Mx^2 + 
          644*Msig^4*Mx^4 - 234*Msig^2*Mx^6 - 9*Mx^8)*Log[mu^2/(4*Msig^2)])/
        (8*Mx^4*(Msig^2 - Mx^2)^2) + Log[mu^2/Msig^2]^2/4 + 
       Log[mu^2/Msig^2]*((-3144 + (63*Msig^2)/(Msig - 2*Mx)^2 + 
           (1077*Msig)/(Msig - 2*Mx) + (896*Msig^4)/Mx^4 + 
           (4*(15*M24^6*(-4*MHc^2*Msig + Msig^3)^2 - 2*M24^4*(1728*MHc^8 + 
                192*MHc^6*Msig^2 + 5632*MHc^4*Msig^4 - 3064*MHc^2*Msig^6 + 
                399*Msig^8) + 4*M24^2*Msig^2*(6912*MHc^8 + 768*MHc^6*Msig^2 + 
                17488*MHc^4*Msig^4 - 9736*MHc^2*Msig^6 + 1281*Msig^8) - 
              8*Msig^4*(6912*MHc^8 + 768*MHc^6*Msig^2 + 20848*MHc^4*Msig^4 - 
                11416*MHc^2*Msig^6 + 1491*Msig^8)))/
            (Msig^2*(M24^2 - 4*Msig^2)^2*(-4*MHc^2 + Msig^2)^2*Mx^2) - 
           (1296*Mx^2)/Msig^2 + (63*Msig^2)/(Msig + 2*Mx)^2 + 
           (1077*Msig)/(Msig + 2*Mx))/288 - 7*Log[mu^2/Mx^2]) + 
       Log[mu^2/MHc^2]*((-362 + (105*MHc^4)/Mx^4 + 
           (-1440*MHc^8*(M24^2 - 4*MHc^2)^2 - 16*MHc^6*(335*M24^4 - 2696*
                M24^2*MHc^2 + 5520*MHc^4)*Msig^2 + 8*MHc^4*(265*M24^4 - 2136*
                M24^2*MHc^2 + 4400*MHc^4)*Msig^4 + MHc^2*(-225*M24^4 + 1816*
                M24^2*MHc^2 - 3760*MHc^4)*Msig^6)/((M24^2 - 4*MHc^2)^2*Msig^2*
             (-4*MHc^2 + Msig^2)^2*Mx^2))/90 + Log[mu^2/Msig^2]/6 - 
         (7*Log[mu^2/Mx^2])/3) + 
       (((-70*M24^6 + 95*M24^4*Mx^2 + 560*M24^2*Mx^4)/(M24^2 - 4*Mx^2)^2 + 
          (-12*Mx^4*(Msig^5 - 5*Msig^3*Mx^2 + 4*Msig*Mx^4)^2 + 
            MHc^2*(-728*Msig^12 + 17189*Msig^10*Mx^2 - 126672*Msig^8*Mx^4 + 
              379562*Msig^6*Mx^6 - 453728*Msig^4*Mx^8 + 185376*Msig^2*Mx^10 - 
              1728*Mx^12))/(MHc^2*Msig^2*(Msig^4 - 5*Msig^2*Mx^2 + 4*Mx^4)^
             2))*Log[mu^2/Mx^2])/(72*Mx^2) + 49*Log[mu^2/Mx^2]^2)*
      UnitStep[-3 + Loops] + alpha4pi*
      (yt^2*(11/(128*Pi^2) + (5*Log[mu^2/MHc^2])/(32*Pi^2))*
        UnitStep[-3 + Loops] + yb^2*(19/(128*Pi^2) + (3*Log[mu^2/MHc^2])/
          (16*Pi^2))*UnitStep[-3 + Loops] + (2/3 - Log[mu^2/MHc^2]/6 - 
		Log[mu^2/Msig^2]/2 + 7*Log[mu^2/Mx^2])*UnitStep[-2 + Loops])];

(* inverse decoupling coefficient for alpha3 *)
izetaAlSU3[alpha34pi_, yt_, yb_, Mx_, MHc_, Msig_, M24_, mu_, Loops_] := Module[{A,B,izeta,alpha4pi},

     A = zetaAlSU3[alpha34pi,yt,yb,Mx,MHc,Msig,M24,mu,2] - 1;
     alpha4pi = (1-UnitStep[-3 + Loops]*A)*alpha34pi;

     A = zetaAlSU3[alpha4pi,yt,yb,Mx,MHc,Msig,M24,mu,2] - 1;
     B = zetaAlSU3[alpha4pi,yt,yb,Mx,MHc,Msig,M24,mu,3] - A - 1;


     izeta = UnitStep[-3 + Loops]*(A^2 - B) - UnitStep[-2 + Loops]*A + 1;
     izeta
];


(*----------------------------------------------------------------------*)
(*------------- B0 loop function for various mass hierarchies -----------*)
(*----------------------------------------------------------------------*)

(B0fin[q2_Real,ma_Real,mb_Real,mu_Real] := Module[
    {outstring},
    If[q2 <= (ma-mb)^2, Return[B0fin1[q2,ma,mb,mu^2]]];
    If[(q2 >= (ma-mb)^2) && (q2 <= (ma+mb)^2),Return[B0fin2[q2,ma,mb,mu^2]]];
    If[q2 >= (ma+mb)^2, Return[B0fin3[q2,ma,mb,mu^2]]];
    outstring = "B0fin["<>ToString[q2]<>","<>ToString[ma]<>
	","<>ToString[mb]<>"]: case not implemented";
    Return[outstring];
]);

(* q2 <= (ma-mb)^2: *)

(B0fin1[q2_,ma_,mb_,mu2_] := 
    (2-Log[ma*mb/mu2]+(ma^2-mb^2)/q2*Log[mb/ma]
     + Sqrt[(ma+mb)^2-q2]*Sqrt[(ma-mb)^2-q2]/q2 *
     Log[  (Sqrt[(ma+mb)^2-q2]+Sqrt[(ma-mb)^2-q2]) /
           (Sqrt[(ma+mb)^2-q2]-Sqrt[(ma-mb)^2-q2]) ]
     ));

(* (ma-mb)^2 <= q2 <= (ma+mb)^2: *)

(B0fin2[q2_,ma_,mb_,mu2_] := 
    (2-Log[ma*mb/mu2]+(ma^2-mb^2)/q2*Log[mb/ma]
     - 2*Sqrt[(ma+mb)^2-q2]*Sqrt[q2-(ma-mb)^2]/q2 *
     ArcTan[ Sqrt[q2-(ma-mb)^2]/Sqrt[(ma+mb)^2-q2] ]
     ));

(* q2 >= (ma+mb)^2: *)

(B0fin3[q2_,ma_,mb_,mu2_] :=
    (2-Log[ma*mb/mu2]+(ma^2-mb^2)/q2*Log[mb/ma]
     - Sqrt[q2-(ma+mb)^2]*Sqrt[q2-(ma-mb)^2]/q2 *
     Log[ - (Sqrt[q2-(ma+mb)^2]+Sqrt[q2-(ma-mb)^2]) /
            (Sqrt[q2-(ma+mb)^2]-Sqrt[q2-(ma-mb)^2]) ]
  ));


B0fin[q2_Real, M2_Real, 0, mu_Real] :=
    2 + Log[mu^2/Abs[q2-M2^2]] + M2^2/q2*Log[Abs[q2-M2^2]/M2^2];
B0fin[q2_Real, 0, M2_Real, mu_Real] :=
    2 + Log[mu^2/Abs[q2-M2^2]] + M2^2/q2*Log[Abs[q2-M2^2]/M2^2];

B0fin[0,M1_,M2_,mu_] :=
    1 - (M1^2*Log[M1^2/mu^2] - M2^2*Log[M2^2/mu^2])/(M1^2-M2^2);


(*----------------------------------------------------------------------*)
(*---- two-loop function from Nucl.Phys.B397:123-142,1993 --------------*)
(*---- and other loop functions that appear in the result --------------*)
(*----------------------------------------------------------------------*)
(*------ thanks to Jens Salomon, who implemented these functions! ------*)

Lambda[x_, y_] := Sqrt[(1 - x - y)^2 - 4*x*y];

Cl2[theta_]:=Im[PolyLog[2,Exp[I*theta]]];
S2 = 4/9/Sqrt[3]*Cl2[Pi/3];

CapitalPhi[1][x_, y_] := (1/Lambda[x, y])*(2*Log[(1 + x - y - Lambda[x, y])/2]*
      Log[(1 - x + y - Lambda[x, y])/2] - Log[x]*Log[y] -
     2*PolyLog[2, (1 + x - y - Lambda[x, y])/2] -
     2*PolyLog[2, (1 - x + y - Lambda[x, y])/2] + (1/3)*Pi^2);

CapitalPhi[2][x_, y_] := (2/Sqrt[-Lambda[x, y]^2])*(Cl2[2*ArcCos[(-1 + x + y)/(2*Sqrt[x*y])]] +
     Cl2[2*ArcCos[(1 + x - y)/(2*Sqrt[x])]] +
     Cl2[2*ArcCos[(1 - x + y)/(2*Sqrt[y])]]);

CapitalPhiLambdaM[M1_Real, M2_Real, M3_Real] :=
Module[{errorstring = "Case not implemented"},
  If[Re[Lambda[M1^2/M3^2, M2^2/M3^2]^2] < 0, Return[M3^2*Lambda[M1^2/M3^2,M2^2/M3^2]^2*CapitalPhi[2][M1^2/M3^2,M2^2/M3^2]]];
  If[M1 + M2 <= M3, Return[M3^2*Lambda[M1^2/M3^2,M2^2/M3^2]^2*CapitalPhi[1][M1^2/M3^2,M2^2/M3^2]]];
  If[M2 + M3 <= M1, Return[M1^2*Lambda[M2^2/M1^2,M3^2/M1^2]^2*CapitalPhi[1][M2^2/M1^2,M3^2/M1^2]]];
  If[M1 + M3 <= M2, Return[M2^2*Lambda[M1^2/M2^2,M3^2/M2^2]^2*CapitalPhi[1][M1^2/M2^2,M3^2/M2^2]]];
  errorstring
]


TfinZEROaux[M2_,M3_] := + 2 M2^2 Log[M2/M3]^2 +
M3^2/2*(2 (1-y) Log[1-y] Log[y] + 2 (1-y) PolyLog[2,y] - (1-y) Pi^2/3)/.{y->M2^2/M3^2};

TfinZEROshort[M2_Real, M3_Real] :=
    Which[M2 < M3, TfinZEROaux[M2, M3],
	  M3 < M2, TfinZEROaux[M3, M2],
	  M2 == M3,0];


d123[M1_, M2_, M3_] := -1/(2*(M1^2*M2^2 + M1^2*M3^2 + M2^2*M3^2) - M1^4 - M2^4 - M3^4);



End[];

EndPackage[];
