Get gradient in another coordinate system by transformation

get_grad.py

from sympy import *
init_printing()
from sympy.diffgeom import *
r, x, y, z = symbols('r x y z', positive=True)
theta, phi = symbols('theta phi')
N = 3
M = Manifold('M', N)
P = Patch('P', M)
rect = CoordSystem('rect', P, ['x', 'y', 'z'])

def get_grad(from_transformation, syms=['r','theta','phi']):
    newcs = CoordSystem('newcs', P, syms)
    newsymbols = symbols(' '.join(syms), positive=True)
    newcs.connect_to(rect, newsymbols, from_transformation(*newsymbols), inverse=False)
    #rect.connect_to(newcs, [x, y, z], to_transformation(x, y, z), inverse=False)
    jac2to1 = newcs.jacobian(rect, newsymbols)
    new_metric = simplify(jac2to1.inv() * jac2to1.inv().T)    
    scale_factors = [sqrt(new_metric[i,i]) for i in range(N)]
    
    # get gradient
    unscaled_grad = new_metric*Matrix([1]*N)  # Matrix(newcs.base_vectors())
    grad = [simplify(unscaled_grad[i]/scale_factors[i])*newcs.base_vector(i) for i in range(N)]
    
    return grad

    # TO BE COMPLETED: get divergence
    # following formula: (8.2.1)
    sqrtabsg = simplify(sqrt((new_metric.det())))
    print sqrtabsg
    def div(f, sqrtabsg):
        print 1/sqrtabsg*f[0]
        sqrtabsg = sqrtabsg.xreplace(dict(zip(newsymbols, newcs.coord_functions())))
        import pdb; pdb.set_trace()
        return sqrtabsg * sum([newcs.base_vector(i)(1/sqrtabsg*f[i]) for i in range(N)])
    return grad  # newcs.coord_functions(), grad, lambda f: div(f, sqrtabsg)


from_spherical = lambda r, theta, phi: (r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta))
grad = get_grad(from_spherical)