# Gummy Worms

Just keep swimming (Dory in Finding Nemo)

Inspired by this post, I decided to create gummy worms like this:

Or these:

When I was young I used to eat them.

Do you want to try? This is the code:

```library(rgl)
library(RColorBrewer)
t=seq(1, 6, by=.04)
f = function(a, b, c, d, e, f, t) exp(-a*t)*sin(t*b+c)+exp(-d*t)*sin(t*e+f)
v1=runif(6,0,1e-02)
v2=runif(6, 2, 3)
v3=runif(6,-pi/2,pi/2)
open3d()
spheres3d(x=f(v1[1], v2[1], v3[1], v1[4], v2[4], v3[4], t),
y=f(v1[2], v2[2], v3[2], v1[5], v2[5], v3[5], t),
z=f(v1[3], v2[3], v3[3], v1[6], v2[6], v3[6], t),
```

# The Coaster Maker by Shiny

The word you invented is well formed and could be used in the Italian language (The Accademia della Crusca regarding to the word “Petaloso”, recently invented by an eight-year-old boy)

Are you tired of your old coasters? Do you like to make things by your own? Do you have a PC and a printer at home? If you answered yes to all these questions, just follow these simple instructions:

• Install R and RStudio in your PC
• Open RStudio and create a new Shiny Web App multiple file (ui.R/server.R)
• Substitute sample code of each file by the code below
• Press Run App
• Press buttom Get your coaster! until you obtain a image you like
• Print the image
• Cut out the image
• Place on the coaster your favorite drinking

These are some examples:

This is the code of ui.R

```#
# This is the user-interface definition of a Shiny web application. You can
# run the application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
library(shiny)
shinyUI(fluidPage(
titlePanel("The coaster maker"),
sidebarLayout(
sidebarPanel(
#helpText(),

# adding the new div tag to the sidebar
tags\$p("This coasters are generated by hypocycloid curves.The curve is formed by the locus of a point,
attached to a circle, that rolls on the inside of another circle.
In the curve's equation the first part denotes the relative position between the two circles,
the second part denotes the rotation of the rolling circle.")),
HTML("

),
),
mainPanel(
plotOutput("HarmPlot")
)
)
))
```

This is the code of server.R

```# This is the server logic of a Shiny web application. You can run the
# application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
library(shiny)
library(ggplot2)
CreateDS = function ()
{
t=seq(-31*pi, 31*pi, 0.002)
a=sample(seq(from=1/31, to=29/31, by=2/31), 1)
b=runif(1, min = 1, max = 3)
data.frame(x=(1-a)*cos(a*t)+a*b*cos((1-a)*t), y=(1-a)*sin(a*t)-a*b*sin((1-a)*t))
}
shinyServer(function(input, output) {
dat<-reactive({if (input\$rerun) dat=CreateDS() else dat=CreateDS()})
output\$HarmPlot<-renderPlot({
ggplot(dat())+
geom_point(data=data.frame(x=0,y=0), aes(x,y), color=rgb(rbeta(1, .5, .5), rbeta(1, .5, .5), rbeta(1, .5, .5)) , shape=19, fill="yellow", size=220)+
geom_polygon(aes(x, y), fill=rgb(rbeta(1, 2, 2), rbeta(1, 2, 2), rbeta(1, 2, 2))) +
theme(legend.position="none",
panel.background = element_rect(fill="white"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
}, height = 500, width = 500)
})
```

# Going Bananas With Hilbert

It seemed that everything is in ruins, and that all the basic mathematical concepts have lost their meaning (Naum Vilenkin, Russian mathematician, regarding to the discovery of Peano’s curve)

Giuseppe Peano found in 1890 a way to draw a curve in the plane that filled the entire space: just a simple line covering completely a two dimensional plane. Its discovery meant a big earthquake in the traditional structure of mathematics. Peano’s curve was the first but not the last: one of these space-filling curves was discovered by Hilbert and takes his name. It is really beautiful:

Hilbert’s curve can be created iteratively. These are the first six iterations of its construction:

As you will see below, R code to create Hilbert’s curve is extremely easy. It is also very easy to play with the curve, altering the order in which points are sorted. Changing the initial `matrix(1)` by some other number, resulting curves are quite appealing:

Let’s go futher. Changing `ggplot` geometry from `geom_path` to `geom_polygon` generate some crazy pseudo-tessellations:

And what if you change the matrix exponent?

And what if you apply polar coordinates?

We started with a simple line and with some small changes we have created fantastical images. And all these things only using black and white. Do you want to add some colors? Try with the following code (if you draw something interesting, please let me know):

```library(reshape2)
library(dplyr)
library(ggplot2)
opt=theme(legend.position="none",
panel.background = element_rect(fill="white"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
hilbert = function(m,n,r) {
for (i in 1:n)
{
tmp=cbind(t(m), m+nrow(m)^2)
m=rbind(tmp, (2*nrow(m))^r-tmp[nrow(m):1,]+1)
}
melt(m) %>% plyr::rename(c("Var1" = "x", "Var2" = "y", "value"="order")) %>% arrange(order)}
# Original
ggplot(hilbert(m=matrix(1), n=1, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=2, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=3, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=4, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=6, r=2), aes(x, y)) + geom_path()+ opt
# Changing order
ggplot(hilbert(m=matrix(.5), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(0), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(tan(1)), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(3), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(-1), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(log(.1)), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(-15), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(-0.001), n=5, r=2), aes(x, y)) + geom_path()+ opt
# Polygons
ggplot(hilbert(m=matrix(log(1)), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(.5), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(tan(1)), n=5, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-15), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-25), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(0), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(1000000), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-1), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-.00001), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
# Changing exponent
gplot(hilbert(m=matrix(log(1)), n=4, r=-1), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(.5), n=4, r=-2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(tan(1)), n=4, r=6), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-15), n=3, r=sin(2)), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-25), n=4, r=-.0001), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(0), n=4, r=200), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(1000000), n=3, r=.5), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-1), n=4, r=sqrt(2)), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-.00001), n=4, r=52), aes(x, y)) + geom_polygon()+ opt
# Polar coordinates
ggplot(hilbert(m=matrix(1), n=4, r=2), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(-1), n=5, r=2), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(.1), n=2, r=.5), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(1000000), n=2, r=.1), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(.25), n=3, r=3), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(tan(1)), n=5, r=1), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(1), n=4, r=1), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(log(1)), n=3, r=sin(2)), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(-.0001), n=4, r=25), aes(x, y)) + geom_polygon()+ coord_polar()+opt
```

# The 2D-Harmonograph In Shiny

If you wish to make an apple pie from scratch, you must first invent the universe (Carl Sagan)

I like Shiny and I can’t stop converting into appsÂ some of my previous experiments:Â today is the turn of the harmonograph. This is a simple application since you only can press a button to generate a random harmonograph-simulated curve. I love how easy is to create a nice interactive appÂ to play with from an existing code. The only trick in this case in to add a `rerun` button in the UI side andÂ transfer the interaction to the server side using a simple `if`. Very easy. This is a screenshot of the application:

Press the button and you will get a new drawing. Most of them are niceÂ scrawls and from time to time you will obtainÂ beautiful shapelyÂ curves.

And no more candyÂ words: It is time to complain. I sayÂ toÂ RStudio with all due respect, you are very cruel.Â You let me to deploy my previous app to your server but you suspended it almost immediatelyÂ for fifteen daysÂ due to “exceeded usage hours”. My only option is paying at least \$440 per year toÂ upgrade my current plan. I tried the ambrosia for an extremely short time. RStudio: Why don’t you launch a cheaper account? Why don’t you launch a free account with just one perpetual alive app at a time? Why don’t you increase the usage hours threshold? I can help you to calculate the return on investmentÂ of these scenarios.

Or, Why don’t you make me a gift for my next birthday? I promise to upload a new app per month to promote your stunning tool. Think about it and please let me know your conclusions.

Meanwhile I will run the app privately. This is the code to do it:

UI.R

```# This is the user-interface definition of a Shiny web application.
# You can find out more about building applications with Shiny here:
#
# http://www.rstudio.com/shiny/

library(shiny)

shinyUI(fluidPage(
titlePanel("Mathematical Beauties: The Harmonograph"),
sidebarLayout(
sidebarPanel(
#helpText(),

# adding the new div tag to the sidebar
tags\$p("A harmonograph is a mechanical apparatus that employs pendulums to create a
geometric image. The drawings created typically are Lissajous curves, or
related drawings of greater complexity. The devices, which began to appear
in the mid-19th century and peaked in popularity in the 1890s, cannot be
conclusively attributed to a single person, although Hugh Blackburn, a professor
of mathematics at the University of Glasgow, is commonly believed to be the official
inventor. A simple, so-called \"lateral\" harmonograph uses two pendulums to control the movement
of a pen relative to a drawing surface. One pendulum moves the pen back and forth along
one axis and the other pendulum moves the drawing surface back and forth along a
perpendicular axis. By varying the frequency and phase of the pendulums relative to
one another, different patterns are created. Even a simple harmonograph as described
can create ellipses, spirals, figure eights and other Lissajous figures (Source: Wikipedia)")),
HTML("<p>Click <a href=\"http://paulbourke.net/geometry/harmonograph/harmonograph3.html\">here</a> to see an image of a real harmonograph</p>")
),
actionButton('rerun','Launch the harmonograph!')
),
mainPanel(
plotOutput("HarmPlot")
)
)
))
```

server.R

```# This is the server logic for a Shiny web application.
# You can find out more about building applications with Shiny here:
#
# http://www.rstudio.com/shiny/
#

library(shiny)

CreateDS = function ()
{

f=jitter(sample(c(2,3),4, replace = TRUE))
d=runif(4,0,1e-02)
p=runif(4,0,pi)
xt = function(t) exp(-d[1]*t)*sin(t*f[1]+p[1])+exp(-d[2]*t)*sin(t*f[2]+p[2])
yt = function(t) exp(-d[3]*t)*sin(t*f[3]+p[3])+exp(-d[4]*t)*sin(t*f[4]+p[4])
t=seq(1, 200, by=.0005)
data.frame(t=t, x=xt(t), y=yt(t))}

shinyServer(function(input, output) {
dat<-reactive({if (input\$rerun) dat=CreateDS() else dat=CreateDS()})
output\$HarmPlot<-renderPlot({
plot(rnorm(1000),xlim =c(-2,2), ylim =c(-2,2), type="n")
with(dat(), plot(x,y, type="l", xlim =c(-2,2), ylim =c(-2,2), xlab = "", ylab = "", xaxt='n', yaxt='n', col="gray10", bty="n"))
}, height = 650, width = 650)
})
```

# 3D-Harmonographs In Motion

I would be delighted to co write a post (Andrew Wyer)

One of the best things about writingÂ a blog is that occasionally you get to know very interesting people.Â Last October 13th I published thisÂ post about the harmonograph, a machine driven by pendulums which creates amazing curves. Two days laterÂ someone called Andrew Wyer made this comment about the post:

Hi, I was fascinated by the harmonographs â€“ I remember seeing similar things done on paper on kids tv in the seventies. I took your code and extended it into 3d so I could experiment with the rgl package. I created some beautiful figures (which I would attach if this would let me).Â In lieu of that here is the code:

I ran his code and I was instantly fascinated: resulting curves were really beautiful. I suggested that weÂ co-write a post and he was delighted with the idea. He proposed to me the following improvement of his own code:

I will try to create an animated gif of one figure

Such a good idea! And no sooner said than done: AndrewÂ rewrote his own code to create stunning animated images of 3D-Harmonograph curves like these:

• Andrew creates 3D curves by adding a third oscillation z generated in the same way as x and y and adds a little colour by setting the colour of each point to a colour in the RGB scale related to its point in 3d space
• Function `spheres3d `to produce an interactive plot that you can drag around to view from different angles; function `spin3d `will rotate the figure around the z axis and at 5 rpm in this case and function `movie3d `renders each frame in a temporary png file and then calls ImageMagick to stitch them into an animated gif file. It is necessary to install ImageMagick separately to create the movie.
• Giving it a duration of 12 seconds at 5 rpm is one rotation which at 12 frames per second results in 144 individual png files but these (by default) are temporary and deleted when the gif is produced.

Although I don’t know Andrew personally, I know he is a good partner to work with. Thanks a lot for sharing this work of art with me and allowingÂ me to share it in Ripples as well.

Here you have the code. I like to imagine these curves as orbits of unexplored planets in a galaxy far, far away …

```library(rgl)
library(scatterplot3d)
#Extending the harmonograph into 3d
#Antonio's functions creating the oscillations
xt = function(t) exp(-d1*t)*sin(t*f1+p1)+exp(-d2*t)*sin(t*f2+p2)
yt = function(t) exp(-d3*t)*sin(t*f3+p3)+exp(-d4*t)*sin(t*f4+p4)
#Plus one more
zt = function(t) exp(-d5*t)*sin(t*f5+p5)+exp(-d6*t)*sin(t*f6+p6)
#Sequence to plot over
t=seq(1, 100, by=.001)
#generate some random inputs
f1=jitter(sample(c(2,3),1))
f2=jitter(sample(c(2,3),1))
f3=jitter(sample(c(2,3),1))
f4=jitter(sample(c(2,3),1))
f5=jitter(sample(c(2,3),1))
f6=jitter(sample(c(2,3),1))
d1=runif(1,0,1e-02)
d2=runif(1,0,1e-02)
d3=runif(1,0,1e-02)
d4=runif(1,0,1e-02)
d5=runif(1,0,1e-02)
d6=runif(1,0,1e-02)
p1=runif(1,0,pi)
p2=runif(1,0,pi)
p3=runif(1,0,pi)
p4=runif(1,0,pi)
p5=runif(1,0,pi)
p6=runif(1,0,pi)
#and turn them into oscillations
x = xt(t)
y = yt(t)
z = zt(t)
#create values for colours normalised and related to x,y,z coordinates
cr = abs(z)/max(abs(z))
cg = abs(x)/max(abs(x))
cb = abs(y)/max(abs(y))
dat=data.frame(t, x, y, z, cr, cg ,cb)
#plot the black and white version
with(dat, scatterplot3d(x,y,z, pch=16,cex.symbols=0.25, axis=FALSE ))
with(dat, scatterplot3d(x,y,z, pch=16, color=rgb(cr,cg,cb),cex.symbols=0.25, axis=FALSE ))
#Set the stage for 3d plots
# clear scene:
clear3d("all")
# white background
bg3d(color="white")
#lights...camera...
light3d()
#action
# draw shperes in an rgl window
#create animated gif (call to ImageMagic is automatic)
movie3d( spin3d(axis=c(0,0,1),rpm=5),fps=12, duration=12 )
#2d plots to give plan and elevation shots
plot(x,y,col=rgb(cr,cg,cb),cex=.05)
plot(y,z,col=rgb(cr,cg,cb),cex=.05)
plot(x,z,col=rgb(cr,cg,cb),cex=.05)
```

# Beautiful Curves: The Harmonograph

Each of us has their own mappa mundi (Gala, my indispensable friend)

The harmonograph is a mechanism which, by means of several pendulums, draws trajectories that can be analyzed not only from a mathematical point of view but also from an artistic one. In its double pendulum version, one pendulum moves a pencil and the other one moves a platform with a piece of paper on it. You can see an example here. The harmonograph is easy to use: you only have to put pendulums into motion and wait for them to stop. The result are amazing undulating drawings like this one:

First harmonographs were built in 1857 by Scottish mathematician Hugh Blackburn, based on the previous work of French mathematician Jean Antoine Lissajous. There is not an unique way to describe mathematically the motion of the pencil. I have implemented the one I found inÂ this sensational blog, where motion in both x and y axis depending on time is defined by:

```x(t)=e-d1tsin(tf1+p1)+e-d2tsin(tf2+p2) y(t)=e-d3tsin(tf3+p3)+e-d4tsin(tf4+p4)```

I initializeÂ parameters randomly so every time you run the script, you obtainÂ a different output. Here isÂ a mosaic with someÂ of mine:

This is the code to simulate the harmonograph (no extra package is required). If you create some niceÂ work of art, I will be very happy to admire it (you can find my email here):

```f1=jitter(sample(c(2,3),1));f2=jitter(sample(c(2,3),1));f3=jitter(sample(c(2,3),1));f4=jitter(sample(c(2,3),1))
d1=runif(1,0,1e-02);d2=runif(1,0,1e-02);d3=runif(1,0,1e-02);d4=runif(1,0,1e-02)
p1=runif(1,0,pi);p2=runif(1,0,pi);p3=runif(1,0,pi);p4=runif(1,0,pi)
xt = function(t) exp(-d1*t)*sin(t*f1+p1)+exp(-d2*t)*sin(t*f2+p2)
yt = function(t) exp(-d3*t)*sin(t*f3+p3)+exp(-d4*t)*sin(t*f4+p4)
t=seq(1, 100, by=.001)
dat=data.frame(t=t, x=xt(t), y=yt(t))
with(dat, plot(x,y, type="l", xlim =c(-2,2), ylim =c(-2,2), xlab = "", ylab = "", xaxt='n', yaxt='n'))
```