Roots of Unity in R

Question

I have a paper on the definition of (i) available here. I was curious about the simultaneous real numbers associated with the roots of unity, and noticed there was not a routine in R.

I have this code and would like to know any efficiency to be had in either the derivation or plotting. Thanks for any suggestions and insights!

Roots of Unity:

Roots.of.unity <- function(nth.root = N, radius=1, simultaneous.reals=TRUE){

  Simultaneous.reals = data.frame()

  Complex.points = numeric()
  Reals = numeric()
  Imaginaries = numeric()

  r = radius^(1/nth.root)


  for(j in 0:(nth.root-1)){
    k= ((360*j)/nth.root)*pi/180

    Simultaneous.reals[(j+1),1] = r*cos(k) - r*sin(k)
    Simultaneous.reals[(j+1),2] = r*cos(k) + r*sin(k)

    Complex.points[j+1] = complex(real = r*cos(k), imaginary = r*sin(k))
    Complex.points = c(Complex.points,Complex.points[j+1])

    Reals[j+1] = r*cos(k)
    Imaginaries[j+1] = r*sin(k)

  }

  circleFun <- function(center = c(0,0),diameter = 2*radius, npoints = 100){
    r = diameter / 2
    tt <- seq(0,2*pi,length.out = npoints)
    xx <- center[1] + r * cos(tt)
    yy <- center[2] + r * sin(tt)
    return(data.frame(x = xx, y = yy))
  }
  dat <- circleFun(c(0,0),2*radius^(1/nth.root),npoints = 100)

  plot(dat,type = 'l',main = substitute(paste("roots of " , "Z"^nth.root ,"= ", Z),list(nth.root=nth.root,Z=radius)),
   xlab = "Real", ylab = "Imaginary",
   xlim = c(min(Simultaneous.reals[,1]),max(Simultaneous.reals[,1])),
   ylim = c(min(Simultaneous.reals[,1]),max(Simultaneous.reals[,1]))
   )
  if(simultaneous.reals==TRUE){
      legend('topleft',c("Complex Points","Simultaneous Reals"),pch = c(19,19), col=c('red','blue'),bty='n')

      points(Simultaneous.reals[,1],rep(0,length(Simultaneous.reals[,1])),
     pch=19,
     col='blue')
      points(Reals,Imaginaries,col='red',pch=19)


  } else {

           legend('topleft',"Complex Points",pch = 19, col='red',bty='n')

  points(Reals,Imaginaries,col='red',pch=19)}

  colnames(Simultaneous.reals) = c("Real 1","Real 2")
  return(cbind(Simultaneous.reals,"Complex Points" = Complex.points[1:nth.root]))


}  

> Roots.of.unity(15)
    Real 1      Real 2        Complex Points
1   1.00000000  1.00000000  1.0000000+0.0000000i
2   0.50680881  1.32028210  0.9135455+0.4067366i
3  -0.07401422  1.41227543  0.6691306+0.7431448i
4  -0.64203952  1.26007351  0.3090170+0.9510565i
5  -1.09905036  0.88999343 -0.1045285+0.9945219i
6  -1.36602540  0.36602540 -0.5000000+0.8660254i
7  -1.39680225 -0.22123174 -0.8090170+0.5877853i
8  -1.18605929 -0.77023591 -0.9781476+0.2079117i
9  -0.77023591 -1.18605929 -0.9781476-0.2079117i
10 -0.22123174 -1.39680225 -0.8090170-0.5877853i
11  0.36602540 -1.36602540 -0.5000000-0.8660254i
12  0.88999343 -1.09905036 -0.1045285-0.9945219i
13  1.26007351 -0.64203952  0.3090170-0.9510565i
14  1.41227543 -0.07401422  0.6691306-0.7431448i
15  1.32028210  0.50680881  0.9135455-0.4067366i