Animations

Generate Simple Animations

Each primitive of the plot library knows how many specifications of type "range" it has to expect. For example, a univariate function graph in 2D such as

plot::Function2d(sin(x), x = 0..2*PI):

expects one plot range for the x coordinate, whereas a bivariate function graph in 3D expects two plot ranges for the x and y coordinate:

plot::Function3d(sin(x^2 + y^2), x = 0..2, y = 0..2):

A contour plot in 2D expects 2 ranges for the x and y coordinate:

plot::Implicit2d(x^2 + y^2 - 1, x = -2..2, y = - 2..2):

A contour plot in 3D expects 3 ranges for the x, y, and z coordinate:

plot::Implicit3d(x^2 + y^2 + z^2 - 1, x = -2..2,
                 y = - 2..2, z = - 2..2):

A line in 2D does not expect any range specification:

plot::Line2d([0, 0], [1, 1]):

    Note:   Whenever a graphical primitive receives a "surplus" range specification by an equation such as a = amin..amax, the parameter a is interpreted as an "animation parameter" assuming values from amin to amax.

Thus, it is very easy indeed to create animated objects: Just pass a "surplus" range equation a = amin..amax to the generating call of the primitive. All other entries and attributes of the primitive that are symbolic expressions of the animation parameter will be animated. In the following call, both the function expression as well as the x range of the function graph depend on the animation parameter. Also, the ranges defining the width and the height of the rectangle as well as the end point of the line depend on it:

plot(
 plot::Function2d(a*sin(x), x = 0..a*PI, a = 0.5..1),
 plot::Rectangle(0..a*PI, 0..a, a = 0.5..1,
                 LineColor = RGB::Black),
 plot::Line2d([0, 0], [PI*a, a], a = 0.5 ..1,
                   LineColor = RGB::Black)
)

Additional range specifications may enter via the graphical attributes. Here is an animated arc whose radius and "angle range" depend on the animation parameter:

plot(plot::Arc2d(1 + a, [0, 0], AngleRange = 0..a*PI, a = 0..1)):

Here, the attribute AngleRange is identified by its attribute name and thus not assumed to be the specification of an animation parameter with animation range.

    Note:   Do make sure that attributes are specified by their correct names. If an incorrect attribute name is used, it may be mistaken for an animation parameter!

In the following examples, we wish to define a static semicircle, using plot::Arc2d with AngleRange = 0..PI. However, AngleRange is spelled incorrectly. A plot is created. It is an animated full circle with the animation parameter AngelRange!

plot(plot::Arc2d(1, [0, 0], AngelRange = 0..PI)):

The animation parameter may be any symbolic parameter (identifier or indexed identifier) that is different from the symbols used for the mandatory range specifications (such as the names of the independent variables in function graphs). The parameter must also be different from any of the protected names of the plot attributes.

    Note:   Animations are created object by object. The names of the animation parameters in different objects need not coincide.

In the following example, different names a, b are used for the animation parameters of the two functions:

plot(plot::Function2d(4*a*x, x = 0..1, a = 0..1),
     plot::Function2d(b*x^2, x = 0..1, b = 1..4)):

An animation parameter is a global symbolic name. It can be used as a global variable in procedures defining the graphical object. The following example features the 3D graph of a bivariate function that is defined by a procedure using the globally defined animation parameter. Further, a fill color functionmycolor is defined that changes the color in the course of the animation. It could use the animation parameter as a global parameter, just as the function f does. Alternatively, the animation parameter may be declared as an additional input parameter. Refer to the help page of FillColorFunction to find out, how many input parameters the fill color function expects and which of the input parameters is fed with the animation parameter. One finds that for plot::Function3d, the fill color function is called with the coordinates x, y, z of the points on the graph. The next input parameter (the 4th argument of mycolor) is the animation parameter:

f := (x, y) -> 4 - (x - a)^2 - (y - a)^2:
mycolor := proc(x, y, z, a) 
  local t;
  begin
     t := sqrt((x - a)^2 + (y - a)^2):
     if t < 0.1 then 
          return(RGB::Red)
     elif t < 0.4 then
          return(RGB::Orange)
     elif t < 0.7 then
          return(RGB::Green)
     else return(RGB::Blue)
     end_if;
end:
plot(plot::Function3d(f, x = -1..1, y = -1..1, a = -1..1, 
                      FillColorFunction = mycolor)):

Play Animations

When an animated plot is created in a MuPAD® notebook, the first frame of the animation appears as a static picture below the input region. To start the animation, double click on the plot. An icon for starting the animation will appear (make sure the item ‘Animation Bar' of the ‘View' menu is enabled):

One can also use the slider to animate the picture "by hand." Alternatively, the ‘Animation' menu provides an item for starting the animation.

The Number of Frames and the Time Range

By default, an animation consists of 50 different frames. The number of frames can be set to be any positive number n by specifying the attribute Frames = n. This attribute can be set in the generating call of the animated primitives, or at some higher node of the graphical tree. In the latter case, this attribute is inherited to all primitives that exist below the node. With a = amin..amax, Frames = n, the i-th frame consists of a snapshot of the primitive with

.

Increasing the number of frames does not mean that the animation runs longer; the renderer does not work with a fixed number of frames per second but processes all frames within a fixed time interval.

In the background, there is a "real time clock" used to synchronize the animation of different animated objects. An animation has a time range measured by this clock. The time range is set by the attributes TimeBegin = t0, TimeEnd = t1 or, equivalently, TimeRange = t0..t1, where t0, t1 are real numerical values representing physical times in seconds. These attribute can be set in the generating call of the animated primitives, or at some higher node of the graphical tree. In the latter case, these attributes are inherited by all primitives that exist below the node.

The absolute value of t0 is irrelevant if all animated objects share the same time range. Only the time difference t1 - t0 matters. It is (an approximation of) the physical time in seconds that the animation will last.

    Note:   The parameter range amin..amax in the specification of the animation parameter a = amin..amax together with Frames = n defines an equidistant time mesh in the time interval set by TimeBegin = t0 and TimeEnd = t1. The frame with a = amin is visible at the time t0, the frame with a = amax is visible at the time t1.

    Note:   With the default TimeBegin = 0, the value of the attribute TimeEnd gives the physical time of the animation in seconds. The default value is TimeEnd = 10, i.e., an animation using the default values will last about 10 seconds. The number of frames set by Frames = n does not influence the time interval, but changes the number of frames displayed in this time interval.

Here is a simple example:

plot(plot::Point2d([a, sin(a)], a = 0..2*PI, 
                   Frames = 100, TimeRange = 0..5)):

The point will be animated for about 5 physical seconds in which it moves along one period of the sine graph. Each frame is displayed for about 0.05 seconds. After increasing the number of frames by a factor of 2, each frame is displayed for about 0.025 seconds, making the animation somewhat smoother:

plot(plot::Point2d([a, sin(a)], a = 0..2*PI, 
                   Frames = 200, TimeRange = 0..5)):

Note that the human eye cannot distinguish between different frames if they change with a rate of more than 25 frames per second. Thus, the number of frames n set for the animation should satisfy

.

Hence, with the default time range TimeBegin = t0 = 0, TimeEnd = t1 = 10 (seconds), it does not make sense to specify Frames = n with n > 250. If a higher frame number is required to obtain a sufficient resolution of the animated object, one should increase the time for the animation by a sufficiently high value of TimeEnd.

What Can Be Animated?

We may regard a graphical primitive as a collection of plot attributes. (Indeed, also the function expression sin(x) in plot::Function2d(sin(x), x = 0..2*PI) is internally realized at the attribute Function = sin(x).) So, the question is:

   "Which attributes can be animated?"

The answer is: "Almost any attribute can be animated!" Instead of listing the attributes that allow animation, it is much easier to characterize the attributes that cannot be animated:

  • None of the canvas attributes can be animated. This includes layout parameters such as the physical size of the picture. See the help page of plot::Canvas for a complete list of all canvas attributes.

  • None of the attributes of 2D scenes and 3D scenes can be animated. This includes layout parameters, background color and style, camera positioning in 3D etc. See the help pages of plot::Scene2d and plot::Scene3d for a complete list of all scene attributes.

    Note that there are camera objects of type plot::Camera that can be placed in a 3D scene. These camera objects can be animated and allow to realize a "flight" through a 3D scene. See section Cameras in 3D for details.

  • None of the attributes of 2D coordinate systems and 3D coordinate systems can be animated. This includes viewing boxes, axes, axes ticks, and grid lines (rulings) in the background. See the help pages of plot::CoordinateSystem2d and plot::CoordinateSystem3d for a complete list of all attributes for coordinate systems.

    Although the ViewingBox attribute of a coordinate system cannot be animated, the user can still achieve animated visibility effects in 3D by clipping box objects of type plot::ClippingBox.

  • None of the attributes that are declared as "Attribute Type: inherited" on their help page can be animated. This includes size specifications such as PointSize, LineWidth etc.

  • RGB and RGBa values cannot be animated. However, it is possible to animate the coloring of lines and surfaces via user defined procedures. See the help pages LineColorFunction and FillColorFunction for details.

  • The texts of annotations such as Footer, Header, Title, legend entries, etc. cannot be animated. The position of titles, however, can be animated.

    There are special text objects plot::Text2d and plot::Text3d that allow to animate the text as well as their position.

  • Fonts cannot be animated.

  • Attributes such as DiscontinuitySearch = TRUE or FillPattern = Solid that can assume only finitely many values from a fixed discrete set cannot be animated.

Nearly all attributes not falling into one of these categories can be animated. You will find detailed information on this issue on the corresponding help pages of primitives and attributes.

Advanced Animations: The Synchronization Model

As already explained in section The Number of Frames and the Time Range, there is a "real time clock" running in the background that synchronizes the animation of different animated objects.

Each animated object has its own separate "real time life span" set by the attributes TimeBegin = t0, TimeEnd = t1 or, equivalently, TimeRange = t0..t1. The values t0, t1 represent seconds measured by the "real time clock."

In most cases, there is no need to bother about specifying the life span. If TimeBegin and TimeEnd are not specified, the default values TimeBegin = 0 and TimeEnd = 10 are used, i.e., the animation will last about 10 seconds. These values only need to be modified

  • if a shorter or longer real time period for the animation is desired, or

  • if the animation contains several animated objects, where some of the animated objects are to remain static while others change.

Here is an example for the second situation. The plot consists of 3 jumping points. For the first 5 seconds, the left point jumps up and down, while the other points remain at their initial position. Then, all points stay static for 1 second. After a total of 6 seconds, the middle point starts its animation by jumping up and down, while the left point remains static in its final position and the right points stays static in its initial position. After 9 seconds, the right point begins to move as well. The overall time span for the animation is the hull of the time ranges of all animated objects, i.e., 15 seconds in this example:

p1 := plot::Point2d(-1, sin(a), a = 0..PI, Color = RGB::Red,
                    PointSize = 5*unit::mm,
                    TimeBegin = 0, TimeEnd = 5):
p2 := plot::Point2d(0, sin(a), a = 0..PI, Color = RGB::Green,
                    PointSize = 5*unit::mm,
                    TimeBegin = 6, TimeEnd = 12):
p3 := plot::Point2d(1, sin(a), a = 0..PI, Color = RGB::Blue,
                    PointSize = 5*unit::mm,
                    TimeBegin = 9, TimeEnd = 15):
plot(p1, p2, p3, PointSize = 3.0*unit::mm, 
     YAxisVisible = FALSE):

Here, all points use the default settings VisibleBeforeBegin = TRUE and VisibleAfterEnd = TRUE which make them visible as static objects outside the time range of their animation. We set VisibleAfterEnd = FALSE for the middle point, so that it disappears after the end of its animation. With VisibleBeforeBegin = FALSE, the right point is not visible until its animation starts:

p2::VisibleAfterEnd := FALSE:
p3::VisibleBeforeBegin := FALSE:
plot(p1, p2, p3, PointSize = 3.0*unit::mm, 
     YAxisVisible = FALSE):

We summarize the synchronization model of animations:

    Note:   The total real time span of an animated plot is the physical real time given by the minimum of the TimeBegin values of all animated objects in the plot to the maximum of the TimeEnd values of all the animated objects.

  • When a plot containing animated objects is created, the real time clock is set to the minimum of the TimeBegin values of all animated objects in the plot. The real time clock is started when pushing the ‘play' button for animations in the graphical user interface.

  • Before the real time reaches the TimeBegin value t0 of an animated object, this object is static in the state corresponding to the begin of its animation. Depending on the attribute VisibleBeforeBegin, it may be visible or invisible before t0.

  • During the time from t0 to t1, the object changes from its original to its final state.

  • After the real time reaches the TimeEnd value t1, the object stays static in the state corresponding to the end of its animation. Depending on the value of the attribute VisibleAfterEnd, it may stay visible or become invisible after t1.

  • The animation of the entire plot ends with the physical time given by the maximum of the TimeEnd values of all animated objects in the plot.

Frame by Frame Animations

There are some special attributes such as VisibleAfter that are very useful to build animations from purely static objects:

    Note:   With VisibleAfter = t0, an object is invisible from the start of the animation until time t0. Then it will appear and remain visible for the rest of the animation.

    Note:   With VisibleBefore = t1, an object is visible from the start of the animation until time t1. Then it will disappear and remain invisible for the rest of the animation.

These attributes should not be combined to define a "visibility range" from t0 to t1. Use the attribute VisibleFromTo instead:

    Note:   With VisibleFromTo = t0..t1, an object is invisible from the start of the animation until time t0. Then it will appear and remain visible until time t1, when it will disappear and remain invisible for the rest of the animation.

We continue the example of the previous section in which we defined the following animated points:

p1 := plot::Point2d(-1, sin(a), a = 0..PI, Color = RGB::Red,
                    PointSize = 5*unit::mm,
                    TimeBegin = 0, TimeEnd = 5):
p2 := plot::Point2d(0, sin(a), a = 0..PI, Color = RGB::Green,
                    PointSize = 5*unit::mm,
                    TimeBegin = 6, TimeEnd = 12):
p3 := plot::Point2d(1, sin(a), a = 0..PI, Color = RGB::Blue,
                    PointSize = 5*unit::mm,
                    TimeBegin = 9, TimeEnd = 15):
p2::VisibleAfterEnd := FALSE:
p3::VisibleBeforeBegin := FALSE:

We add a further point p4 that is not animated. We make it invisible at the start of the animation via the attribute VisibleFromTo. It is made visible after 7 seconds to disappear again after 13 seconds:

p4 := plot::Point2d(0.5, 0.5, Color = RGB::Black,
                    PointSize = 5*unit::mm,
                    VisibleFromTo = 7..13):

The start of the animation is determined by p1 which bears the attribute TimeBegin = 0, the end of the animation is determined by p3 which has set TimeEnd = 15:

plot(p1, p2, p3, p4, PointSize = 3.0*unit::mm, 
     YAxisVisible = FALSE):

Although a typical MuPAD animation is generated object by object, each animated object taking care of its own animation, we can also use the attributes VisibleAfter, VisibleBefore, VisibleFromTo to build up an animation frame by frame:

    Note:   "Frame by frame animations": Choose a collection of (typically static) graphical primitives that are to be visible in the i-th frame of the animation. Set VisibleFromTo = t[i]..t[i+1] for these primitives, where t[i]..t[i+1] is the real time life span of the i-th frame (in seconds). Finally, plot all frames in a single plot command.

Here is an example. We let two points wander along the graphs of the sine and the cosine function, respectively. Each frame is to consist of a picture of two points. We use plot::Group2d to define the frame; the group forwards the attribute VisibleFromTo to all its elements:

for i from 0 to 101 do
    t[i] := i/10;
end_for:
for i from 0 to 100 do
  x := i/100*PI;
  myframe[i] := plot::Group2d(
        plot::Point2d([x, sin(x)], Color = RGB::Red),
        plot::Point2d([x, cos(x)], Color = RGB::Blue),
        VisibleFromTo = t[i]..t[i + 1]);
end_for:
plot(myframe[i] $ i = 0..100, PointSize = 5.0*unit::mm):

This "frame by frame" animation certainly needs a little bit more coding effort than the equivalent objectwise animation, where each of the points is animated:

delete i:
plot(plot::Point2d([i/100*PI, sin(i/100*PI)], i = 0..100,
                   Color = RGB::Red),
     plot::Point2d([i/100*PI, cos(i/100*PI)], i = 0..100,
                   Color = RGB::Blue),
     Frames = 101, TimeRange = 0..10,
     PointSize = 5.0*unit::mm):

There is, however, a special kind of plot where "frame by frame" animations are very useful. Note that in the present version of the graphics, new plot objects cannot be added to a scene that is already rendered. With the special "visibility" animations for static objects, however, one can easily simulate a plot that gradually builds up: Fill the frames of the animation with static objects that are visible for a limited time only. The visibility can be chosen very flexibly by the user. For example, the static objects can be made visible only for one frame (VisibleFromTo) so that the objects seem to move.

In the following example, we use VisibleAfter to fill up the plot gradually. We demonstrate the caustics generated by sunlight in a tea cup. The rim of the cup, regarded as a mirror, is given by the function , x ∈ [- 1, 1] (a semicircle). Sun rays parallel to the y-axis are reflected by the rim. After reflection at the point (x, f(x)) of the rim, a ray heads into the direction if x is positive. It heads into the direction if x is negative. Sweeping through the mirror from left to right, the incoming rays as well as the reflected rays are visualized as lines. In the animation, they become visible after the time 5 x, where x is the coordinate of the rim point at which the ray is reflected:

f := x -> -sqrt(1 - x^2):
plot(// The static rim:
     plot::Function2d(f(x), x = -1..1, Color = RGB::Black),
     // The incoming rays:
     plot::Line2d([x, 2], [x, f(x)], VisibleAfter = 5*x
                 ) $ x in [-1 + i/20 $ i = 1..39],
     // The reflected rays leaving to the right:
     plot::Line2d([x, f(x)], 
                  [1, f(x) + (1-x)*(f'(x) - 1/f'(x))/2],
                  Color = RGB::Orange, VisibleAfter = 5*x
                 ) $ x in [-1 + i/20 $ i =  1..19],
     // The reflected rays leaving to the left:
     plot::Line2d([x, f(x)], 
                  [-1, f(x) - (x+1)*(f'(x) - 1/f'(x))/2],
                  Color = RGB::Orange, VisibleAfter = 5*x
                 ) $ x in [-1 + i/20 $ i = 21..39],
     ViewingBox = [-1..1, -1..1]):

Compare the spherical mirror with a parabolic mirror that has a true focal point:

f := x -> -1 + x^2:
plot(// The static rim:
     plot::Function2d(f(x), x = -1..1, Color = RGB::Black),
     // The incoming rays:
     plot::Line2d([x, 2], [x, f(x)], VisibleAfter = 5*x
                 ) $ x in [-1 + i/20 $ i = 1..39],
     // The reflected rays leaving to the right:
     plot::Line2d([x, f(x)], 
                  [1, f(x) + (1-x)*(f'(x) - 1/f'(x))/2],
                  Color = RGB::Orange, VisibleAfter = 5*x
                 ) $ x in [-1 + i/20 $ i =  1..19],
     // The reflected rays leaving to the left:
     plot::Line2d([x, f(x)], 
                  [-1, f(x) - (x+1)*(f'(x) - 1/f'(x))/2],
                  Color = RGB::Orange, VisibleAfter = 5*x
                 ) $ x in [-1 + i/20 $ i = 21..39],
     ViewingBox = [-1..1, -1..1]):

Examples

Example 1

We build a 2D animation that displays a function f(x) together with the integral . The area between the graph of f and the x-axis is displayed as an animated hatch object. The current value of F(x) is displayed by an animated text:

DIGITS := 2:
// the function:
f := x -> cos(x^2):
// the anti-derivative:
F := x -> numeric::int(f(y), y = 0..x):
// the graph of f(x):
g := plot::Function2d(f(x), x = 0..6, Color = RGB::Blue):
// the graph of F(x):
G := plot::Function2d(F(x), x = 0..6, Color = RGB::Black):
// a point moving along the graph of F(x):
p := plot::Point2d([a, F(a)], a = 0..6, Color = RGB::Black):
// hatched region between the origin and the moving point p:
h := plot::Hatch(g, 0, 0 ..a, a = 0..6, Color = RGB::Red):
// the right border line of the hatched region:
l := plot::Line2d([a, 0], [a, f(a)], a = 0..6, 
                  Color = RGB::Red):
// a dashed vertical line from f to F:
L1 := plot::Line2d([a, f(a)], [a, F(a)], a = 0..6, 
                  Color = RGB::Black, LineStyle = Dashed):
// a dashed horizontal line from the y axis to F:
L2 := plot::Line2d([-0.1, F(a)], [a, F(a)], a = 0..6, 
                  Color = RGB::Black, LineStyle = Dashed):
// the current value of F at the moving point p:
t := plot::Text2d(a -> F(a), [-0.2, F(a)], a = 0..6,
                  HorizontalAlignment = Right):
plot(g, G, p, h, l, L1, L2, t, 
     YTicksNumber = None, YTicksAt = [-1, 1]):
delete DIGITS:

Example 2

We build two 3D animations. The first starts with a rectangular strip that is deformed to an annulus in the x, y plane:

c := a -> 1/2 *(1 - 1/sin(PI/2*a)):
mycolor := (u, v, x, y, z) -> [(u - 0.8)/0.4, 0, (1.2 - u)/0.4]:
rectangle2annulus := plot::Surface(
   [c(a) + (u - c(a))*cos(PI*v), (u - c(a))*sin(PI*v), 0],
   u = 0.8..1.2, v = -a..a, a = 1/10^10..1,
   FillColorFunction = mycolor, Mesh = [3, 40], Frames = 40):
plot(rectangle2annulus, Axes = None,
     CameraDirection = [-11, -3, 3]):

The second animation twists the annulus to become a Moebius strip:

annulus2moebius := plot::Surface(
   [((u - 1)*cos(a*v*PI/2) + 1)*cos(PI*v),
    ((u - 1)*cos(a*v*PI/2) + 1)*sin(PI*v),
     (u - 1)*sin(a*v*PI/2)],
   u = 0.8..1.2, v = -1..1, a = 0..1,
   FillColorFunction = mycolor, Mesh = [3, 40], Frames = 20):
plot(annulus2moebius, Axes = None,
     CameraDirection = [-11, -3, 3]):

Note that the final frame of the first animation coincides with the first frame of the second animation. To join the two separate animations, we can set appropriate visibility ranges and plot them together. After 5 seconds, the first animation object vanishes and the second takes over:

rectangle2annulus::VisibleFromTo := 0..5:
annulus2moebius::VisibleFromTo := 5..7:
plot(rectangle2annulus, annulus2moebius, Axes = None,
     CameraDirection = [-11, -3, 3]):

Example 3

In this example, we consider the planar celestial 3 body problem. We solve the system of differential equations

,

,

,

,

,

,

which is nothing but the equations of motions for two planets with masses m1, m2 at positions (x1, y1), (x2, y2) revolving in the x, y plane around a sun of mass ms positioned at (xs, ys). We specify the mass ratios: The first planet is a giant with a mass m1 that is 4% of the sun's mass. The second planet is much smaller:

ms := 1: m1 := 0.04: m2 := 0.0001:

As we will see, the motion of the giant is nearly undisturbed by the small planet. The small one, however, is heavily disturbed by the giant and, finally, kicked out of the system after a near collision.

We solve the ODEs via the MuPAD numerical ODE solve numeric::odesolve2 that provides a solution vector

.

The initial conditions are chosen such that the total momentum vanishes, i.e., the total center of mass stays put (at the origin):

Y := numeric::odesolve2(numeric::ode2vectorfield(
 {xs''(t) = 
   -m1*(xs(t)-x1(t))/sqrt((xs(t)-x1(t))^2 + (ys(t)-y1(t))^2)^3
   -m2*(xs(t)-x2(t))/sqrt((xs(t)-x2(t))^2 + (ys(t)-y2(t))^2)^3,
  ys''(t) = 
   -m1*(ys(t)-y1(t))/sqrt((xs(t)-x1(t))^2 + (ys(t)-y1(t))^2)^3
   -m2*(ys(t)-y2(t))/sqrt((xs(t)-x2(t))^2 + (ys(t)-y2(t))^2)^3,
  x1''(t) = 
   -ms*(x1(t)-xs(t))/sqrt((x1(t)-xs(t))^2 + (y1(t)-ys(t))^2)^3
   -m2*(x1(t)-x2(t))/sqrt((x1(t)-x2(t))^2 + (y1(t)-y2(t))^2)^3,
  y1''(t) =
   -ms*(y1(t)-ys(t))/sqrt((x1(t)-xs(t))^2 + (y1(t)-ys(t))^2)^3
   -m2*(y1(t)-y2(t))/sqrt((x1(t)-x2(t))^2 + (y1(t)-y2(t))^2)^3,
  x2''(t) = 
   -ms*(x2(t)-xs(t))/sqrt((x2(t)-xs(t))^2 + (y2(t)-ys(t))^2)^3
   -m1*(x2(t)-x1(t))/sqrt((x2(t)-x1(t))^2 + (y2(t)-y1(t))^2)^3,
  y2''(t) = 
   -ms*(y2(t)-ys(t))/sqrt((x2(t)-xs(t))^2 + (y2(t)-ys(t))^2)^3
   -m1*(y2(t)-y1(t))/sqrt((x2(t)-x1(t))^2 + (y2(t)-y1(t))^2)^3,
  xs(0)  = -m1   ,   x1(0)  = ms,     x2(0)  =  0,
  ys(0)  = 0.7*m2,   y1(0)  = 0,      y2(0)  = -0.7*ms,
  xs'(0) = -1.01*m2, x1'(0) = 0,      x2'(0) =  1.01*ms,
  ys'(0) = -0.9*m1,  y1'(0) = 0.9*ms, y2'(0) =  0},
 [xs(t), xs'(t), ys(t), ys'(t),
  x1(t), x1'(t), y1(t), y1'(t), 
  x2(t), x2'(t), y2(t), y2'(t)]
)):

The positions [xs(t), ys(t)] = [Y(t)[1], Y(t)[3]], [x1(t), y1(t)] = [Y(t)[5], Y(t)[7]], [x2(t), y2(t)] = [Y(t)[9], Y(t)[11]] are computed on an equidistant time mesh with dt = 0.05. The animation is built up "frame by frame" by defining static points with suitable values of VisibleFromTo and static line segments with suitable values of VisibleAfter.

Setting VisibleFromTo = t..t + 0.99*dt, each solution point is visible only for a short time (the factor 0.99 makes sure that not two points can be visible simultaneously on each orbit). The orbits of the points are realized as line segments from the positions at time t - dt to the positions at time t. The line segments become visible at time t and stay visible for the rest of the animation (VisibleAfter = t), thus leaving a "trail" of the moving points. We obtain the following graphical solution (the computation takes about two minutes on a 1 GHz computer):

dt := 0.05: imax := 516:
plot(// The sun:
     plot::Point2d(Y(t)[1], Y(t)[3], Color = RGB::Orange,
                   VisibleFromTo = t..t + 0.99*dt,
                   PointSize = 4*unit::mm
                  ) $  t in [i*dt $ i = 0..imax],
     // The giant planet:
     plot::Point2d(Y(t)[5], Y(t)[7], Color = RGB::Red,
                   VisibleFromTo = t..t + 0.99*dt,
                   PointSize = 3*unit::mm 
                  ) $  t in [i*dt $ i = 0..imax],
     // The orbit of the giant planet:
     plot::Line2d([Y(t - dt)[5], Y(t - dt)[7]], 
                  [Y(t)[5], Y(t)[7]], Color = RGB::Red,
                   VisibleAfter = t
                 ) $  t in [i*dt $ i = 1..imax],
     // The small planet:
     plot::Point2d(Y(t)[9], Y(t)[11], Color = RGB::Blue,
                   VisibleFromTo = t..t + 0.99*dt,
                   PointSize = 2*unit::mm
                  ) $  t in [i*dt $ i = 0..imax],
     // The orbit of the small planet:
     plot::Line2d([Y(t - dt)[9], Y(t - dt)[11]], 
                  [Y(t)[9], Y(t)[11]], Color = RGB::Blue,
                  VisibleAfter = t
                 ) $  t in [i*dt $ i = 1..imax]
):

Was this topic helpful?