Search code examples
pythonmathmatplotlibpolar-coordinates

How to draw a curve in Cartesian coordinates in a semicircle tube in polar coordinates?


import numpy as np
import matplotlib.pylab as plt

def tube():
    theta = np.linspace(0, np.pi/2, 30)

    x = np.cos(theta)
    y = np.sin(theta)
    z = x*0.8
    w = y*0.8

    plt.plot(z,w)
    plt.plot(x,y)
    plt.axis("equal")
    plt.show()

print plt.figure(1);tube()

enter image description here

def euler():
    A, B, a = 40, 10, 2

    t  = 10  # time
    dt = 1e-3 # interval

    nbpt = int(t/dt)

    n = 1
    s = 1. # sign of the derivative, initially chosen
    y = [0]*nbpt # result

    while n < nbpt:
        yp2 = B - A*y[n-1]**a
        if yp2 < 0:
            s = -s
            n -= 1 # recalculating the previous value
        else:
            y[n] = y[n-1] + dt*s*np.sqrt(yp2)
            n += 1

    plt.plot(np.linspace(0,t,nbpt),y)
    plt.show()

print plt.figure(2);euler()

enter image description here

I want to draw the curve made with euler() in the tube made with tube(). I guess I must go from cartesian coordinates to polar coordinates but is there anyway to make the process easier with Python ?


Solution

  • There are many ways to do it since the question does not fully pin down what transformation you are looking for. However, assuming any transformation will do so long as the resultant curve oscillates between the boundary lines of the tube, you could use:

    def polarmap(x, y):
        # normalize x and y from 0 to 1
        x = (x-x.min())/(x.max()-x.min())
        y = (y-y.min())/(y.max()-y.min())
    
        # make theta go from 0 to pi/2
        theta = np.pi*x/2
    
        # make r go from 0.8 to 1.0 (the min and max tube radius)
        r = 0.2*y + 0.8
    
        # convert polar to cartesian
        x = r*np.cos(theta)
        y = r*np.sin(theta)
        plt.plot(x, y)
    

    For example,

    import numpy as np
    import matplotlib.pylab as plt
    
    
    def tube():
        theta = np.linspace(0, np.pi/2, 30)
    
        x = np.cos(theta)
        y = np.sin(theta)
        z = x*0.8
        w = y*0.8
    
        plt.plot(z,w)
        plt.plot(x,y)
    
    def euler():
        A, B, a = 40, 10, 2
    
        t  = 10  # time
        dt = 1e-3 # interval
    
        nbpt = int(t/dt)
    
        n = 1
        s = 1. # sign of the derivative, initially chosen
        y = [0]*nbpt # result
    
        while n < nbpt:
            yp2 = B - A*y[n-1]**a
            if yp2 < 0:
                s = -s
                n -= 1 # recalculating the previous value
            else:
                y[n] = y[n-1] + dt*s*np.sqrt(yp2)
                n += 1
    
        x = np.linspace(0,t,nbpt)
        y = np.array(y)
        return x, y
    
    def polarmap(x, y):
        # normalize x and y from 0 to 1
        x = (x-x.min())/(x.max()-x.min())
        y = (y-y.min())/(y.max()-y.min())
    
        # make theta go from 0 to pi/2
        theta = np.pi*x/2
    
        # make r go from 0.8 to 1.0 (the min and max tube radius)
        r = 0.2*y + 0.8
    
        # convert polar to cartesian
        x = r*np.cos(theta)
        y = r*np.sin(theta)
        plt.plot(x, y)
    
    fig, ax = plt.subplots()
    tube()
    x, y = euler()
    polarmap(x, y)
    plt.axis("equal")
    plt.show()
    

    which yields enter image description here


    Notice that in polarmap the first step was to normalize both x and y so that they both ranged from 0 to 1. You can think of them as parameters on equal footing. If you swap the two parameters before passing them to polarmap, e.g.:

    x, y = euler()
    x, y = y, x    # swap x and y
    polarmap(x, y)
    

    then you get

    enter image description here