import sympy as sp
c2,KL,G,theta=sp.symbols("c2 KL G theta", nonzero=True)
pd,T=sp.symbols("phidot Theta")
# momentum: pd-c2*T/2=4*pi*G*theta
# alpha at O(v^3): c2*T=KL*pd
lhs=sp.simplify(pd-c2*(KL*pd/c2)/2)
pdsol=sp.solve(sp.Eq(lhs,4*sp.pi*G*theta),pd)[0]
# integrated Phi equation: lap zeta=3 pd -3c2 T/2
lapz=sp.simplify(3*pdsol-sp.Rational(3,2)*c2*(KL*pdsol/c2))
assert sp.simplify(lapz-12*sp.pi*G*theta)==0
print("phidot =",pdsol)
print("laplacian(zeta) =",lapz)
print("c2 cancels from C_L exactly at O(v^3)")
print("PASS")