1# evolve the RGEs of the standard model from electroweak scale up
2# by dpgeorge
3
4import math
5
6
7class RungeKutta(object):
8    def __init__(self, functions, initConditions, t0, dh, save=True):
9        self.Trajectory, self.save = [[t0] + initConditions], save
10        self.functions = [lambda *args: 1.0] + list(functions)
11        self.N, self.dh = len(self.functions), dh
12        self.coeff = [1.0 / 6.0, 2.0 / 6.0, 2.0 / 6.0, 1.0 / 6.0]
13        self.InArgCoeff = [0.0, 0.5, 0.5, 1.0]
14
15    def iterate(self):
16        step = self.Trajectory[-1][:]
17        istep, iac = step[:], self.InArgCoeff
18        k, ktmp = self.N * [0.0], self.N * [0.0]
19        for ic, c in enumerate(self.coeff):
20            for if_, f in enumerate(self.functions):
21                arguments = [(x + k[i] * iac[ic]) for i, x in enumerate(istep)]
22                try:
23                    feval = f(*arguments)
24                except OverflowError:
25                    return False
26                if abs(feval) > 1e2:  # stop integrating
27                    return False
28                ktmp[if_] = self.dh * feval
29            k = ktmp[:]
30            step = [s + c * k[ik] for ik, s in enumerate(step)]
31        if self.save:
32            self.Trajectory += [step]
33        else:
34            self.Trajectory = [step]
35        return True
36
37    def solve(self, finishtime):
38        while self.Trajectory[-1][0] < finishtime:
39            if not self.iterate():
40                break
41
42    def solveNSteps(self, nSteps):
43        for i in range(nSteps):
44            if not self.iterate():
45                break
46
47    def series(self):
48        return zip(*self.Trajectory)
49
50
51# 1-loop RGES for the main parameters of the SM
52# couplings are: g1, g2, g3 of U(1), SU(2), SU(3); yt (top Yukawa), lambda (Higgs quartic)
53# see arxiv.org/abs/0812.4950, eqs 10-15
54sysSM = (
55    lambda *a: 41.0 / 96.0 / math.pi ** 2 * a[1] ** 3,  # g1
56    lambda *a: -19.0 / 96.0 / math.pi ** 2 * a[2] ** 3,  # g2
57    lambda *a: -42.0 / 96.0 / math.pi ** 2 * a[3] ** 3,  # g3
58    lambda *a: 1.0
59    / 16.0
60    / math.pi ** 2
61    * (
62        9.0 / 2.0 * a[4] ** 3
63        - 8.0 * a[3] ** 2 * a[4]
64        - 9.0 / 4.0 * a[2] ** 2 * a[4]
65        - 17.0 / 12.0 * a[1] ** 2 * a[4]
66    ),  # yt
67    lambda *a: 1.0
68    / 16.0
69    / math.pi ** 2
70    * (
71        24.0 * a[5] ** 2
72        + 12.0 * a[4] ** 2 * a[5]
73        - 9.0 * a[5] * (a[2] ** 2 + 1.0 / 3.0 * a[1] ** 2)
74        - 6.0 * a[4] ** 4
75        + 9.0 / 8.0 * a[2] ** 4
76        + 3.0 / 8.0 * a[1] ** 4
77        + 3.0 / 4.0 * a[2] ** 2 * a[1] ** 2
78    ),  # lambda
79)
80
81
82def drange(start, stop, step):
83    r = start
84    while r < stop:
85        yield r
86        r += step
87
88
89def phaseDiagram(system, trajStart, trajPlot, h=0.1, tend=1.0, range=1.0):
90    tstart = 0.0
91    for i in drange(0, range, 0.1 * range):
92        for j in drange(0, range, 0.1 * range):
93            rk = RungeKutta(system, trajStart(i, j), tstart, h)
94            rk.solve(tend)
95            # draw the line
96            for tr in rk.Trajectory:
97                x, y = trajPlot(tr)
98                print(x, y)
99            print()
100            # draw the arrow
101            continue
102            l = (len(rk.Trajectory) - 1) / 3
103            if l > 0 and 2 * l < len(rk.Trajectory):
104                p1 = rk.Trajectory[l]
105                p2 = rk.Trajectory[2 * l]
106                x1, y1 = trajPlot(p1)
107                x2, y2 = trajPlot(p2)
108                dx = -0.5 * (y2 - y1)  # orthogonal to line
109                dy = 0.5 * (x2 - x1)  # orthogonal to line
110                # l = math.sqrt(dx*dx + dy*dy)
111                # if abs(l) > 1e-3:
112                #    l = 0.1 / l
113                #    dx *= l
114                #    dy *= l
115                print(x1 + dx, y1 + dy)
116                print(x2, y2)
117                print(x1 - dx, y1 - dy)
118                print()
119
120
121def singleTraj(system, trajStart, h=0.02, tend=1.0):
122    tstart = 0.0
123
124    # compute the trajectory
125
126    rk = RungeKutta(system, trajStart, tstart, h)
127    rk.solve(tend)
128
129    # print out trajectory
130
131    for i in range(len(rk.Trajectory)):
132        tr = rk.Trajectory[i]
133        print(" ".join(["{:.4f}".format(t) for t in tr]))
134
135
136# phaseDiagram(sysSM, (lambda i, j: [0.354, 0.654, 1.278, 0.8 + 0.2 * i, 0.1 + 0.1 * j]), (lambda a: (a[4], a[5])), h=0.1, tend=math.log(10**17))
137
138# initial conditions at M_Z
139singleTraj(
140    sysSM, [0.354, 0.654, 1.278, 0.983, 0.131], h=0.5, tend=math.log(10 ** 17)
141)  # true values
142