Circles, animations, and sunflowers

Capturing useful R concepts in the context of fooling around with bits of trigonometry and numerology. Nothing beats a good \(\pi\). You know what I golden mean.

Though disuse, laziness and sloth I often forget R concepts and methods, so I mess about when the weather is bad. I aim for ad hoc clarity in my code rather than generalization. This makes it a little easier for me to remember what I did. There is nothing like picking up your old code and having to reverse engineer the thing. Even more horrid for someone else. At least it aint Perl. Or SAS.

Set up the R environment

The libraries used.

require(tidyverse)    # for graphing capabilities, data manipulation
require(MASS)         # for converting decimal fractions to common fractions
require(plotly)       # for animating graphs

Circling points

The ratio of any circle’s circumference \(C\) to its diameter \(d\) is \(\pi\), an irrational number with the approximate value of \(3.1415926535897931\). \(\pi = \frac{C}{d}\). Since \(d\) is equal to twice the circle’s radius \(r\), \(\pi = \frac{C}{2r}\)

The equation can be rearranged to calculate \(C\) from \(r\) and \(r\) from \(C\). \(C = 2\pi r\), \(r = \frac{C}{2\pi}\).

The relationship is used to define angles. The circumference of the circle is a series of arcs of length \(r\). The length of the entire circle (ie. \(C\)) is \(2\pi\) times \(r\). In addition to being the radius of the circle, \(r\) measures the length of an arc on the circle. An arc length \(r\) defines an angle formed by the radius intersecting the end of the arc. This angle is the unit radian \(rad\).

An angle \(\theta\) is the ration of the length of the arc \(s\) to the radius \(r\): \(\theta = \frac{s}{r}\). As a circle grows or contracts, \(r\) and \(s\) also grow and contract. An angle defined by a radian on any sized circle is always the same value of \(\theta\).

\[\frac{s\ units}{r\ units} = \frac{s}{r} = \theta\] The units cancel out, making \(\theta\) dimensionless.

Traveling around the circumference of the circle requires \(2\pi\ rad\), so values of \(\theta\) varies from 0 to \(2\pi\) around that path. Half way around the circle \(\theta = \pi\). A quarter way around the circle \(\theta = \frac{\pi}{2}\), and an eighth around the circle \(\theta = \frac{\pi}{4}\). These values correspond to \(360^\circ\), \(180^\circ\), \(90^\circ\), and \(45^\circ\).

Polar coordinates

A polar coordinate defines a point by distance from the center of the circle (the length of \(r\)) and by the angle, \(\theta\), formed between the radius and a base location (where \(\theta = 0\). The polar coordinate is in the form \((r,\theta)\).

Relating polar to Cartesian coordinates

The trigonometric functions \(cosine\) and \(sine\) can be used to convert polar coordinates to \((x,y)\) Cartesian coordinates. \[x = r\times cos(\theta)\] \[y = r\times sin(\theta)\]

Plotting points on the circle

Here points on a unit circle are generated and graphed on a plot. We want points around the entire circle, but we want the points separated so the pattern is clear.

First we generate a series of fractions from 0 to 2 in steps of 1/8. These will be multiplied by \(pi\) to create angles from \(\frac{1}{8}\pi\) to \(2\pi\), a complete turn around the circle in steps of \(\frac{1}{8}\pi\).

These are multiplied by \(\pi\) to generate a sequence of values of \(\theta\) around the circle. The radius of the unit circle \(r = 1\) is combined with the angles to create a data frame of polar coordinates.

r <- 1
steps <- seq(from = 0, to = 2, by = 1/8)
thetas <- steps * pi
unit_circle <- data.frame("r" = 1,         # the unit circle, where radius is 1
                          "theta" = thetas
                      )
unit_circle  # a collection of polar coordinates on a circle

The data frame is expanded to add the Cartesian coordinate of each of the polar coordinates.

unit_circle$x <- unit_circle$r * cos(unit_circle$theta)
unit_circle$y <- unit_circle$r * sin(unit_circle$theta)
unit_circle

Graph (plot) the points on a circle

The defined Cartesian \((x,y)\) coordinate points of unit_circle are graphed on a plot.

ggplot(unit_circle, aes(x, y))+
  theme_minimal() +
  geom_point(color = "blue") +
  coord_fixed()

We can do some things to create labels as fractions of pi.

steps_asfractions <- fractions(steps)
unit_circle$label <- paste(as.character(steps_asfractions), "pi")
unit_circle

label_scaler <- 1.17    # this is used to move the label locations further out,
                        # i.e., puts the labels on a slightly larger cicle
unit_circle$label_x <- label_scaler * unit_circle$x
unit_circle$label_y <- label_scaler * unit_circle$y
unit_circle

Now the graph can be drawn with labels for the angles. The labels help relate pi (\(\pi\)) to angles in degrees. \(2\pi = 360^\circ\), \(\pi=180^\circ\), \(\frac{\pi}{2}=90^\circ\), \(\frac{\pi}{4}=45^\circ\), \(\frac{3\pi}{2} = 270^\circ\).

ggplot(unit_circle, aes(x, y))+
  theme_minimal() +
  geom_point(color = "blue") +
  geom_text(aes(x=label_x,
                y=label_y,
                label = label),
            size = 4
            ) +
  coord_fixed()

Changing circles

Changing position

A circle’s position is changed by adding or subtracting a constant value to the x and y coordinates.

steps <- seq(from = 0, to = 2, by = 1/32)
thetas <- steps * pi
unit_circle <- data.frame(r = 1,
                          theta = thetas,
                          x = 1 * cos(thetas),
                          y = 1 * sin(thetas)
                          )

unit_plot <- ggplot(unit_circle,aes(x, y)) +
  theme_minimal() +
  geom_point(color = "#7570b3") +
  xlim(-4, 4) +
  ylim(-4, 4) +
  coord_fixed()

unit_plot  # display the plot

The unit circle is moved on the plot by adjusting its coordinates by addition or substraction. Here the circle is moved three units right and three units left.

ggplot(unit_circle, 
       aes(x + 3, y + 3)) +
  theme_minimal() +
  geom_point(color = "#7570b3") +
  xlim(-4, 4) +
  ylim(-4, 4) +
  coord_fixed()

Show the circle at multiple new positions at once

The unit_plot graph is built on by adding new unit circles shifted to various locations on the plane.

unit_plot +
  theme_minimal() +
  xlim(-4, 4) +
  ylim(-4, 4) +
  geom_point(aes(x = x + 3, y = y),  # circle to right by adding to x coordinates
             color = "#1b9e77") +
  geom_point(aes(x = x, y = y + 3),  # circle up by adding to y coordinates
             color = "#1b9e77") +
  geom_point(aes(x = x + 3, y = y + 3),  # circle right and up
             color = "#1b9e77") +
  geom_point(aes(x = x - 3, y = y - 3),  # circle down and left
             color = "#1b9e77")

Scale for 'x' is already present. Adding another scale for 'x', which will replace the
existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the
existing scale.

Changing size

The unit circle is used as a template for circles or other size. The plot assigned to unit_plot is built on by adding new sets of x,y coordinates. The new coordinate sets are created by scaling the original unit circle coordinates, which are part of unit_plot.

unit_plot +  
  theme_minimal() +
  xlim(-4, 4) +
  ylim(-4, 4) +
  geom_point(aes(x = 2 * x, y = 2 * y),      # add a bigger circle
             color = "#d95f02") +
  geom_point(aes(x = x / 2, y = y / 2),      # add a smaller circle
             color = "#d95f02")

Scale for 'x' is already present. Adding another scale for 'x', which will replace the
existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the
existing scale.

Distort

Multiplication and division are changes to the distance from the center of the plot: changes to the radius.

varr <- (1:nrow(unit_circle))
unit_plot +
  xlim(-4, 4) +
  ylim(-4, 4) +
  theme_minimal() +
  geom_point(aes(x = 2*x, y),  # expand on the x axis
             color = "#d95f02") +
  geom_point(aes(x = x, 2*y),  # expand on the x axis
             color = "#d95f02")

Scale for 'x' is already present. Adding another scale for 'x', which will replace the
existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the
existing scale.

Concentric circles

Concentric circles share the same center. Points on concentric circles can be generating by combining a collection of radii of different length with angles around the circle.

The Cartesian cross product is an outcome of combining all members of one set with all members of another set. A cross product of a set of radii and a set of angles creates multiple polar coordinates.

#  create a collection of radii starting from the 0 origin
r = seq(from = 0, to = 5, by = 1)


#  create a set of angles from 0 to pi (equivalent to 0 to 180 degrees)
theta = seq(from = 0, to = 2, by = 1/6) * pi


#  The application of the function crossing maps each radius 
# (there are 6 radii in the set) to 
#  all rows in the angles data frame (there are 9 rows).  
# The product has 54 rows. 
#  Each row defines a polar coordinate (r,theta)
half_circles <- crossing(r, theta)
half_circles

#  Convert polar coordinates to Cartesian coordinates for plotting in x,y
half_circles$x <- half_circles$r * cos(half_circles$theta)
half_circles$y <- half_circles$r * sin(half_circles$theta) 

#  Add the equivalent angle in degrees to each theta
half_circles$degrees <- half_circles$theta * 180 / pi
half_circles

NA

Plot the points. The point is labeled with the radius of its polar coordinate. The point is colored by the equivalent degrees of theta.

ggplot(half_circles, aes(x, y)) +
  theme_minimal() +
  geom_point(aes(color = as.factor(degrees))) +
  geom_text(aes(label = r, color = as.factor(degrees)), nudge_y = 0.3) +
  geom_path(aes(color= as.factor(r))) +
  geom_path(aes(color= as.factor(degrees))) +
  coord_fixed()

Spiral

A spiral is created by increasing the radius for each point generated. Using seq_along(theta) creates a vector from 1 to the number that is the length of theta in increments of 1. We want the spriral to start at the origin, so 1 is subtracted from the generated vector, which efectively shifts the sequence to the left by 1.

theta <- seq(0, 6*pi, 1/4)  # three times around in 1/4 steps
r <- seq_along(theta)       # ever increasing radius
r <- r - 1                  # shift to start at 0

spiral <- data.frame(r = r,
                     theta = theta)
# show the polar coordinates
spiral

# derive the cartesian coordinates
spiral$x <- r * cos(theta)
spiral$y <- r * sin(theta)

# plot the spiral
ggplot(spiral, aes(x, y)) +
  theme_minimal() +
  geom_point(color = "#d95f02") +
  coord_fixed()

Animating the spiral, a bit of a tangent

Fun with spiral animation. Focus is on the creation of plot frames.

Rather than rendering the points all at once, the plot of the spiral can be animated using ggplotly() from the package plotly. This function animates the plot by drawing a series of plot frames, which are identified in the data frame.

The mechanics of generating the points of a spiral guarantees that the data frame is sequenced by radius. This makes it possible to reason about the process. Without the guarantee, some extra work (lazily and slovenly left out of this work) needs to be done. Don’t get complacent.

The spiral as a moving point

In this section, each plot frame identifies a single point. The radius of the point is used as the plot frame variable. The mechanics of the generation of the points making

Note that each frame (which is a point indexed by spiral$r) is a new plot. As the animation plays, a single point seems to move on the path of the spiral. In reality, each point is separate from the others. The preceding plots are replaced, not added to.

ggplotly(                           # the new function
  ggplot(spiral, aes(x, y)) +
  theme_minimal() +
  geom_point(aes(frame = r),        # the frame variable, which is part of spiral
             color = "#d95f02") +
  coord_fixed()
)

Warning: Ignoring unknown aesthetics: frame

The growing spiral

Rather than an apparently moving point on a spiral, I want to show a growing spiral by retaining rather than losing previous points.

In the version of plotly used, retaining points throughout the animation is not an option. We have to perform a bit of finagling.

During animation, each frame replaces the next. Retaining plotted points requires that each frame has its own point and all of the preceding points.

The process:

1. create a vector, framenumbers, that indexes each row in spiral
2. iterate over framenumbers capturing each element as framenumber
3. at each, copy the rows of spiral from 1 to framenumber, add a new column, plotframe to the copied rows and set it to framenumber, and add the copied rows as a data frame item to a growing list
4. when complete, collapse the list of data frame items into a single data frame

The new data frame contains plot frames with progressively more rows. The plot frame rows are identified in the new data frame as plotframe, which is used to step through the data to animate the plot using the function ggplotly.

framenumbers <- seq_along(spiral$r)

listofexpandedframes <- 
  lapply(framenumbers,           # elements of framenumbers are fed one at a time to the function
         function(framenumber){
           cbind(
             spiral[c(1:framenumber), ],
             plotframe = framenumber
           )
         })

growingspiral <- dplyr::bind_rows(listofexpandedframes) 

head(growingspiral,20)

ggplotly(
  ggplot(growingspiral, aes(x, y)) +
    theme_minimal() +
    geom_point(aes(frame = plotframe), 
               color = "#d95f02") +
    coord_fixed()
)

Warning: Ignoring unknown aesthetics: frame

Sunflower head patterns

See https://archive.bridgesmathart.org/2010/bridges2010-483.pdf and https://math.okstate.edu/people/segerman/talks/sunflower_spiral_talk.pdf for a discussion on the sunflower spiral. The algorithms captured here are from these sources.

Trips around the circle

A set of angles created by multiplying an angle by a set of integers results in a set of coterminal angles. In the case below, the angle is \(2\pi\) and the integers are the set from 1 to 50. \(2\pi\) is equal to \(360^\circ\). Although the fifty angles are all different angles, they all terminate with the same vertex, but with different trips around the circle to get there.

Te coterminal angles are later transformed such that they are no longer coterminal.

The radii are transformed by the square root function. The has the effect of making the points more evenly distributed. Taking the square root has more impact to large numbers than to small. For example, the square root of 100 is 10 and the square root of 1 is 1. Replacing the radii of 1 to n with the square root of 1 to n makes the radii more uniformly distributed in space, resulting in better point spacing when the points are plotted.

n <- seq(1,50)

theta <- 2*pi * n

r <- sqrt(n)

sunflower <- data.frame(n = n,
                        r = r,
                        theta = theta,
                        trips = theta / (2 * pi), 
                        x = r * cos(theta),
                        y = r * sin(theta)
                   )
sunflower

ggplot(sunflower, aes(x, round(y))) + 
  theme_minimal() +
  geom_point(color = "#d95f02") +
  ggtitle("Coterminal angles with different radius lengths")

Henry Segerman’s model

More interesting is when theta is adjusted by the golden ratio. This changes each angle theta by multiplying theta by phi, the golden ratio.
\[\theta = \theta * \phi\] \[\phi = \frac{1 + \sqrt 5}{2} \approx 1.618034\] After the transformation, the angles are multiples of the golden ratio. This have been found to have an arrangement found in nature, including the way florets and later seeds are packed into the head of a sunflower. This is the algorithm used by Henry Segerman.

The labels are the values of n in the data, which is the ordering of the points in the data.

n <- seq(1,90)

phi <- (1 + sqrt(5)) / 2.0
theta <- (2*pi * n) * phi

r <- sqrt(n)

sunflower <- data.frame(n = n,
                        r = r,
                        theta = theta,
                        x = r * cos(theta),
                        y = r * sin(theta)
                   )
sunflower

ggplot(sunflower, aes(x, y, label = n)) + 
  theme_minimal() +
  geom_point(color = "#d95f02") +
  geom_text() +
  ggtitle("Henry Segerman sunflower model") +
  coord_fixed()

Even more interesting when we greatly increase the number of angles. The result is a set of points that is a fairly good model of a sunflower head.

n <- seq(1,959)

r <- sqrt(n) 

phi <- (1 + sqrt(5)) / 2.0    # phi is the golden ratio
phi

[1] 1.618034

theta <- 2*pi * n
theta <- phi * theta          # adjust angle  by the golden ratio
                              

ssunflower <- data.frame(n = n,
                        r = r,
                        theta = theta,
                        phi_count = (theta / (2 * pi)) / phi,
                        x = r * cos(theta),
                        y = r * sin(theta),
                        style = "segerman"
                   )
ggplot(ssunflower, aes(x, y, color = style)) + 
  ggtitle("Henry Segerman sunflower model") +
  theme_minimal() +
  geom_point() +
  coord_fixed()

Helmut Vogel’s model

Proposed by Vogel in 1979. https://en.wikipedia.org/wiki/Fibonacci_number#Nature and reference to paper.
\[\theta = \frac{2\pi}{\phi^2}n,\ r = c\sqrt n\]

n <- seq(1,959)

r <- sqrt(n)

phi <- (1 + sqrt(5)) / 2.0    # phi is the golden ratio
phi

[1] 1.618034

theta <- ((2*pi) / phi^2) * n


vsunflower <- data.frame(n = n,
                         r = r,
                         theta = theta,
                         phi_count = (theta / (2 * pi)) / phi,
                         x = r * cos(theta),
                         y = r * sin(theta),
                         style = "vogel"
                   )
ggplot(vsunflower, aes(x, y, color = style)) + 
  theme_minimal() +
  geom_point() +
  ggtitle("Helmut Vogel sunflower model") +
  coord_fixed()

Comparing the two results visually

The two models are combined to plot both at the same time. This is done by the rbind function, which combines the data frame rows into a single data frame. The combined data frame is the object processed by the ggplot function.

gg <- ggplot(rbind(ssunflower, vsunflower), 
       aes(x, y, color = style)) + 
  theme_minimal() +
  geom_point() +
  ggtitle("Overlay models for comparison") +
  coord_fixed()

gg

ggplotly for richer user experience on the web

Wrapping a ggplot with ggplotly makes it possible to interact with graphs on the web even where there is no animation.

ggplotly(gg)