Sage Notebook "cycloid_animations"

lecture Theoretical Mechanics IPSP

Universität Leipzig, winter term 2020

Author: Jürgen Vollmer (2020)

lizence: Creative Commons Attribution-ShareAlike 4.0 International CC BY-SA 4.0

Sage is an OpenSource Project that implements a broach range of computer-based mathematics and computer algebra in a Jupyter notebook.

Documentation and information on the installation of Sage can be found at https://sagemath.org

A very nice introduction to working with Sage is given in the book Paul Zimmermann, u.a.: "Computational Mathematics with SageMath"

import NumPy functions

In [1]:
import numpy as np

declare variables

In [2]:
theta, phi, d = var("theta,phi,d")

position of the reflectors

is provided here as a sum of two vectors $\mathbf M$ and $\mathbf D$ that depend on $\theta$, $\phi$ and $d$

In [3]:
M =      vector([-theta, 1])
D =  d * vector([-sin(theta+phi), cos(theta+phi)])

we draw functions in the range $x=0$ to $x=$xMax
and draw resolution$=50$ points in each interval of length $1$
i.e. we evaluate the cycloid for $\theta$-values in the liste thetaValues

In [4]:
xMax       = 10.
resolution = 50.
thetaValues = -np.arange((1.+xMax)*resolution)/resolution

you can look at variables by using show:

In [5]:
show(xMax)

watch out however: the output might be long!

In [6]:
show(thetaValues)

you can comment a line by #

In [7]:
# show(thetaValues)

you get help by preceding a function name by a ?

?[function_name]

like e.g.

In [8]:
?np.arange

start interactive plot function

the slider _d provides values for $d \in [0,3]$

the slider phi_pi provides values for $\frac{\varphi}{\pi} \in [-1,1]$

the slider position provides the position of the wheel and the reflector

In [9]:
##   @interact  indicates that we are using interactive slider elements in the plot
##   this allows us to interactively explore the impact of different parameters
@interact

# the sliders and their values are introduced in square brackets
# the selected values are available in the program by the names   *_d*, *phi_pi* und *position* 
def _( _d=(0,3,0.1), phi_pi=(-1,1,0.1), position=(0,xMax,0.1) ) :

    # size of the plot window
    _xMin = -1.
    _xMax =  xMax
    if (_d<1) :
        _yMin = 0.
        _yMax = 2.
    else :
        _yMin = 1-_d
        _yMax = 1+_d
    
    # lines for trail and the height of the axis
    #  the  alpha-value  specifies the transparency of lines, ranging from 0 for  invisible  ->  1 for solid
    #  the image is stored in the variable  plt  
    plt  = plot( 0, xmin=_xMin, xmax=_xMax, ymin=_yMin, ymax=_yMax, color="black", alpha=0.1)

    #  one can add more elements by  "+="
    plt += plot( 1, xmin=_xMin, xmax=_xMax, ymin=_yMin, ymax=_yMax, color="gray" , alpha=0.1)

    # add cycloide
    #  the square brackets contain a list of  points 
    #  .subs()  substitutes numbers into the mathematical expression
    #  line()  draws the line, specified by connecting the points
    path = [(M+D).subs([d==_d, phi==phi_pi*pi, theta==_theta]) \
                 for _theta in thetaValues]
    plt += line(path, xmin=_xMin, xmax=_xMax, ymin=_yMin, ymax=_yMax )

    # add position of axis and wheel
    #  the position of the asis is  M,  evaluated for  theta=-position
    #  the minus sign is necessary in order to make the wheel go right
    #  point()  adds points  at the positions indicated by its argument
    #  circle(P, R, ...)  draws a circle with center  P  and Radius  R
    axis = M.subs([theta==-position])
    plt += point(axis, color="black", size=50)
    plt += circle(axis, 1, color="lightgray", thickness=3)
    plt += circle(axis, 1, color="gray", thickness=1, fill=True, alpha=0.1)
    
    # position of  reflector  and  spork  from  the axis to the reflector
    #  the reflektor is at the position  M+D,  evaluated for  theta=-position
    reflector = (M+D).subs([d==_d, phi==phi_pi*pi, theta==-position])
    plt += point(reflector, color="red", size=50)
    plt += line((axis, reflector), color="red" )
    
    #  display plot
    show(plt)
Widget Javascript not detected.  It may not be installed or enabled properly.

Cykloids: plot and save

In [10]:
def Cycloid( _d, phi_pi, thetaMax, *args, **kwargs ):
    """"
    this function makes a plot of a cycloid
    determined by the parameters  _d, phi_pi, thetaMax,  that are provided by the function arguments
    The optional arguments in  ars  and  kwargs  will be passed to the 'plot' function.
    This allows us to specify colors, line thickness, and other plot parameters
    """
    
    # thetaValues  -- analogous to interactive plot 
    resolution = 50.
    thetaValues = -np.arange(thetaMax*resolution)/resolution

    # size of plot-window   -- analogous to interactive plot 
    _xMin = -1.
    _xMax =  xMax
    if (_d<1) :
        _yMin = 0.
        _yMax = 2.
    else :
        _yMin = 1-_d
        _yMax = 1+_d
    
    # lines for trail and the height of the axis
    plt  = plot( 0, xmin=_xMin, xmax=_xMax, ymin=_yMin, ymax=_yMax, color="black", alpha=0.1)
    plt += plot( 1, xmin=_xMin, xmax=_xMax, ymin=_yMin, ymax=_yMax, color="gray" , alpha=0.1)

    # cycloid
    bahnkurve = [(M+D).subs([d==_d, phi==phi_pi*pi, theta==_theta]) \
                 for _theta in thetaValues]
    plt += line(bahnkurve, xmin=_xMin, xmax=_xMax, ymin=_yMin, ymax=_yMax, **kwargs )

    # return plot as result of the fuction call
    return plt
In [11]:
# initial condition for wheel and spork
phi_pi   = 0.0

# range of  theta-values
thetaMax = 3*pi+1

# Cykloid provides plot of cycloide with specified parameters
cyclo  = Cycloid( 1.0, phi_pi, thetaMax, color="red",   thickness=4 )

# ... add further cycloids with other colors and line thickness
cyclo += Cycloid( 0.8, phi_pi, thetaMax, color="blue",  thickness=3 )
cyclo += Cycloid( 0.6, phi_pi, thetaMax, color="green", thickness=2 )
cyclo += Cycloid( 0.4, phi_pi, thetaMax, color="gray",  thickness=1.5 )
cyclo += Cycloid( 0.2, phi_pi, thetaMax, color="yellow"  )

# show plot
show(cyclo)
In [12]:
# store image as  svg-file
cyclo.save_image('p04_animate-cycloids_tracks.svg', figsize=[9,4])

make and store an animation

In [13]:
def Wheel( _d, phi_pi, position ):
    """"
    make a plot of the wheel and the reflector
    parameter
     _d       distance of reflector from the fulcrum 
     phi_pi   initial position of the reflector
     position of the wheel, i.e. the angle  theta  that specifies how the the wheel has gone
    """
    
    # determine  position of axis
    axis = M.subs([theta==-position])
    # ... and plot the axis and the wheel
    plt  = point(axis, color="black", size=50)
    plt += circle(axis, 1, color="lightgray", thickness=3)
    
    # determine  position of reflector
    reflector = (M+D).subs([d==_d, phi==phi_pi*pi, theta==-position])
    # ... and plot the reflector and the spork from the fulcrum to the reflector
    plt += point(reflector, color="red", size=50)
    plt += line((axis, reflector), color="red" )

    # return plot
    return plt
In [14]:
# select paramters
_d       = 1.0
phi_pi   = 0.0
thetaMax = 3*pi

# theta-values where the position of the reflector will be evaluated
resolution = 10.
thetaValues = np.arange((1.+thetaMax)*resolution)/resolution

# plot cycloid
cyclo = Cycloid( _d, phi_pi, thetaMax, thickness=3 )

# plot wheel and reflector for different positions
# ... and make the animation
movie = animate([cyclo + Wheel( _d, phi_pi, position ) for position in thetaValues])
In [15]:
#  further information on the use of  animate()  is given when uncommenting and executing the following line:
# ?animate
In [16]:
#  ...or still more comprehensive -- as suggested in  ?animate
# ?sage.plot.animate
In [17]:
# show amimation  --  this will take a moment
movie.show()
In [18]:
#  save animation  as  animated-gif-file
movie.save("p04_animate-cycloids_d1.0.gif")

#  you might try to add the option, i.e. call
#     movie.save("p04_animate-cycloids_d1.0.gif", use_ffmpeg=True)
#  at times this results in smoother running animations

flag of values

In [19]:
def WheelFlag( phi_pi, position ):
    """
    make a plot with the wheel and the trail of several reflectors at different distance from the fulcrum
     
    parameter
       phi_pi   initial position of the reflector
       position of the wheel, i.e. the angle  theta  that specifies how far the wheel has gone
    """
    
    # determine  position  of axis and wheel
    axis = M.subs([theta==-position])
    # ... and plot the axis and the wheel
    plt  = point(axis, color="black", size=50, figsize=10, )
    plt += circle(axis, 1, axes=False, color="lightgray", thickness=3)
    plt += circle(axis, 1, axes=False, color="gray", thickness=1, fill="True", alpha=0.1)

    # determine  position  of  reflector
    reflector = (M+D).subs([d==1, phi==phi_pi*pi, theta==-position])
    # ... and plot the reflector and the spork from the fulcrum to the reflector
    plt += point(reflector, axes=False, color="red", size=50)
    plt += line((axis, reflector), axes=False, hue=sin(1) )
    
    # determine  position  of  reflectors further inside
    for _d in [0.1,0.2,..,0.9] :
        reflector = (M+D).subs([d==_d, phi==phi_pi*pi, theta==-position])
        # ... and add to the figure
        #   hue=...  selects rainbow colors for the different functions
        plt += point(reflector, axes=False, hue=sin(_d), size=40)

    # cycloids
    #   determine the tracks of the curtate trochoids
    #   from the intial position of the wheel till the present position  position
    # ... and add the tracks as lines to the plot
    plt += Cycloid( 1.0, phi_pi, position, axes=False, color="red", thickness=4 )
    for _d in [0.1,0.2,..,0.9] :
        plt += Cycloid( _d, phi_pi, position, axes=False, hue=sin(_d),  thickness=3 )

    # return resulting plot
    return plt
In [20]:
# select the resolution:  how many positions of the wheel do we consider?
#   for high resolution the wheel goes more smoothly
#   but files will be bigger and the calculation takes longer
resolution = 10.

# maximum position:  how far will the wheel go?
thetaMax = 3*pi

# thetaValues, used to draw the wheel
#   specifies the positions of the sporks in the images
thetaValues = np.arange((1.+thetaMax)*resolution)/resolution

# make animation
movie2 = animate([WheelFlag( 0, position ) for position in thetaValues])
In [21]:
# show animation
movie2.show()

#   for figures with a lot of detail the option   use_ffmpeg=True   results in smaller animations
#   when  ffmpeg  is not installed on your machine, you might try
# a.show(use_ffmpeg=True)
In [22]:
#  save animation as  animated-gif-file
movie2.save("p04_animate-cycloids_flag.gif")

# also here you may add ffmpeg-compression  as follows:
#  (however, this does not make much of a difference in the present case)
# a.save("03_4_Zyklioden_animation.gif", use_ffmpeg=True)