[Python] Unable to collect terms to reach an expression with Sympy

I have the expression

Delta_T*(-S_xx*r**4 - S_xx*r**2*ri**2 - S_xx*r**2*ro**2 + S_xx*ri**2*ro**2 + S_yy*r**4 + S_yy*r**2*ri**2 + S_yy*r**2*ro**2 - S_yy*ri**2*ro**2)*cos(2*theta)/(4*r**2*(ri**2 + ro**2)*log(ro/ri)) - f(r) + g(theta) + (Delta_T*(S_xx + S_yy)*(ri**2 + ro**2)*log(r) - r**2*(Delta_T*S_xx*sin(theta)**2 - Delta_T*S_xx/2 - Delta_T*S_yy*sin(theta)**2 + Delta_T*S_yy/2) - (-2*Delta_T*S_xx*ri**2*ro**2*sin(theta)**2 + Delta_T*S_xx*ri**2*ro**2 + 2*Delta_T*S_yy*ri**2*ro**2*sin(theta)**2 - Delta_T*S_yy*ri**2*ro**2)/(2*r**2))/(2*ri**2*log(ro/ri) + 2*ro**2*log(ro/ri))


which should be equivalent to g(theta) - f(r) + (2 * Delta_T * sp.log(r) * (S_xx + S_yy) + Delta_T * sp.cos(2 * theta) * (S_yy - S_xx)) / (4 * sp.log(ro / ri)) (there might be a minus sign difference, I am not sure, but it doesn't matter because these expressions are supposed to be worth 0)

This means I should be able to separate r from theta. Or, in Sympy's words, r from sin(theta)**2, cos(2*theta), but I wasn't able to progress. I tried an approach with collect() to no avail.

For example expr_collected = sp.collect(V_r_integral - V_theta_integral, [r, sp.cos(theta)**2, sp.sin(theta)**2]) yields

Delta_T*(S_xx*ri**2*ro**2 - S_yy*ri**2*ro**2 + r**4*(-S_xx + S_yy) + r**2*(-S_xx*ri**2 - S_xx*ro**2 + S_yy*ri**2 + S_yy*ro**2))*cos(2*theta)/(4*r**2*(ri**2 + ro**2)*log(ro/ri)) - f(r) + g(theta) + (Delta_T*(S_xx + S_yy)*(ri**2 + ro**2)*log(r) + r**2*(Delta_T*S_xx/2 - Delta_T*S_yy/2 - (Delta_T*S_xx - Delta_T*S_yy)*sin(theta)**2) + (-Delta_T*S_xx*ri**2*ro**2/2 + Delta_T*S_yy*ri**2*ro**2/2 - (-2*Delta_T*S_xx*ri**2*ro**2 + 2*Delta_T*S_yy*ri**2*ro**2)*sin(theta)**2/2)/r**2)/(2*ri**2*log(ro/ri) + 2*ro**2*log(ro/ri))

and we see terms with r squared multiplying terms with cos(2*theta). Same with expr_collected = sp.collect(V_r_integral - V_theta_integral, [r**2, sp.cos(theta)**2, sp.sin(theta)**2, sp.cos(2*theta)])

Here's a working example, the last line being the expression I want to recast:

import sympy as sp
r, theta, H0, ri, ro, rho, T_0, Delta_T, S_xx, S_yy = sp.symbols('r theta H0 ri ro rho T_0 Delta_T S_xx S_yy')
J_r = -2 * H0 * (1/r - 1/(ri**2 + ro**2)*(r + ri**2*ro**2/r**3)) * sp.cos(2*theta)
J_theta = -H0 * (1/(ri**2 + ro**2)*(2*r - 2*ri**2*ro**2/r**3)) * sp.sin(2*theta)
J_r = J_r.subs(H0, -(S_xx - S_yy) / (4 * rho) * Delta_T / sp.log(ro/ri))
J_theta = J_theta.subs(H0, -(S_xx - S_yy) / (4 * rho) * Delta_T / sp.log(ro/ri))
T = T_0 - Delta_T * sp.log(r/ro) / sp.log(ro/ri)
grad_T_r = sp.diff(T, r)
grad_T_theta = sp.diff(T, theta)
R = sp.Matrix([[sp.cos(theta), sp.sin(theta)],[-sp.sin(theta), sp.cos(theta)]])
S_cartesian = sp.Matrix([[S_xx, 0], [0, S_yy]])
S_polar = R*S_cartesian*R.T
S_rr = S_polar[0,0]
S_rtheta = S_polar[0,1]
S_thetar= S_polar[1,0]
S_thetatheta= S_polar[1,1]
V = sp.Function('V')(r, theta)
grad_V_r = sp.diff(V, r)
grad_V_theta = 1/r * sp.diff(V, theta)
sigma = 1/rho
eq_r = sp.Eq(J_r, -sigma * grad_V_r -sigma * (S_rr * grad_T_r + S_rtheta * grad_T_theta))
eq_theta = sp.Eq(J_theta, -sigma * grad_V_theta - sigma * (S_thetar * grad_T_r + S_thetatheta * grad_T_theta))
sol_grad_V_r = sp.solve(eq_r, grad_V_r)[0]
sol_grad_V_theta = sp.solve(eq_theta, grad_V_theta)[0]
f = sp.Function('f')(r)
g = sp.Function('g')(theta)
V_r_integral = sp.integrate(sol_grad_V_r, r) + g
V_theta_integral = sp.integrate(sol_grad_V_theta * r, theta) + f
print(V_r_integral - V_theta_integral)

Continue reading...