###

polarCurve <- function(n, rphi = function(phi,A) A*phi,
                       turns = 2.5, sigma.x = 1, sigma.y = sigma.x,
                       ...)
{
    ## Purpose: return 2-column matrix points on a Polar coordinate curve

    ## Arguments: n: #{points};  rphi: r(phi) curve defining function
    ##		turns: phi will be in [0, turns*2*pi]
    ##		sigma.x, sigma.y: N(0,sigma) noise will be added to x & y
    ##          ... : further arguments passed to rphi(.)
    ## -------------------------------------------------------------------------
    ## Author: Martin Maechler, Date:  6 Dec 2002, 12:17
    if(0 > (n <- as.integer(n)))
        stop("`n' must be non-negative integer")
    if(0 == n)
        return(list(x=numeric(0), y = numeric(0)))
    ## else n >= 1
    phi <- seq(0, turns*2*pi, length = n)
    r   <- rphi(phi, ...)
    cbind(x = r * cos(phi) + rnorm(n, 0, sigma.x),
          y = r * sin(phi) + rnorm(n, 0, sigma.y))
}

plot((m1 <- polarCurve(20, A = 1)), type = 'o')

splineCurve <- function(x, y=NULL, n = 8*length(x),  method = "fmm")
{
  ## Purpose: Interpolating spline *curve* (not function)
  ## -------------------------------------------------------------------------
  ## Arguments: same as for spline()
  ## -------------------------------------------------------------------------
  ## Author: Martin Maechler, Date:  6 Dec 2002, 12:34
   xy <- xy.coords(x,y)
   nP <- length(xy$x)
   if(nP == 0)
       return(list(x=numeric(0), y = numeric(0)))
   ## else :
   i <- 1:nP
   Sx <- splinefun(i, xy$x, method = method)
   Sy <- splinefun(i, xy$y, method = method)
   ti <- seq(1, nP, length = n)
   cbind(x = Sx(ti), y = Sy(ti))
}

# plot(m1)
lines(sm1 <- splineCurve(m1), col = "blue")
lines(sm2 <- splineCurve(m1, method = "natural"), col = "dark red")
