#!/usr/bin/env python
# coding: utf-8

# ## The differential equation for the two-loop Ice-cream cone integral in general dimensions

# \begin{equation}
#   \omega(t)= {\textbf{U}\over
#     \textbf{F}^2} \left(\textbf{U}^3\over  \textbf{F}^2\right)^{\epsilon}
# \end{equation}
# 
# with
# \begin{align}
#   \textbf{U}&:=(y_1+y_2)(x_1+z)+zx_1,\cr
#                         \textbf{V}&:=k_2^2y_1y_2(z+x_1)+zx_1(k_1^2y_1+k_3^2y_2),\\
# \nonumber  \textbf{F}&:=  (\mu_1^2y_1+\mu_2^2y_2^2+m_1^2x_1+m_2^2z)\textbf{U}-t \textbf{V}.
# \end{align}
# 
# In this worksheet we present some numerical cases for the $\epsilon$-dependent differential operator with respect to $t$ for the ice-cream cone graph

# This is the generic two-loop Ice-cream cone  differential equation with $\epsilon=0$ derived in [arXiv:2209.10962](https://arxiv.org/abs/2209.10962) and [arXiv:2302.14840](https://arxiv.org/abs/2302.14840)
# The differential operator is
# $$\textrm{PFIceCreamCoefficient}[2]\,Dt^2+\textrm{PFIceCreamCoefficient}[1]\,Dt+\textrm{PFIceCreamCoefficient}[0]$$

# In[1]:


PFIceCreamCoefficient = dict()
with open("PFIcecream-Coefficients.txt") as f:
    for e in f.read().replace("\n", "").split(";"):
        if e:
            name, expr = e.split(":=")
            name = name.strip()
            expr = SR(expr)
            if name[0] == "c":
                PFIceCreamCoefficient[int(name[1:])] = expr


# In[2]:


from ore_algebra import *
OA, t, Dt = DifferentialOperators(QQ, 't')


# In[3]:


Dt,t,epsilon= var('Dt t epsilon')


# <b>The all equal parameters $m_1=m_2=\mu_1=\mu_2=k_1^2=k_2^2=k_3^2=1$</b>

# In[4]:


PFAllEqual=dict()
PFAllEqual[0]=(-24*t^3+38*t^4-16*t^5+2*t^6)*Dt^3+(-60*t^2+206*t^3-132*t^4+22*t^5)*Dt^2+(178*t^2-232*t^3+58*t^4)*Dt+(-64*t^2+32*t^3)
PFAllEqual[1]=(-12*t^3+7*t^4+11*t^5-7*t^6+t^7)*Dt^3+(18*t^2+95*t^3-26*t^4-37*t^5+10*t^6)*Dt^2+(48*t+53*t^2-242*t^3+5*t^4+24*t^5)*Dt+(-12-48*t-112*t^2+64*t^3+12*t^4)
PFAllEqual[2]=(24*t^2+29*t^3-17*t^4-17*t^5+5*t^6)*Dt^2+(-71*t^2-130*t^3-23*t^4+28*t^5)*Dt+(-18-64*t-48*t^2+56*t^3+26*t^4)
PFAllEqual[3]=(-12*t-40*t^2-36*t^3+8*t^5)*Dt+(-6-12*t+12*t^2+36*t^3+18*t^4)
PFAllEqual[4]=4*t + 12*t^2 + 12*t^3 + 4*t^4


# In[5]:


facPFAllEqual0=OA(PFAllEqual[0]).factor()
[x.order() for x in facPFAllEqual0]


# In[7]:


PFnum=(PFIceCreamCoefficient[2]*Dt^2+PFIceCreamCoefficient[1]*Dt+PFIceCreamCoefficient[0]).subs(m1=1,m2=1,m3=1,m4=1,p12=1,p22=1,p32=1)


# In[8]:


facPFnum=OA(PFnum).factor()
[x.order() for x in facPFnum]


# In[9]:


facPFAllEqual0[2]-9*facPFnum[1]


# In[10]:


18*OA((t-1)^2)*OA((t-3)^2)*facPFAllEqual0[1]-facPFnum[0]


# <b>The all equal mass case $m_1=m_2=\mu_1=\mu_2=1$ with $|k_1|=3$, $|k_2|=83$ and $|k_3|=71$</b>

# In[11]:


PF2IceCreamAllEqualMass=dict()
with open("PFIceCream-All-Equal-Mass.txt") as f:
    for e in f.read().replace("\n", "").replace(" PF","PF").split(";"):
            if e:
                name, expr = e.split(" := ")
                if name[0] == "P":
                    _,num=name.split("PF")
                    print("loading epsilon order", num)
                    PF2IceCreamAllEqualMass[int(num)] =  SR(expr)


# In[12]:


[OA(PF2IceCreamAllEqualMass[i]).order() for i in range(5)]


# In[13]:


facPF2IceCreamAllEqualMass=OA(PF2IceCreamAllEqualMass[0]).factor()
[x.order() for x in facPF2IceCreamAllEqualMass]


# In[14]:


exp1=facPF2IceCreamAllEqualMass[1].leading_coefficient()


# In[15]:


exp2=PFIceCreamCoefficient[2].subs(m1=1,m2=1,m3=1,m4=1,p12=3^2,p22=83^2,p32=71^2)


# In[15]:


norm=(exp1/exp2).simplify_rational()


# In[16]:


PFnum=(PFIceCreamCoefficient[2]*Dt^2+PFIceCreamCoefficient[1]*Dt+PFIceCreamCoefficient[0]).subs(m1=1,m2=1,m3=1,m4=1,p12=3^2,p22=83^2,p32=71^2)


# In[90]:


facPF2IceCreamAllEqualMass[1]-OA(norm*PFnum)


# <b>The all equal scoop mass case $m_1=m_2=1$ with $\mu_1=59$,  $\mu_2=5$,  $|k_1|=89$, $|k_2|=2$ and $|k_3|=11$</b>

# In[33]:


PF2IceCreamScoopEqualMass=dict()
with open("PFIceCream-Scoop-Equal-Mass.txt") as f:
    for e in f.read().replace("\n", "").replace(" PF","PF").split(";"):
            if e:
                name, expr = e.split(" := ")
                if name[0] == "P":
                    _,num=name.split("PF")
                    print("loading epsilon order", num)
                    PF2IceCreamScoopEqualMass[int(num)] =  SR(expr)


# In[34]:


[OA(PF2IceCreamScoopEqualMass[i]).order() for i in range(5)]


# In[93]:


facPF2IceCreamScoopEqualMass=OA(PF2IceCreamScoopEqualMass[0]).factor()
[x.order() for x in facPF2IceCreamScoopEqualMass]


# In[104]:


exp1=facPF2IceCreamScoopEqualMass[1].leading_coefficient()
exp2=PFIceCreamCoefficient[2].subs(m3=1,m4=1,m1=sqrt(59),m2=sqrt(5),p12=89^2,p22=2^2,p32=11^2)
norm=(exp1/exp2).simplify_rational(); norm


# In[109]:


PFnum=(PFIceCreamCoefficient[2]*Dt^2+PFIceCreamCoefficient[1]*Dt+PFIceCreamCoefficient[0]).subs(m3=1,m4=1,m1=sqrt(59),m2=sqrt(5),p12=89^2,p22=2^2,p32=11^2).expand()


# In[110]:


facPF2IceCreamScoopEqualMass[1]-OA(norm*PFnum)


# <b>Generic mass configuration  $|k_1|=43$, $|k_2|=5$, $|k_3|=89$,  $\mu_1=53$,  $\mu_2=23$, $m_1=7$ and $m_2=47$</b>

# In[35]:


PF2IceCreamGeneric=dict()
with open("PFIceCream-Generic.txt") as f:
    for e in f.read().replace("\n", "").replace(" PF","PF").split(";"):
            if e:
                name, expr = e.split(" := ")
                if name[0] == "P":
                    _,num=name.split("PF")
                    print("loading epsilon order", num)
                    PF2IceCreamGeneric[int(num)] =  SR(expr)


# In[36]:


[OA(PF2IceCreamGeneric[i]).order() for i in range(5)]


# <b>This is the factorisation of the $\epsilon^0$ term</b>

# In[7]:


facPF2IceCreamGeneric=OA(PF2IceCreamGeneric[0]).factor(verbose=True)


# In[8]:


[x.order() for x in facPF2IceCreamGeneric]


# In[128]:


#19, 7, 5, 61, 31, 61, 31


# In[134]:


exp1=facPF2IceCreamGeneric[2].leading_coefficient()
exp2=PFIceCreamCoefficient[2].subs(m3=61,m4=31,m1=sqrt(61),m2=sqrt(31),p12=19^2,p22=7^2,p32=5^2)
norm=(exp1/exp2).simplify_rational(); norm


# In[137]:


PFnum=(PFIceCreamCoefficient[2]*Dt^2+PFIceCreamCoefficient[1]*Dt+PFIceCreamCoefficient[0]).subs(m3=61,m4=31,m1=sqrt(61),m2=sqrt(31),p12=19^2,p22=7^2,p32=5^2).expand()


# In[139]:


facPF2IceCreamGeneric[2]-OA(norm*PFnum)