**MuPAD® notebooks will be removed in a future release. Use MATLAB® live scripts instead.**

**MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.**

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]):

Whenever a graphical primitive receives a “surplus”
range specification by an equation such as * a =
amin..amax*, the parameter

`a`

`amin`

`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.

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.

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 function`

`mycolor`

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)):

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.

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.

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`

.

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`

.

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.

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:

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.

There are some special attributes such as `VisibleAfter`

that
are very useful to build animations from purely static objects:

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.

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:

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:

*“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]):

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:

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]):

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 *m*_{1}, *m*_{2} at
positions (*x*_{1}, *y*_{1}), (*x*_{2}, *y*_{2}) revolving
in the *x*, *y* plane
around a sun of mass *m*_{s} positioned
at (*x*_{s}, *y*_{s}).
We specify the mass ratios: The first planet is a giant with a mass *m*_{1} 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 [*x*_{s}(*t*), *y*_{s}(*t*)] = ```
[Y(t)[1],
Y(t)[3]]
```

, [*x*_{1}(*t*), *y*_{1}(*t*)] = ```
[Y(t)[5],
Y(t)[7]]
```

, [*x*_{2}(*t*), *y*_{2}(*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] ):