from sympy import *

print("This determines the function exp(psi) associated to the cubic")
print("Lefschetz pencil on CP^2, and the other coefficients in the")
print("differential equation. It then proceeds to solve the equation.")
print()

q = symbols('q')
f = (1-q**3)*(1-q**6)*(1-q**9)*(1-q**12)*(1-q**15)*(1-q**18)*(1-q**21)*(1-q**24)*(1-q**27)*(1-q**30)*(1-q**33)*(1-q**36)
exppsi = series(f**4,q,n=39)
print("exp(psi) is")
print(exppsi)
print()

eta = series(-diff(exppsi, q)/exppsi,q,n=38)
print("eta = -d psi/dq is")
print(eta)
print()

lambd = (-1 - 7*q**3 - 35*q**6 - 103*q**9 - 312*q**12 - 819*q**15 - 1996*q**18 - 4591*q**21 - 10122*q**24 - 21299*q**27 - 43501*q**30 - 86090*q**33 + O(q**36))/q

zq0 = series(q*q*(lambd*lambd - diff(lambd/exppsi,q))/4,q,n=36)/(q*q)
zq = expand(zq0,q,n=34)
print("the count of rulings, z^{(2)}, is")
print(zq)
print()

def iterate(listofterms):
  a = listofterms[0]
  b = listofterms[1]
  c = listofterms[2]
  d = listofterms[3]
  a0 = expand(-b*exppsi,q)
  b0 = expand(-4*exppsi*zq*a-eta*b,q)
  c0 = expand(-d*exppsi,q)
  d0 = expand(-4*exppsi*zq*c-eta*d,q)
  a1 = 1 + integrate(a0,q)
  b1 = integrate(b0,q)
  c1 = integrate(c0,q)
  d1 = 1 + integrate(d0,q)
  return([a1,b1,c1,d1])

print("We now start solving the ODE iteratively")
print()
s = [1,0,0,1]
for i in range(0,34):
  s0 = iterate(s)
  print("iteration", i)
  s = s0
print("Solution:")
print(s)
print()

qu = series(s[2]/s[0],q,n=30)
aa = -series(q/qu,q,n=28)
print("Where the fibre rho = 0 is in terms of q:")
print()
print(aa/q)
print()

jinvariant = aa**3*(aa**3-24*q**3)**3/(aa**3 - 27*q**3)
jhesse = series(jinvariant,q,n=28)/q**9
print("Plug it into the j-invariant of the mirror of P^2:")
print(jhesse)