drosofff
12/26/2016 - 10:34 PM

From http://menugget.blogspot.fr/2012/12/spirograph-with-r.html#more

gcd <- function(a, b)
{
  a=as.integer(a)
  b=as.integer(b)
  stopifnot(!(a==0 && b==0))
  if ( a == 0 )
    return(b)

  if ( b == 0 )
    return(a)

  return(abs(gcd(b, a-b*floor(a/b))))
}

lcm <- function(a, b)
{
  a=as.integer(a)
  b=as.integer(b)
  return(abs(a*b)/gcd(a,b))
}
RES <- vector(mode="list", 100)
PARAM <- vector(mode="list",100)
LIM <- c()
set.seed(1937) # remove or change seed for new randomness !
for(i in seq(RES)){
 #i=1
 a.rad=sample(seq(1,10, 0.1),1) 
 b.rad=sample(seq(-3,10, 0.1),1)
 bc=sample(seq(-10,10),1)
 a.rev=lcm(abs(a.rad), abs(b.rad))
 a.rev=a.rev+1
 tmp <- runif(2, min=-100, max=100)
 a.cen=list(x=tmp[1], y=tmp[2])
 PARAM[[i]]$x <-  tmp[1]
 PARAM[[i]]$y <-  tmp[2]
 PARAM[[i]]$P <- c(a.rad, b.rad, bc, a.rev)
 LIM <- range(c(LIM, unlist(a.cen)))
 RES[[i]] <- spirographR(A.RADIUS=a.rad, B.RADIUS=b.rad, BC=bc, A.REV=a.rev, A.CEN=a.cen)
}
# png("spirograph.png", width=4, height=4, units="in", res=600)
par(mar=c(0,0,0,0), bg=1)
for(i in seq(RES)){
 if(i == 1){
  plot(RES[[i]], t="n", xlim=LIM, ylim=LIM, asp=1)
 }
 lines(RES[[i]], col=rgb(runif(1,.6,1), 1, 1, runif(1,0.4,0.8)), lwd=0.4)
# text (PARAM[[i]]$x, PARAM[[i]]$y, label=toString(PARAM[[i]]$P), col="white", cex=.5)
}
#spirographR() will produce either a hypotrochoid or an epitrochoid. 
#'A' is a circle of radius 'A.RADIUS'
#'B' is a circle of radius 'B.RADIUS' travelling around 'A'
#'C' is a point relative to the center of 'B' which rotates with the turning of 'B'.
#'BC' is the distance from the center of 'B' to 'C'
#'A.REV' is the number of revolutions that 'B' should travel around 'A'
#'N.PER.A.REV' is the number of radial increments to be calculated per revolution
#'A.CEN' is the position of the center of 'A'
#
#NOTE: A positive value for 'B' will result in a epitrochoid, while a negative value will result in a hypotrochoid.
#
spirographR <- function(
A.RADIUS=1, 
B.RADIUS=-4, 
BC=-2, 
A.REV=4, 
N.PER.A.REV=360,
A.CEN=list(x=0, y=0)){
 
 B.CEN.START <- list(x=0, y=A.CEN$y + A.RADIUS + B.RADIUS) #starting position of B circle
 A.ANGLE <- seq(0, 2*pi*A.REV,, A.REV*N.PER.A.REV) #Radians around A for calculation
 A.CIR <- 2*pi*A.RADIUS #Circumference of A
 B.CIR <- 2*pi*B.RADIUS #Circumference of B
 
 ###Find position of B.CEN
 B.CEN <- c()
 HYP <- A.RADIUS + B.RADIUS
 ADJ <- sin(A.ANGLE) * HYP
 OPP <- cos(A.ANGLE) * HYP
 B.CEN$x <- A.CEN$x + ADJ 
 B.CEN$y <- A.CEN$y + OPP 
 
 ###Find position of C.POINT
 C.POINT <- c()
 A.CIR.DIST <- A.CIR * A.ANGLE / (2*pi)
 B.POINT.ANGLE <- A.CIR.DIST / B.CIR * 2*pi
 HYP <- BC
 ADJ <- sin(B.POINT.ANGLE) * HYP
 OPP <- cos(B.POINT.ANGLE) * HYP
 C.POINT$x <- B.CEN$x + ADJ 
 C.POINT$y <- B.CEN$y + OPP 
 
 ###Return trajectory of C
 C.POINT
}
RES <- vector(mode="list", 100)
PARAM <- vector(mode="list",100)
LIM <- c()
set.seed(7)
for(i in seq(RES)){
 #i=1
 a.rad=sample(seq(1,10, 0.1),1) 
 b.rad=sample(seq(-3,10, 0.1),1)
 bc=sample(seq(-10,10),1)
 a.rev=lcm(abs(a.rad), abs(b.rad))
 if (a.rev==0) {a.rev=1}
 tmp <- runif(2, min=-100, max=100)
 a.cen=list(x=tmp[1], y=tmp[2])
 PARAM[[i]]$x <-  tmp[1]
 PARAM[[i]]$y <-  tmp[2]
 PARAM[[i]]$P <- c(a.rad, b.rad, bc, a.rev)
 LIM <- range(c(LIM, unlist(a.cen)))
 RES[[i]] <- spirographR(A.RADIUS=a.rad, B.RADIUS=b.rad, BC=bc, A.REV=a.rev, A.CEN=a.cen)
}
# png("spirograph.png", width=4, height=4, units="in", res=600)
par(mar=c(0,0,0,0), bg=1)
for(i in seq(RES)){
 if(i == 1){
  plot(RES[[i]], t="n", xlim=LIM, ylim=LIM, asp=1)
 }
 lines(RES[[i]], col=rgb(runif(1), runif(1), runif(1), runif(1,0.4,0.8)), lwd=0.4)
 text (PARAM[[i]]$x, PARAM[[i]]$y, label=i, col="white", cex=1)
}

## then
Spiros <- vector(mode="list", 100)
Spiros[[1]] <- PARAM[[1]] # to pick up nice spirographs (up to 100)
gcd <- function(a, b)
{
  a=as.integer(a)
  b=as.integer(b)
  stopifnot(!(a==0 && b==0))
  if ( a == 0 )
    return(b)

  if ( b == 0 )
    return(a)

  return(abs(gcd(b, a-b*floor(a/b))))
}

lcm <- function(a, b)
{
  a=as.integer(a)
  b=as.integer(b)
  return(abs(a*b)/gcd(a,b))
}
 
RES <- vector(mode="list", 100)
PARAM <- vector(mode="list",100)
LIM <- c()
set.seed(7)
for(i in seq(RES)){
 #i=1
 a.rad=sample(seq(1,10, 0.1),1) 
 b.rad=sample(seq(-3,10, 0.1),1)
 bc=sample(seq(-10,10),1)
 a.rev=lcm(abs(a.rad), abs(b.rad))
 if (a.rev==0) {a.rev=1}
 tmp <- runif(2, min=-100, max=100)
 a.cen=list(x=tmp[1], y=tmp[2])
 PARAM[[i]]$x <-  tmp[1]
 PARAM[[i]]$y <-  tmp[2]
 PARAM[[i]]$P <- c(a.rad, b.rad, bc, a.rev)
 LIM <- range(c(LIM, unlist(a.cen)))
 RES[[i]] <- spirographR(A.RADIUS=a.rad, B.RADIUS=b.rad, BC=bc, A.REV=a.rev, A.CEN=a.cen)
}
# png("spirograph.png", width=4, height=4, units="in", res=600)
par(mar=c(0,0,0,0), bg=1)
for(i in seq(RES)){
 if(i == 1){
  plot(RES[[i]], t="n", xlim=LIM, ylim=LIM, asp=1)
 }
 lines(RES[[i]], col=rgb(runif(1), runif(1), runif(1), runif(1,0.4,0.8)), lwd=0.4)
 text (PARAM[[i]]$x, PARAM[[i]]$y, label=toString(PARAM[[i]]$P), col="white", cex=.5)
}
# dev.off()
# Spiros is a list of selected spirographs 

RES <- vector(mode="list", 100)
PARAM <- vector(mode="list",100)
LIM <- c()
Coordinates <-vector(mode="list",100)
for(i in seq(Spiros)){
 a.rad=Spiros[[i]]$P[1] 
 b.rad=Spiros[[i]]$P[2]
 bc=Spiros[[i]]$P[3]
 a.rev=Spiros[[i]]$P[4]
 tmp <- runif(2, min=-100, max=100)
 a.cen=list(x=tmp[1], y=tmp[2])
 Coordinates[[i]]$x = tmp[1]
 Coordinates[[i]]$y = tmp[2]
 PARAM[[i]]$P <- c(a.rad, b.rad, bc, a.rev)
 LIM <- range(c(LIM, unlist(a.cen)))
 RES[[i]] <- spirographR(A.RADIUS=a.rad, B.RADIUS=b.rad, BC=bc, A.REV=a.rev, A.CEN=a.cen)
}
# png("spirograph.png", width=4, height=4, units="in", res=600)
par(mar=c(0,0,0,0), mfrow=c(10,10), oma = c(0,0,0,0), bg=1)
for(i in seq(RES)){
 plot(RES[[i]], t="n", axes=F,xlim=c(min(RES[[i]]$x),max(RES[[i]]$x)), ylim=c(min(RES[[i]]$y),max(RES[[i]]$y)), asp=1)
 lines(RES[[i]], col=rgb(runif(1,.6,1), 1, 1, runif(1,0.4,0.8)), lwd=0.4)
 text (mean(RES[[i]]$x),mean(RES[[i]]$y), label=toString(PARAM[[i]]$P), col="red", cex=1)
}