In [44]:

    

from sympy import symbols, Function, Rational,Matrix,cse

def zparam():
    px,py,pz,x,y,z,s = symbols('px py pz x y z s')
    n = Function('n')(x,y,z)
    nx = n.diff(x)
    ny = n.diff(y)
    nz = n.diff(z)
    sdot = n/pz
    pxdot = nx*n/pz
    pydot = ny*n/pz
    pzdot = nz*n/pz
    xdot = px/pz
    ydot = py/pz
    zdot = Rational(1)
    euVec = Matrix([pxdot,pydot,pzdot,xdot,ydot,zdot,sdot]).T
    jac = Matrix([euVec.diff(px)[:],
         euVec.diff(py)[:],
         euVec.diff(pz)[:],
         euVec.diff(x)[:],
         euVec.diff(y)[:],
         euVec.diff(z)[:],
         euVec.diff(s)[:]])
    cseFunc = cse(jac.T,optimizations='basic')
    print(cseFunc)
    
def sparam():
    px,py,pz,x,y,z,s = symbols('px py pz x y z s')
    n = Function('n')(x,y,z)
    nx = n.diff(x)
    ny = n.diff(y)
    nz = n.diff(z)
    sdot = Rational(1)
    pxdot = nx
    pydot = ny
    pzdot = nz
    xdot = px/n
    ydot = py/n
    zdot = pz/n
    euVec = Matrix([pxdot,pydot,pzdot,xdot,ydot,zdot,sdot]).T
    jac = Matrix([euVec.diff(px)[:],
         euVec.diff(py)[:],
         euVec.diff(pz)[:],
         euVec.diff(x)[:],
         euVec.diff(y)[:],
         euVec.diff(z)[:],
         euVec.diff(s)[:]])
    cseFunc = cse(jac.T,optimizations='basic')
    print(cseFunc)
    
if __name__=="__main__":
    zparam()
    sparam()


    

    




([(x0, n(x, y, z)), (x1, Derivative(x0, x)), (x2, pz**(-2)), (x3, x0*x2), (x4, 1/pz), (x5, Derivative(x0, y)), (x6, x4*(x0*Derivative(x0, x, y) + x1*x5)), (x7, Derivative(x0, z)), (x8, x4*(x0*Derivative(x0, x, z) + x1*x7)), (x9, x4*(x0*Derivative(x0, y, z) + x5*x7))], [Matrix([
[ 0,  0, -x1*x3, x4*(x0*Derivative(x0, x, x) + x1**2),                                   x6,                                   x8, 0],
[ 0,  0, -x3*x5,                                   x6, x4*(x0*Derivative(x0, y, y) + x5**2),                                   x9, 0],
[ 0,  0, -x3*x7,                                   x8,                                   x9, x4*(x0*Derivative(x0, z, z) + x7**2), 0],
[x4,  0, -px*x2,                                    0,                                    0,                                    0, 0],
[ 0, x4, -py*x2,                                    0,                                    0,                                    0, 0],
[ 0,  0,      0,                                    0,                                    0,                                    0, 0],
[ 0,  0,    -x3,                                x1*x4,                                x4*x5,                                x4*x7, 0]])])
([(x0, n(x, y, z)), (x1, Derivative(x0, x, y)), (x2, Derivative(x0, x, z)), (x3, Derivative(x0, y, z)), (x4, 1/x0), (x5, Derivative(x0, x)), (x6, x0**(-2)), (x7, px*x6), (x8, Derivative(x0, y)), (x9, Derivative(x0, z)), (x10, py*x6), (x11, pz*x6)], [Matrix([
[ 0,  0,  0, Derivative(x0, x, x),                   x1,                   x2, 0],
[ 0,  0,  0,                   x1, Derivative(x0, y, y),                   x3, 0],
[ 0,  0,  0,                   x2,                   x3, Derivative(x0, z, z), 0],
[x4,  0,  0,               -x5*x7,               -x7*x8,               -x7*x9, 0],
[ 0, x4,  0,              -x10*x5,              -x10*x8,              -x10*x9, 0],
[ 0,  0, x4,              -x11*x5,              -x11*x8,              -x11*x9, 0],
[ 0,  0,  0,                    0,                    0,                    0, 0]])])

In [ ]: