curvilinear coordinates

[Staffan Lundberg]

I am working with curvilinear coordinates, especiallt how to determine scale factors for spherical coordinates and how to express the unit vectors i terms of i, j, k (cartesian base vectors). I submit some code

#

# curvilinear.py

#from sympy import *
x, y, z = symbols("x y z")
rho= symbols("rho",positive=True,real=True)
theta = symbols("theta",positive=True,real=True)

phi = symbols("phi",positive=True,real=True)
# spherical coord
x=rho*sin(phi)*cos(theta)
y=rho*sin(phi)*sin(theta)
z=rho*cos(phi)

A=Matrix([[x],[y],[z]])
r=A
dr=diff(r,rho)
hr=dr.norm()
hr=simplify(hr)

dphi=diff(r,phi)
hphi=dphi.norm()
hphi=simplify(hphi)

dtheta=diff(r,theta)
htheta=dtheta.norm()
htheta=simplify(htheta)

r_hat=dr/hr
fi_hat=dphi/hphi
th_hat=dtheta/htheta
print(fi_hat)
print(th_hat)
print(r_hat)

#

#

Problems with th_hat. Despite telling sympy that phi is positive,  I get

Abs(sin(phi))  instead of sin(phi).  Thus python does not cancel the factor sin(phi).

Has anyone some hints how to solve this issue. Maybe a bug in sympy?

/Staffan L

[Arthur Ryman]

Staffan,

Just a guess, but sin(phi) goes negative even for positive values of phi. You’d have limit the range to 0 <= phi  <= pi. 

— Arthur 

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/717dfd84-7b24-4a28-83c5-19ba7ca53e97n%40googlegroups.com.

[Chris Smith]

If you want to ignore `Abs`, replace it: `expr.replace(Abs, Id)`

/c

[Alan Bromborsky]

The following code may be useful -

def square_root_of_expr(expr):
    """
    If expression is product of even powers then every power is divided
    by two and the product is returned.  If some terms in product are
    not even powers the sqrt of the absolute value of the expression is
    returned.  If the expression is a number the sqrt of the absolute
    value of the number is returned.
    """
    if expr.is_number:
        if expr > 0:
            return sqrt(expr)
        else:
            return sqrt(-expr)
    else:
        expr = trigsimp(expr)
        coef, pow_lst = sqf_list(expr)
        if coef != S.One:
            if coef.is_number:
                coef = square_root_of_expr(coef)
            else:
                coef = sqrt(abs(coef))  # Product coefficient not a number
        for p in pow_lst:
            f, n = p
            if n % 2 != 0:
                return sqrt(abs(expr))  # Product not all even powers
            else:
                coef *= f ** (n / S(2))  # Positive sqrt of the square of an expression
        return coef

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/ad627534-58eb-4001-a76f-2aad0d6f16e8n%40googlegroups.com.