Ikeda map (fractal)

2018-07-09

Generating an Ikeda Map

Today in my ongoing quest to generate fractals in Rstats using every available avenue: generating an Ikeda map.

Methods for generating fractals so far:

Ikeda map

An Ikeda map is a model of light going acound a nonlinear optical resonator. To be honest, the wiki page and the original paper by Ikeda are impenetrable to me!

However, creating the map is relatively straightforward.

  1. Start at a random point
  2. Update the position of the point
  3. Continue updating the point’s position for N steps
  4. Draw the path joining all the positions of this point
  5. Repeat steps 1-4
library(tidyverse)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#' Computing a single path from the given starting point
#'
#' @param u ikeda parameter
#' @param x,y starting point
#' @param N is the number of iterations
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
compute_ikeda_trajectory <- function(u=0.918, x1, y1, N) {
  x    <- y <- double(N)
  x[1] <- x1
  y[1] <- y1

  for (i in 2:N) {
    t    <- 0.4 - 6/(1 + x[i-1]*x[i-1] + y[i-1]*y[i-1])
    x[i] <- 1 + u * (x[i-1] * cos(t) - y[i-1] * sin(t))
    y[i] <-     u * (x[i-1] * sin(t) + y[i-1] * cos(t))
  }

  data_frame(x=x, y=y)
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Generate P random starting points
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P <- 2000  # Number of starting  points
N <- 1000  # Number of iterations of each point
set.seed(1)
x <- rnorm(P) * 10
y <- rnorm(P) * 10


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Generate the P paths from the random starting points
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
res <- seq(P) %>%
  map(~compute_ikeda_trajectory(u=0.928, x[.x], y[.x], N))
Warning: `data_frame()` is deprecated, use `tibble()`.
This warning is displayed once per session.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Plot all the paths
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
par(pty='s',bg='black')
plot(0, type='n', as=1, xlim=c(-10, 10), ylim=c(-10, 10))
res %>% walk(~lines(.x, col='#ffffff08'))

Tweetable Ikeda map

This one was a hell of a squeeze to get it into a tweet:

  • Use 1e3 instead of 1000 to save a single character.
  • Drop leading 0 from 0.9
  • Use rep(0,N) instead of double(N) to save a single character
  • use lst instead of tibble or data.frame to save at least 3 chars
  • Dropping from 10 to 9 as a multiplier to trip a single digit each time.
  • Came out at exactly 280 characters!
  • Update Thanks to Jonathan Carroll for pointing out a shorter way to call map2 it’s down to 272 characters which means I can include the #rstats hashtag! (which then brings it back up to 280 characters)
#rstats
library(tidyverse)
N=1e3
x=y=rep(0,N)
F=function(X,Y){x[1]=X
y[1]=Y
for(i in 2:N){X=x[i-1];Y=y[i-1]
t=.4-6/(1+X^2+Y^2)
S=sin(t);C=cos(t)
x[i]=1+.9*(X*C-Y*S)
y[i]=.9*(X*S+Y*C)}
lst(x,y)}
r=rnorm
plot(0,as=1,yli=c(-9,9))
walk(map2(r(N)*9,r(N)*9,F),~lines(.x,co='#00000010'))