Vector Calculus functions for rectangular, cylindrical, and spherical coordinates

[hamorabi]

I'm sharing a couple of the most common functions in vector calculus. As far as I'm aware, the prime already has these functions but works only for rectangular coordinates. As most people that have studied electromagnetics know, cylindrical and spherical coordinates are just as widely used.

There are currently 7 functions included. I chose to included only rectangular, cylindrical, and spherical coordinate systems. I could have easily included other coordinate systems like parabolic, ellipsoidal, etc, but I didn't use them during my studies, so didn't include them. In case anyone is wondering, the code is based on the general curvilinear coordinate systems.

Just a word on general usage.

The first parameter is always either a vector-valued function presented in a list or just a scalar function. Functions with "vlist" take a vector-valued function, while those with ff take scalar functions.

The second parameter, "vars", is the list of variables used.

The third parameter, "csys", is the coordinate system used. To use cylindrical coordinates, you can write "cylinder", or since no one wants to type that, just enter 2. 1 is for rectangular, 2 is for cylindrical, and 3 is for spherical.

The fourth parameter, "uu", is optional. You include it if you want want to evaluate the result at some point.

For example, let's say we want to find the curl in cylindrical coordinates of the vector field with components 2*z^2 in the r direction, 5*r in the theta direction, and -3*r^2 in the z direction, we enter:

PHP Code:

curl2([2*z^2,5*r,-3*r^2],[r,t,z],2) 


If we want to evaluate the result at some point, say (1,1,1):

PHP Code:

curl2([2*z^2,5*r,-3*r^2],[r,t,z],2,[1,1,1]) 


I will probably try to add other functions to do arc length, line integral, green and stoke's theorem ...

PHP Code:

#cas
def metric(str,vars):
    
    if (type(str) == DOM_STRING):
        csys = LOWER(str)
    elif (type(str) == DOM_INT):
        if (str==1):
            csys = "rectangular"
        elif (str==2):
            csys = "cylinder"
        elif (str==3):
            csys = "spherical"
    
    if (csys=="rectangular"):
        return (seq(1,ii,0,len(vars)-1))
    elif (csys=="cylinder"):
        return ([1,(vars[0]),1])
    elif (csys=="spherical"):
        return ([1,(vars[0]),(vars[0])*sin(vars[1])])    

def curl2(vlist,vars,csys,uu=[]):
    """
    Example: curl2([2*z^2,5*r,-3*r^2],[r,t,z],2)
    """
    csys = metric(csys,vars)
    tmp = product(csys)
    tmp1 = (csys[0])/tmp*(diff(csys[2]*vlist[2],vars[1])-(diff(csys[1]*vlist[1],vars[2])))
    tmp2 = (csys[1])/tmp*(diff(csys[0]*vlist[0],vars[2])-(diff(csys[2]*vlist[2],vars[0]))) 
    tmp3 = (csys[2])/tmp*(diff(csys[1]*vlist[1],vars[0])-(diff(csys[0]*vlist[0],vars[1])))
    res = [tmp1,tmp2,tmp3]
    if (len(uu) == len(vars)):
        res = substituteVal(res,vars,uu)
    return simplify(res)

def dirderiv2(vlist,vars,csys,uu):
    return dot(grad2(vlist,vars,csys),mat2list(uu/ABS(uu)))

def div2(vlist,vars,csys,uu=[]):
    """
    Example: div2([3*y,5-2x,z^2-2],[x,y,z],1)
    """
    csys = metric(csys,vars)
    res = 1/(product(csys))*sum(diff(product(csys)/(csys[ii])*vlist[ii],vars[ii]),ii,0,len(vlist)-1)
    if (len(uu) == 3):
        res = substituteVal(res,vars,uu)
    return simplify(res)

def grad2(ff,vars,csys,uu=[]):
    """
    Example: grad2(r,[r,t,p],3)
    """
    csys = metric(csys,vars)
    res = seq((diff(ff,vars[ii]))/(csys[ii]),ii,0,len(csys)-1)
    if (len(uu) == len(vars)):
        res = substituteVal(res,vars,uu)
    return simplify(res)

def jacobian2(vlist,vars,uu=[]):
    """
    Example: jacobian2([r*cos(t), r*sin(t), z],[r,t,z])
    """
    jac = []
    for ii in range (0,len(vars)):
        jac = concat(jac,TRN(diff(vlist,vars[ii])))
    if (len(uu) == len(vars)):
        jac = substituteVal(jac,vars,uu)
    return simplify(jac)

def laplacian2(vlist,vars,csys,uu=[]):
    """
    Example (scalar): laplacian2(5*r*cos(t)+3*z*r^2,[r,t,z],2)
    Example (vector): laplacian2([2*r*cos(t), -4*r*sin(t), 3],[r,t,z],2)
    """
    if (type(vlist) == DOM_LIST):
        res = grad2(div2(vlist,vars,csys),vars,csys)-curl2(curl2(vlist,vars,csys),vars,csys)
    else:
        res = div2(grad2(vlist,vars,csys),vars,csys)
    if (len(uu) == len(vars)):
        res = substituteVal(res,vars,uu)
    return simplify(res)

def unitnormal2(ff,vars,csys,uu=[]):
    """
    Example: unitnormal2(z-a*x-b*y,[x,y,z],1)
    """
    gradVect = grad2(ff,vars,csys)
    if (len(uu) == len(vars)):
        gradVect = substituteVal(gradVect,vars,uu)
    if (csys == 2 or csys == "cylinder"):
        gradVect[1] = 0
    if (csys == 3 or csys == "spherical"):
        gradVect[1] = 0
        gradVect[2] = 0
    return simplify(gradVect/ABS(gradVect))

def substituteVal(res,vars,uu):
    for ii in range(0, len(uu)):
            res = subst(res,equal(vars[ii],uu[ii]))
    return res
#end