.packageName <- "circular"
###############################################################
#                                                             #
#       R port: Claudio Agostinelli  <claudio@unive.it>       #
#                                                             #
#       Date: January, 14, 2003                               #
#       Version: 0.1-6                                        #
#                                                             #
###############################################################

A1 <- function(kappa) {
    result <- besselI(kappa, nu=1, expon.scaled = TRUE)/besselI(kappa, nu=0, expon.scaled = TRUE)
    return(result)
}
###############################################################
#                                                             #
#       R port: Claudio Agostinelli  <claudio@unive.it>       #
#                                                             #
#       Date: January, 14, 2003                               #
#       Version: 0.1-1                                        #
#                                                             #
###############################################################

A1inv <- function(x) {
   ifelse (0 <= x & x < 0.53, 2 * x + x^3 + (5 * x^5)/6,
           ifelse (x < 0.85, -0.4 + 1.39 * x + 0.43/(1 - x), 1/(x^3 - 4 * x^2 + 3 * x)))
}

I.0 <- function(x) {
    besselI(x=x, nu=0, expon.scaled = FALSE)
}

I.1 <- function(x) {
    besselI(x=x, nu=1, expon.scaled = FALSE)
}

I.p <- function(p, x) {
    besselI(x=x, nu=p, expon.scaled = FALSE)
}
#############################################################
#                                                           #
#   as.circular function                                    #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

as.circular <- function (x, ...) {
    if (is.circular(x)) return(x)
    else if(!is.null(xcircularp <- circularp(x))) circular(x, type=xcircularp$type, units=xcircularp$units, template=xcircularp$template, modulo=xcircularp$modulo, zero=xcircularp$zero, rotation=xcircularp$rotation)
    else circular(x, ...)
}

#############################################################
#                                                           #
#   as.data.frame.circular function                         #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: September, 22, 2003                               #
#   Version: 0.1-2                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

as.data.frame.circular <- function(x, row.names=NULL, optional=FALSE) {
    if (is.matrix(x)) {
        if (!is.null(xcircularp <- circularp(x))) {
            typep <- xcircularp$type
            unitsp <- xcircularp$units
            templatep <- xcircularp$template
            modulop <- xcircularp$modulo
            zerop <- xcircularp$zero
            rotationp <- xcircularp$rotation
        } else {
            typep <- "angles"
            unitsp <- "radians"
            templatep <- "none"
            modulop <- "asis"
            zerop <- 0
            rotationp <- "counter"
        }

        d <- dim(x)
        nrows <- d[1]; ir <- seq(length = nrows)
        ncols <- d[2]; ic <- seq(length = ncols)
        dn <- dimnames(x)
        row.names <- dn[[1]]
        collabs <- dn[[2]]
        if (any(empty <- nchar(collabs)==0))
        collabs[empty] <- paste("Circular", ic, sep = "")[empty]
        value <- vector("list", ncols)
    for(i in ic)
        value[[i]] <- as.circular(x[,i], type=typep, units=unitsp, modulo=modulop, zero=zerop, rotation=rotationp)
        if (length(row.names) != nrows)
        row.names <- if(optional) character(nrows) else as.character(ir)
        if (length(collabs) == ncols)
        names(value) <- collabs
        else if(!optional)
            names(value) <- paste("Circular", ic, sep="")
        attr(value, "row.names") <- row.names
        class(value) <- "data.frame"
        return(value)

    } else
    return(as.data.frame.vector(x, row.names, optional))
}
 

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rcardioid function                                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: April, 29, 2003                                   #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

rcardioid <- function(n, mu=0, rho=0, units=c("radians", "degrees"), ...) {
    units <- match.arg(units)
    if (units=="degrees") {
        mu <- mu/180*pi
    }

    if (rho < -0.5 | rho > 0.5)
        stop("rho must be between -0.5 and 0.5")        
    i <- 1
    result <- rep(0, n)
    while(i <= n) {
        x <- runif(1, 0, 2 * pi)
        y <- runif(1, 0, (1 + 2 * rho)/(2 * pi))
        f <- (1 + 2 * rho * cos(x - mu))/(2 * pi)
        if(y <= f) {
            result[i] <- x
            i <- i + 1
        }
    }    
    if (units=="degrees") result <- result/pi*180
    result <- circular(result, units=units, ...)
    return(result)
}

#############################################################
#                                                           #
#   dcardioid function                                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 23, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-1                                           #
#############################################################

dcardioid <- function(x, mu=0, rho=0) {
  
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    attr(x, "circularp") <-  NULL
    
    if (units=="degrees") {
        mu <- mu/180*pi
    }  

    if (rho < -0.5 | rho > 0.5)
        stop("rho must be between -0.5 and 0.5")
    d <- (1 + 2 * rho * cos(x - mu))/(2 * pi)
    d <- unclass(d)
    return(d)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   change.point function                                   #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

change.point <- function(x) {
        x <- as.circular(x)
        xcircularp <- circularp(x)
        units <- xcircularp$units
        x <- conversion.circular(x, units="radians")
    phi <- function(x) {
        arg <- A1inv(x)
        if(besselI(x=arg, nu=0, expon.scaled = FALSE) != Inf)
            result <- x * A1inv(x) - log(besselI(x=arg, nu=0, expon.scaled = FALSE))
        else result <- x * A1inv(x) - (arg + log(1/sqrt(2 * pi * arg) * (1 + 1/(8 * arg) + 9/(128 * arg^2) + 225/(1024 * arg^3))))
        result
    }
    n <- length(x)
    rho <- rho.circular(x)
    R1 <- c(1:n)
    R2 <- c(1:n)
    V <- c(1:n)
    for(k in 1:(n - 1)) {
        R1[k] <- rho.circular(x[1:k]) * k
        R2[k] <- rho.circular(x[(k + 1):n]) * (n - k)
        if(k >= 2 & k <= (n - 2)) {
            V[k] <- k/n * phi(R1[k]/k) + (n - k)/n * phi(R2[k]/(n - k))
        }
    }
    R1[n] <- rho * n
    R2[n] <- 0
    R.diff <- R1 + R2 - rho * n
    rmax <- max(R.diff)
    rave <- mean(R.diff)
    k.r <- (1:n)[R.diff == max(R.diff)]
    V <- V[2:(n - 2)]
    if(n > 3) {
        tmax <- max(V)
        tave <- mean(V)
        k.t <- (1:(n - 3))[V == max(V)] + 1
    }
    else stop("Sample size must be at least 4")
    return(list(n=n, rho=rho, rmax=rmax, k.r=k.r, rave=rave, tmax=tmax, k.t=k.t, tave=tave))
}
#############################################################
#                                                           #
#   circular function                                       #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: September, 22, 2003                               #
#   Version: 0.6-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

circular <- function(x, type=c("angles", "directions"), units=c("radians", "degrees"), template=c("none", "geographics"), modulo=c("asis", "2pi", "pi"), zero=0, rotation=c("counter", "clock"), names) {

    type <- match.arg(type)
    units <- match.arg(units)
    template <- match.arg(template) 
    modulo <- match.arg(modulo)
    rotation <- match.arg(rotation)
  
    if (template=="geographics") {
        zero <- pi/2
        rotation <- "clock"
    }

    if (is.data.frame(x)) x <- as.matrix(x)

    if (is.matrix(x)) {
    nseries <- ncol(x)
    ndata <- nrow(x)
    if (missing(names)) {
            names <- if(!is.null(dimnames(x))) colnames(x) else paste("Circular", seq(nseries), sep="")
    }
        dimnames(x) <- list(NULL, names)
    } else {
        nseries <- 1
    ndata <- length(x)
    }
####    if (ndata == 0) stop("circular object must have one or more observations")
    if (modulo!="asis") {
    if (modulo=="2pi") {
        ang <- 2
    } else {
        ang <- 1
    }
    if (units=="radians") {
        x <- x %% (ang*pi)
    } else {
        x <- x %% (ang*180)
    }
    }

    attr(x, "circularp") <- list(type=type, units=units, template=template, modulo=modulo, zero=zero, rotation=rotation) #-- order is fixed
    attr(x, "class") <- "circular"
    return(x)
}

#############################################################
#                                                           #
#       conversion.circular function                        #
#   Author: Claudio Agostinelli                         #
#   E-mail: claudio@unive.it                            #
#   Date: July, 31, 2003                                #
#   Version: 0.1-1                                      #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################

conversion.circular <- function(x, units=c("radians", "degrees")) {
    units <- match.arg(units)
    x <- as.circular(x)
    value <- attr(x, "circularp")
    unitsp <- value$units
          
    if (unitsp=="degrees" & units=="radians") {
    x <- x/180*pi
    } else if (unitsp=="radians" & units=="degrees") {
               x <- x/pi*180
    }
    value$units <- units
    circularp(x) <- value 
    
    return(x)
}

#############################################################
#                                                           #
#   circularp function                                  #
#   Author: Claudio Agostinelli                         #
#   E-mail: claudio@unive.it                            #
#   Date: March, 7, 2003                                #
#   Version: 0.1                                        #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################
 
circularp <- function(x) attr(x, "circularp")

#############################################################
#                                                           #
#   circularp<- function                                    #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 18, 2003                                #
#   Version: 0.2-2                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################
 
"circularp<-" <- function(x, value) {
    cl <- class(x)
    if (length(value)!=6) stop("value must have six elements")

    type <- value$type
    units <- value$units
    template <- value$template
    modulo <- value$modulo
    zero <- value$zero
    rotation <- value$rotation
 
    if (type!="angles" & type!="directions") stop("type (value[1]) must be 'angles', 'directions' or 'geographics'")

    if (units!="radians" & units!="degrees") stop("units (value[2]) must be 'radians' or 'degrees'")

    if (template!="none" & template!="geographics") stop("template (value[3]) must be 'none' or 'geographics'")
    
    if (modulo!="asis" & modulo!="2pi" &  modulo!="pi") stop("modulo (value[4]) must be 'asis' or 'pi' or '2pi'")

    if (rotation!="clock" & rotation!="counter") stop("rotation (value[6]) must be 'clock' or 'counter'")

    attr(x, "circularp") <- value
    if (inherits(x, "circular") && is.null(value))
        class(x) <- cl["circular" != cl]
    return(x)
}

#############################################################
#                                                           #
#   is.circular function                                #
#   Author: Claudio Agostinelli                         #
#   E-mail: claudio@unive.it                            #
#   Date: March, 7, 2003                                #
#   Version: 0.1                                        #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################
 
is.circular <- function (x) inherits(x, "circular")

#############################################################
#                                                           #
#   [.circular function                                 #
#   Author: Claudio Agostinelli                         #
#   E-mail: claudio@unive.it                            #
#   Date: March, 7, 2003                                #
#   Version: 0.1                                        #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################
 
"[.circular" <- function(x, i, ...) {
    y <- NextMethod("[", ...)
    class(y) <- class(x)
    attr(y, "circularp") <- attr(x, "circularp")
    return(y)
}

#############################################################
#                                                           #
#   print.circular function                             #
#   Author: Claudio Agostinelli                         #
#   E-mail: claudio@unive.it                            #
#   Date: June, 21, 2003                                #
#   Version: 0.2                                        #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################
 
print.circular <- function(x, info=TRUE, ...) {
    x.orig <- x
    x <- as.circular(x)
    if (info) {
    xcircularp <- attr(x, "circularp")
    type <- xcircularp$type
    units <- xcircularp$units
        template <- xcircularp$template
    modulo <- xcircularp$modulo
        zero <- xcircularp$zero
        rotation <- xcircularp$rotation

        cat("Circular Data: \nType =", type,
               "\nUnits =", units,
               "\nTemplate =", template, 
               "\nModulo =", modulo,
               "\nZero =", zero,
               "\nRotation =", rotation, "\n")
    }       
    attr(x, "class") <- attr(x, "circularp") <- attr(x, "na.action") <- NULL
    NextMethod("print", x, quote = FALSE, right = TRUE, ...)
    invisible(x.orig)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

###############################################################
#                                                             #
#       R port: Claudio Agostinelli  <claudio@unive.it>       #
#                                                             #
#       Date: January, 14, 2003                               #
#       Version: 0.1-6                                        #
#                                                             #
###############################################################

cor.circular <- function(x, y=NULL, test = FALSE) {
        ncx <- NCOL(x)
        ncy <- NCOL(y)        
        n <- NROW(x)
        x <- as.circular(x)
        x <- conversion.circular(x, units="radians")
        if (!is.null(y)) { 
            y <- as.circular(y)
            y <- conversion.circular(y, units="radians")
        }
        if (!is.null(y) & NROW(x)!=NROW(y)) 
            stop("x and y must have the same number of observations")
        if (is.null(y) & ncx<2)
            stop("supply both x and y or a matrix-like x")

        cor.internal <- function(x, y, test=FALSE) {
       x.bar <- mean.circular(x)
       y.bar <- mean.circular(y)
       num <- sum(sin(x - x.bar) * sin(y - y.bar))
       den <- sqrt(sum(sin(x - x.bar)^2) * sum(sin(y - y.bar)^2))
       result <- num/den
           if (test) {
           l20 <- mean.default(sin(x - x.bar)^2)
           l02 <- mean.default(sin(y - y.bar)^2)
           l22 <- mean.default((sin(x - x.bar)^2) * (sin(y - y.bar)^2))
           test.stat <- sqrt((n * l20 * l02)/l22) * result
           p.value <- 2 * (1 - pnorm(abs(test.stat)))
           result <- c(result, test.stat, p.value)
           }
       return(result)
        }

        if (is.null(y)) {
            result <- matrix(1, ncol=ncx, nrow=ncx)
            if (test) {
                test.stat <- matrix(0, ncol=ncx, nrow=ncx)
                p.value <- matrix(0, ncol=ncx, nrow=ncx)
            }
            for (i in 1:ncx) {
                 for (j in i:ncx) {
                      res <- cor.internal(x=x[,i], y=x[,j], test=test)
                      result[i,j] <- result[j,i] <- res[1]
                      if (test) {
                          test.stat[i,j] <- test.stat[j,i] <- res[2] 
                          p.value[i,j] <- p.value[j,i] <- res[3]
                      }
                 }
            }            
        } else {
            attributes(x) <- c(attributes(x), list(dim=c(n, ncx)))
            attributes(y) <- c(attributes(y), list(dim=c(n, ncy)))
            result <- matrix(1, ncol=ncy, nrow=ncx)
            if (test) {
                test.stat <- matrix(0, ncol=ncy, nrow=ncx)
                p.value <- matrix(0, ncol=ncy, nrow=ncx)
            }
            for (i in 1:ncx) {
                 for (j in 1:ncy) {
                      res <- cor.internal(x=x[,i], y=y[,j], test=test)
                      result[i,j] <- res[1]
                      if (test) {
                          test.stat[i,j] <- res[2] 
                          p.value[i,j] <- res[3]
                      }
                 }
            }
       }
    if (ncx==1 | (!is.null(y) & ncy==1)) {
        result <- c(result)
        if (test) {
            test.stat <- c(test.stat)
            p.value <- c(p.value)
        }
    }


    if (test) {
        result <- list(cor=result, statistic=test.stat, p.value=p.value)
    }

    return(result)

}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

###############################################################
#                                                             #
#       R port: Claudio Agostinelli  <claudio@unive.it>       #
#                                                             #
#       Date: January, 14, 2003                               #
#       Version: 0.1-6                                        #
#                                                             #
###############################################################

deg <- function(x) {
    (x * 180)/pi
}

density <- function(x, ...) UseMethod("density")

density.default <- function(x, bw = "nrd0", adjust = 1, kernel = c("gaussian",
    "epanechnikov", "rectangular", "triangular", "biweight",
    "cosine", "optcosine"), window = kernel, width, give.Rkern = FALSE,
    n = 512, from, to, cut = 3, na.rm = FALSE, ...) {
###    base::density(x=x, bw=bw, adjust=adjust, kernel=kernel, window=window, width=width, give.Rkern=give.Rkern, n=n, from=from, to=to, cut=cut, na.rm=na.rm)
### The function is copied from version R 1.9.0 of 2003/10/27 (unstable) since it is not possible to use the above redirection for problems in passing the parameter 'bw'

  if(!missing(window) && missing(kernel))
        kernel <- window
    kernel <- match.arg(kernel)
    if(give.Rkern)
        ##-- sigma(K) * R(K), the scale invariant canonical bandwidth:
        return(switch(kernel,
                      gaussian = 1/(2*sqrt(pi)),
                      rectangular = sqrt(3)/6,
                      triangular  = sqrt(6)/9,
                      epanechnikov= 3/(5*sqrt(5)),
                      biweight    = 5*sqrt(7)/49,
                      cosine      = 3/4*sqrt(1/3 - 2/pi^2),
                      optcosine   = sqrt(1-8/pi^2)*pi^2/16
                      ))

    if (!is.numeric(x))
        stop("argument must be numeric")
    name <- deparse(substitute(x))
    x <- as.vector(x)
    x.na <- is.na(x)
    if (any(x.na)) {
        if (na.rm) x <- x[!x.na]
        else stop("x contains missing values")
    }
    N <- nx <- length(x)
    x.finite <- is.finite(x)
    if(any(!x.finite)) {
        x <- x[x.finite]
        nx <- sum(x.finite)
    }
    n.user <- n
    n <- max(n, 512)
    if (n > 512) n <- 2^ceiling(log2(n)) #- to be fast with FFT

    if (missing(bw) && !missing(width)) {
        if(is.numeric(width)) {
            ## S has width equal to the length of the support of the kernel
            ## except for the gaussian where it is 4 * sd.
            ## R has bw a multiple of the sd.
            fac <- switch(kernel,
                          gaussian = 4,
                          rectangular = 2*sqrt(3),
                          triangular = 2 * sqrt(6),
                          epanechnikov = 2 * sqrt(5),
                          biweight = 2 * sqrt(7),
                          cosine = 2/sqrt(1/3 - 2/pi^2),
                          optcosine = 2/sqrt(1-8/pi^2)
                          )
            bw <- width / fac
        }
        if(is.character(width)) bw <- width
    }
    if (is.character(bw)) {
        if(length(x) < 2)
            stop("need at least 2 points to select a bandwidth automatically")
        bw <- switch(tolower(bw),
                     nrd0 = bw.nrd0(x),
                     nrd = bw.nrd(x),
                     ucv = bw.ucv(x),
                     bcv = bw.bcv(x),
                     sj = , "sj-ste" = bw.SJ(x, method="ste"),
                     "sj-dpi" = bw.SJ(x, method="dpi"),
                     stop("unknown bandwidth rule"))
    }
    if (!is.finite(bw)) stop("non-finite `bw'")
    bw <- adjust * bw
    if (bw <= 0) stop("`bw' is not positive.")

    if (missing(from))
        from <- min(x) - cut * bw
    if (missing(to))
	to   <- max(x) + cut * bw
    if (!is.finite(from)) stop("non-finite `from'")
    if (!is.finite(to)) stop("non-finite `to'")
    lo <- from - 4 * bw
    up <- to + 4 * bw
    y <- .C("massdist",
	    x = as.double(x),
	    nx = nx,
	    xlo = as.double(lo),
	    xhi = as.double(up),
	    y = double(2 * n),
	    ny = as.integer(n),
	    PACKAGE = "base")$y * (nx/N)
    kords <- seq(0, 2*(up-lo), length = 2 * n)
    kords[(n + 2):(2 * n)] <- -kords[n:2]
    kords <- switch(kernel,
		    gaussian = dnorm(kords, sd = bw),
                    ## In the following, a := bw / sigma(K0), where
                    ##	K0() is the unscaled kernel below
		    rectangular = {
                        a <- bw*sqrt(3)
                        ifelse(abs(kords) < a, .5/a, 0) },
		    triangular = {
                        a <- bw*sqrt(6) ; ax <- abs(kords)
                        ifelse(ax < a, (1 - ax/a)/a, 0) },
		    epanechnikov = {
                        a <- bw*sqrt(5) ; ax <- abs(kords)
                        ifelse(ax < a, 3/4*(1 - (ax/a)^2)/a, 0) },
		    biweight = { ## aka quartic
                        a <- bw*sqrt(7) ; ax <- abs(kords)
                        ifelse(ax < a, 15/16*(1 - (ax/a)^2)^2/a, 0) },
		    cosine = {
                        a <- bw/sqrt(1/3 - 2/pi^2)
                        ifelse(abs(kords) < a, (1+cos(pi*kords/a))/(2*a),0)},
		    optcosine = {
                        a <- bw/sqrt(1-8/pi^2)
                        ifelse(abs(kords) < a, pi/4*cos(pi*kords/(2*a))/a, 0)}
                    )
    kords <- fft( fft(y)* Conj(fft(kords)), inv=TRUE)
    kords <- Re(kords)[1:n]/length(y)
    xords <- seq(lo, up, length = n)
#    keep <- (xords >= from) & (xords <= to)
    x <- seq(from, to, length = n.user)
    structure(list(x = x, y = approx(xords, kords, x)$y, bw = bw, n = N,
		   call=match.call(), data.name=name, has.na = FALSE),
	      class="density")
  }

#############################################################
#                                                           #
#   density.circular function                               #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: November, 19, 2003                                #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-5                                           #
#                                                           #
#############################################################

density.circular <- function(x, z, bw, adjust = 1, type = c("K", "L"), kernel = c("vonmises", "wrappednormal"), na.rm = FALSE, from=0, to=2*pi, n=512, K=10, ...) {

    name <- deparse(substitute(x))
    xx <- x
    x <- as.circular(x)

    xcircularp <- circularp(x)
    type <- xcircularp$type
    units <- xcircularp$units
    template <- xcircularp$template
    zero <- xcircularp$zero
    rotation <- xcircularp$rotation
    
    x <- conversion.circular(x, units="radians")
  
    kernel <- match.arg(kernel)

    if (!is.numeric(n))
        stop("argument must be numeric")
    n <- round(n)
    if (n <=0)
         stop("argument must be integer and positive")     

    if (!is.numeric(from))
        stop("argument must be numeric")      
    if (!is.numeric(to))
        stop("argument must be numeric")      
    if (!is.finite(from)) 
        stop("non-finite `from'")
    if (!is.finite(to)) 
        stop("non-finite `to'")
    
    if (!is.numeric(x)) 
        stop("argument must be numeric")
    x <- as.vector(x)
    x.na <- is.na(x)
    if (any(x.na)) {
        if (na.rm) 
            x <- x[!x.na]
        else stop("x contains missing values")
    }
    
    nx <- length(x)
    x.finite <- is.finite(x)
    if (any(!x.finite)) {
        x <- x[x.finite]
        nx <- sum(x.finite)
    }

    if (missing(z)) {
        z <- zz <- circular(seq(from=from, to=to, length=n), units="radians", type=type, template=template, zero=zero, rotation=rotation)
        
    } else {
        z <- zz <- as.circular(z)
        z <- conversion.circular(z, units="radians")
        if (!is.numeric(z))
            stop("argument 'z' must be numeric")
        namez <- deparse(substitute(z))
        z <- as.vector(z)
        z.na <- is.na(z)
        if (any(z.na)) {
            if (na.rm) {
                z <- z[!z.na]
            } else {
                stop("z contains missing values")
            }
        }
    
        nz <- length(z)
        z.finite <- is.finite(z)
        if (any(!z.finite)) {
            z <- z[z.finite]
            nz <- sum(z.finite)
        }
    }

    bw <- adjust * bw
    if (!is.numeric(bw))
        stop("argument must be numeric")        
    if (!is.finite(bw)) 
        stop("non-finite `bw'")
    if (bw <= 0) 
        stop("`bw' is not positive.")

    if (kernel=="vonmises") {
        y <- sapply(z, dvonmises, mu=x, kappa=bw)
    } else if (kernel=="wrappednormal") {
        y <- sapply(z, dwrappednormal, mu=x, sd=bw, K=K)
    } else {
        stop("other kernels not implemented yet")
    }
    y <- apply(y, 2, sum)/nx

    if (units=="degrees") xx <- conversion.circular(xx, units="degrees")

    structure(list(data = xx, x = zz, y = y, bw = bw, n = nx, kernel=kernel, call = match.call(), data.name=name, has.na = FALSE), class = "density.circular")
} 

#############################################################
#                                                           #
#   plot.density.circular function                          #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: August, 01, 2003                                  #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.2-1                                           #
#                                                           #
#############################################################

plot.density.circular <- function(x, main = NULL, xlab = NULL, ylab = "Density circular", type = "l", zero.line = TRUE, points.plot=FALSE, points.col=1, points.pch=1, plot.type = c("circle", "line"), axes=TRUE, ticks=TRUE, bins, shrink=1, tcl=0.025, tol = 0.04, uin, xlim=c(-1, 1), ylim=c(-1, 1), ...) {

    x$x <- conversion.circular(x$x, units="radians")
    x$data <- conversion.circular(x$data, units="radians")
  
    plot.type <- match.arg(plot.type)
    if (missing(bins)) {
	bins <- NROW(x)
    } else {
	bins <- round(bins)
	if (bins<=0) stop("bins must be non negative")
    }
    
    if (is.null(xlab)) 
        xlab <- paste("N =", x$n, "  Bandwidth =", formatC(x$bw))
    if (is.null(main)) 
        main <- deparse(x$call)

    if (plot.type == "line") {
        xorder <- order(x$x)
        x$x <- x$x[xorder]
        x$y <- x$y[xorder] 
        plot.default(x, main = main, xlab = xlab, ylab = ylab, type = type, ...)
        if (zero.line) 
            abline(h = 0, lwd = 0.1, col = "gray")
        if (points.plot)
            points(x$data, rep(min(x$y),length(x$data)), col=points.col, pch=points.pch)
    } else {
        x$x <- as.circular(x$x)
        xcircularp <- attr(x$x, "circularp")
        xtype <- xcircularp$type
        units <- xcircularp$units
        template <- xcircularp$template
        modulo <- xcircularp$modulo
        zero <- xcircularp$zero
        rotation <- xcircularp$rotation
    
        x$x <- conversion.circular(x$x, units="radians")

        xlim <- shrink * xlim
        ylim <- shrink * ylim
        midx <- 0.5 * (xlim[2] + xlim[1])
        xlim <- midx + (1 + tol) * 0.5 * c(-1, 1) * (xlim[2] - xlim[1])
        midy <- 0.5 * (ylim[2] + ylim[1])
        ylim <- midy + (1 + tol) * 0.5 * c(-1, 1) * (ylim[2] - ylim[1])
        oldpin <- par("pin")
        xuin <- oxuin <- oldpin[1]/diff(xlim)
        yuin <- oyuin <- oldpin[2]/diff(ylim)
       if (missing(uin)) {
           if (yuin > xuin) yuin <- xuin
           else xuin <- yuin
       } else {
           if (length(uin) == 1) uin <- uin * c(1, 1)
           if (any(c(xuin, yuin) < uin)) stop("uin is too large to fit plot in")
           xuin <- uin[1]; yuin <- uin[2]
       }
       xlim <- midx + oxuin/xuin * c(-1, 1) * diff(xlim) * 0.5
       ylim <- midy + oyuin/yuin * c(-1, 1) * diff(ylim) * 0.5
       plot(cos(seq(0, 2 * pi, length = 1000)), sin(seq(0, 2 * pi, length = 1000)), axes = FALSE, xlab = "", ylab = "", main = "", type = "l", xlim=xlim, ylim=ylim, xaxs="i", yaxs="i")
 
        if (rotation=="clock") x$x <- -x$x
        x$x <- x$x + zero
            
        if (axes) {
	    axis.circular(units = units, template=template, modulo = modulo, zero=zero, rotation=rotation)
        }
 
        if (ticks) {
            at <- (0:bins)/bins*2*pi
            if (rotation=="clock") at <- -at
            at <- at + zero

            ticks.circular(circular(x=at, type="angles", units="radians", modulo="asis", zero=zero, rotation=rotation), tcl=tcl)
        }

        z <- (x$y+1)*cos(x$x)
        y <- (x$y+1)*sin(x$x)
        xorder <- order(x$x)
        z <- z[xorder]
        y <- y[xorder] 
        lines(x=z, y=y, type = type, ...)
        if (points.plot) {
            points.circular(x$data, col=points.col, pch=points.pch)
        }
    }
}

#############################################################
#                                                           #
#   lines.density.circular function                         #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 23, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#                                                           #
#############################################################

lines.density.circular <- function(x, type = "l", zero.line = TRUE, points.plot=FALSE, points.col=1, points.pch=1, plot.type = c("circle", "line"), bins, shrink=1, tcl=0.025, ...) {

    x$x <- conversion.circular(x$x, units="radians")
    x$data <- conversion.circular(x$data, units="radians")
  
    plot.type <- match.arg(plot.type)
    if (missing(bins)) {
	bins <- NROW(x)
    } else {
	bins <- round(bins)
	if (bins<=0) stop("bins must be non negative")
    }
    
    if (plot.type == "line") {
        xorder <- order(x$x)
        x$x <- x$x[xorder]
        x$y <- x$y[xorder] 
        lines.default(x, type = type, ...)
        if (zero.line) 
            abline(h = 0, lwd = 0.1, col = "gray")
        if (points.plot)
            points(x$data, rep(min(x$y),length(x$data)), col=points.col, pch=points.pch)
    } else {
        x$x <- as.circular(x$x)
        xcircularp <- attr(x$x, "circularp")
        xtype <- xcircularp$type
        units <- xcircularp$units
        template <- xcircularp$template
        modulo <- xcircularp$modulo
        zero <- xcircularp$zero
        rotation <- xcircularp$rotation
    
        x$x <- conversion.circular(x$x, units="radians")

        if (rotation=="clock") x$x <- -x$x
        x$x <- x$x + zero
            
        z <- (x$y+1)*cos(x$x)
        y <- (x$y+1)*sin(x$x)
        xorder <- order(x$x)
        z <- z[xorder]
        y <- y[xorder] 
        lines(x=z, y=y, type = type, ...)
        if (points.plot) {
            points.circular(x$data, col=points.col, pch=points.pch)
        }
    }
}


#############################################################
#                                                           #
#   print.density.circular function                         #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 23, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#                                                           #
#############################################################

print.density.circular <- function(x, digits=NULL, ...)
{
    cat("\nCall:\n\t",deparse(x$call),
	"\n\nData: ",x$data.name," (",x$n," obs.);",
	"\tBandwidth 'bw' = ",formatC(x$bw,digits=digits), "\n\n",sep="")
    print(summary(as.data.frame(x[c("x","y")])), digits=digits, ...)
    invisible(x)
}


###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   kuiper.test function                                    #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 25, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

kuiper.test <- function(x, alpha=0) {
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    if (!any(c(0, 0.01, 0.025, 0.05, 0.1, 0.15)==alpha)) stop("'alpha' must be one of the following values: 0, 0.01, 0.025, 0.05, 0.1, 0.15")
    x <- sort(x %% (2 * pi))/(2 * pi)
    n <- length(x)
    i <- 1:n
    D.P <- max(i/n - x)
    D.M <- max(x - (i - 1)/n)
    V <- (D.P + D.M) * (sqrt(n) + 0.155 + 0.24/sqrt(n))
    result <- list(statistic=V, alpha=alpha)
    class(result) <- "kuiper.test"
    return(result)
}

#############################################################
#                                                           #
#   print.kuiper.test function                              #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.kuiper.test <- function(x, digits=4, ...) {
    V <- x$statistic
    alpha <- x$alpha
    kuiper.crits <- cbind(c(0.15, 0.1, 0.05, 0.025, 0.01), c(1.537, 1.62, 1.747, 1.862, 2.001))
    cat("\n", "      Kuiper's Test of Uniformity", "\n", "\n")
    cat("Test Statistic: ", round(V, digits=digits), "\n")
    if (alpha == 0) {
        if (V < 1.537)
        cat("P-value > 0.15", "\n", "\n")
    else if (V < 1.62)
         cat("0.10 < P-value < 0.15", "\n", "\n")
    else if (V < 1.747)
         cat("0.05 < P-value < 0.10", "\n", "\n")
    else if (V < 1.862)
         cat("0.025 < P-value < 0.05", "\n", "\n")
    else if (V < 2.001)
         cat("0.01 < P-value < 0.025", "\n", "\n")
    else cat("P-value < 0.01", "\n", "\n")
    } else {
    Critical <- kuiper.crits[(1:5)[alpha == c(kuiper.crits[, 1])],2]
    cat("Level", alpha, "Critical Value:", round(Critical, 4), "\n")
    if (V > Critical)
        cat("Reject Null Hypothesis", "\n", "\n")
    else cat("Do Not Reject Null Hypothesis", "\n", "\n")
    }
    invisible(x)
}


###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   lm.circular function                                    #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: August, 01, 2003                                  #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

lm.circular <- function(y, x, order = 1, level = 0.05) {
    x <- as.circular(x)
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    y <- as.circular(y)
    y <- conversion.circular(y, units="radians")
    attr(y, "circularp") <- attr(y, "class") <- NULL

    n <- length(x)
    cy <- cos(y)
    sy <- sin(y)
    order.matrix <- t(matrix(rep(c(1:order), n), ncol = n))
    cos.x <- cos(x * order.matrix)
    sin.x <- sin(x * order.matrix)
    cos.lm <- lm(cy ~ cos.x + sin.x)
    sin.lm <- lm(sy ~ cos.x + sin.x)
    cos.fit <- cos.lm$fitted
    sin.fit <- sin.lm$fitted
    g1.sq <- t(cos.fit) %*% cos.fit
    g2.sq <- t(sin.fit) %*% sin.fit
    rho <- sqrt((g1.sq + g2.sq)/n)
    y.fitted <- atan(sin.fit, cos.fit)
    Y1 <- cy
    Y2 <- sy
    ones <- matrix(1, n, 1)
    X <- cbind(ones, cos.x, sin.x)
    W <- cbind(cos((order + 1) * x), sin((order + 1) * x))
    M <- X %*% solve(t(X) %*% X) %*% t(X)
    I <- diag(n)
    H <- t(W) %*% (I - M) %*% W
    N <- W %*% solve(H) %*% t(W)
    cc <- n - (2 * order + 1)
    N1 <- t(Y1) %*% (I - M) %*% N %*% (I - M) %*% Y1
    D1 <- t(Y1) %*% (I - M) %*% Y1
    T1 <- cc * (N1/D1)
    N2 <- t(Y2) %*% (I - M) %*% N %*% (I - M) %*% Y2
    D2 <- t(Y2) %*% (I - M) %*% Y2
    T2 <- cc * (N2/D2)
    p1 <- 1 - pchisq(T1, 2)
    p2 <- 1 - pchisq(T2, 2)
    pvalues <- cbind(p1, p2)
    circ.lm <- list()
    circ.lm$call <- match.call()
    circ.lm$rho <- rho
    circ.lm$fitted <- y.fitted %% (2 * pi)
    circ.lm$x <- cbind(x, y)
    circ.lm$residuals <- (y - y.fitted) %% (2 * pi)
    circ.lm$coefficients <- cbind(cos.lm$coefficients, sin.lm$coefficients)
    circ.lm$p.values <- pvalues
    circ.lm$A.k <- mean(cos(circ.lm$residuals))
    circ.lm$kappa <- A1inv(circ.lm$A.k)
    if (pvalues[1] > level & pvalues[2] > level)
        circ.lm$message <- paste("Higher order terms are not significant at the ", level, " level", sep = "")
    else circ.lm$message <- paste("Higher order terms are significant at the ", level, " level", sep = "")
    return(circ.lm)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   mean.circular function                                  #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.2-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

mean.circular <- function(x, na.rm=FALSE, ...) {
  if (na.rm) 
       x <- x[!is.na(x)]
   x <- as.circular(x)
   xcircularp <- attr(x, "circularp")
   unitsp <- xcircularp$units

   if (any(is.na(x))) {
       circmean <- NA
   } else {
       if (unitsp=="degrees") 
           x <- x/180*pi

       sinr <- sum(sin(x))
       cosr <- sum(cos(x))

       if (sqrt((sinr^2 + cosr^2))/length(x) > .Machine$double.eps) {
           circmean <- atan(sinr, cosr)
       } else {
           circmean <- NA
       }
   }
       
   if (unitsp=="degrees") 
       circmean <- circmean/pi*180

   attr(circmean, "circularp") <- xcircularp
   attr(circmean, "class") <- "circular"

   return(circmean)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   mle.vonmises function                                   #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: September, 22, 2003                               #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.2-3                                           #
#############################################################

mle.vonmises <- function(x, mu, kappa, bias=FALSE) {

    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")

    n <- length(x)
    sinr <- sum(sin(x))
    cosr <- sum(cos(x))
    est.mu <- FALSE 
    if (missing(mu)) {  
        mu <- atan(sinr, cosr)
        est.mu <- TRUE
    } else {
        if (units=="degrees") mu <- mu/180*pi
    }
    est.kappa <- FALSE
    if (missing(kappa)) {
        V <- mean.default(cos(x - mu))
        if (V > 0) {
            kappa <- A1inv(V)
        } else {
            kappa <- 0
        }
        if (bias == TRUE) {
            if (kappa < 2) {
                kappa <- max(kappa - 2 * (n * kappa)^-1, 0)
            } else {
                kappa <- ((n - 1)^3 * kappa)/(n^3 + n)
            }
        }
        est.kappa <- TRUE
    }

    A1temp <- A1(kappa)
    se.mu <- se.kappa <- 0
    if (est.mu) se.mu <- 1/(n*kappa*A1temp)
    if (est.kappa) se.kappa <- 1/(n*(1-A1temp/kappa-A1temp^2))
    result <- list()

    if (units=="degrees") {
        mu <- mu/pi*180
    }
    
    attr(mu, "circularp") <- xcircularp
    attr(mu, "class") <- "circular"
    
    result$call <- match.call()
    result$mu <- mu
    result$kappa <- kappa
    result$se.mu <- sqrt(se.mu)
    result$se.kappa <- sqrt(se.kappa)
    result$est.mu <- est.mu
    result$est.kappa <- est.kappa
    class(result) <- "mle.vonmises"
    return(result)
}

#############################################################
#                                                           #
#   print.mle.vonmises function                             #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-2                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.mle.vonmises <- function(x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("mu: ")
    cat(format(x$mu, digits=digits), " (", format(x$se.mu, digits=digits), ")\n")
    cat("\n")
    cat("kappa: ")    
    cat(format(x$kappa, digits=digits), " (", format(x$se.kappa, digits=digits), ")\n")
    cat("\n")    
    if (!x$est.mu) cat("mu is known\n")
    if (!x$est.kappa) cat("kappa is known\n")
    invisible(x)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   mle.vonmises.bootstrap.ci function                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: August, 1, 2003                                   #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-1                                           #
#############################################################

mle.vonmises.bootstrap.ci <- function(x, mu, bias = FALSE, alpha = 0.05, reps = 1000) {

  if (require(boot)) {

      x <- as.circular(x)
      xcircularp <- circularp(x)
      units <- xcircularp$units
      x <- conversion.circular(x, units="radians")

      if (missing(mu)) {
          sinr <- sum(sin(x))
          cosr <- sum(cos(x))
          mu <- atan(sinr, cosr)
      } else {
          attr(mu, "circularp") <- xcircularp
          attr(mu, "class") <- "circular"
          mu <- conversion.circular(mu, units="radians")
      }
    
      mle.vonmises.mu <- function(x, i) {
          sinr <- sum(sin(x[i]))
          cosr <- sum(cos(x[i]))
          mu <- atan(sinr, cosr)
          return(mu)
      }

      mle.vonmises.kappa <- function(x, i, mu, bias) {
          n <- length(x[i])
          V <- mean(cos(x[i] - mu))
          if (V > 0) {
              kappa <- A1inv(V)
          } else {
              kappa <- 0
          }
          if (bias == TRUE) {
              if (kappa < 2) {
                  kappa <- max(kappa - 2 * (n * kappa)^-1, 0)
              } else {
                  kappa <- ((n - 1)^3 * kappa)/(n^3 + n)
              }
          }
          return(kappa)
      }
      
      mean.bs <- boot(data = x, statistic = mle.vonmises.mu, R = reps, stype="i")

      mean.reps <- mean.bs$t
      mean.reps <- sort(mean.reps %% (2 * pi))
      B <- reps
      spacings <- c(diff(mean.reps), mean.reps[1] - mean.reps[B] + 2 * pi)
      max.spacing <- (1:B)[spacings == max(spacings)]
      off.set <- 2 * pi - mean.reps[max.spacing + 1]
      if (max.spacing != B)
      mean.reps2 <- mean.reps + off.set
      else mean.reps2 <- mean.reps
      mean.reps2 <- sort(mean.reps2 %% (2 * pi))
      mean.ci <- quantile(mean.reps2, c(alpha/2, 1 - alpha/2))
      if (max.spacing != B)
      mean.ci <- mean.ci - off.set
    
      
      kappa.bs <- boot(data = x, statistic = mle.vonmises.kappa, R = reps, stype="i", mu=mu, bias = bias)
      kappa.reps <- kappa.bs$t

      kappa.ci <- quantile(kappa.reps, c(alpha/2, 1 - alpha/2))

      if (units=="degrees") {
          mean.reps <- mean.reps/pi*180
          mean.ci <- mean.ci/pi*180
      }
      attr(mean.reps, "circularp") <- xcircularp
      attr(mean.reps, "class") <- "circular"
      attr(mean.ci, "circularp") <- xcircularp
      attr(mean.ci, "class") <- "circular"

      result <- list()
      result$call <- match.call()
      result$mu.ci <- mean.ci
      result$mu <- c(mean.reps)
      result$kappa.ci <- kappa.ci
      result$kappa <- c(kappa.reps)
      result$alpha <- alpha
      class(result) <- "mle.vonmises.bootstrap.ci"
      return(result)

   } else {
       stop("To use this function you have to install the package 'boot' \n")
   }

}


#############################################################
#                                                           #
#   print.mle.vonmises.bootstrap.ci function                #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: September, 17, 2003                               #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

print.mle.vonmises.bootstrap.ci <- function(x, ...) {
    cat("Bootstrap Confidence Intervals for Mean Direction and Concentration", "\n")
    cat("Confidence Level:  ", round(100 * (1 - x$alpha),2), "%", "\n")
    cat("Mean Direction:           ", "Low =", round(x$mu.ci[1], 2), "  High =", round(x$mu.ci[2], 2), "\n")
    cat("Concentration Parameter:  ", "Low =", round(x$kappa.ci[1], 2), "  High =", round(x$kappa.ci[2], 2), "\n")

}
        

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   mle.wrappedcauchy function                              #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 31, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-1                                           #
#############################################################

mle.wrappedcauchy <- function(x, mu, rho, tol = 1e-015, max.iter = 100) {
    if (length(tol)==1) tol <- rep(tol, 2)
    if (length(tol) > 2) stop("'tol' must have less than 2 elements")
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    if (missing(mu)) {
        mu <- mean.circular(x)
    } else {
        if (units=="degrees") mu <- mu/180*pi
    }
    
    if (missing(rho)) rho <- rho.circular(x)
    if (rho < 0 | rho > 1) stop("'rho' must be between 0 and 1")
    
    mu1.old <- (2 * rho * cos(mu))/(1 + rho^2)
    mu2.old <- (2 * rho * sin(mu))/(1 + rho^2)
    w.old <- 1/(1 - mu1.old * cos(x) - mu2.old * sin(x))
    flag <- TRUE
    iter <- 0 
    while (flag & iter <= max.iter) {
           iter <- iter + 1
       mu1.new <- sum(w.old * cos(x))/sum(w.old)
       mu2.new <- sum(w.old * sin(x))/sum(w.old)
       diff1 <- abs(mu1.new - mu1.old)
       diff2 <- abs(mu2.new - mu2.old)
       if ((diff1 < tol[1]) && (diff2 < tol[2]))
           flag <- FALSE
       else {
           mu1.old <- mu1.new
           mu2.old <- mu2.new
           w.old <- 1/(1 - mu1.old * cos(x) - mu2.old * sin(x))
       }
    }
    mu.const <- sqrt(mu1.new^2 + mu2.new^2)
    rho <- (1 - sqrt(1 - mu.const^2))/mu.const
    mu <- atan(mu2.new, mu1.new) %% (2 * pi)
    if (units=="degrees") {
        mu <- mu/pi*180
    }
    
    attr(mu, "circularp") <- xcircularp
    attr(mu, "class") <- "circular"

    result <- list()
    result$call <- match.call()
    result$mu <- mu 
    result$rho <- rho
    result$convergence <- TRUE
    if (iter > max.iter) {
        result$convergence <- FALSE
    }
    class(result) <- "mle.wrappedcauchy"
    return(result)
} 

#############################################################
#                                                           #
#   print.mle.wrappednormal function                    #
#   Author: Claudio Agostinelli                         #
#   E-mail: claudio@unive.it                            #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                        #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################

print.mle.wrappedcauchy <- function(x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("mu: ")
    cat(format(x$mu, digits=digits), "\n")
    cat("\n")
    cat("rho: ")    
    cat(format(x$rho, digits=digits), "\n")
    if (!x$convergence) cat("\n The convergence is not achieved after the prescribed number of iterations \n")

    invisible(x)
}
#############################################################
#                                                           #
#   mle.wrappednormal function                              #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 31, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-2                                           #
#############################################################

mle.wrappednormal <- function(x, mu, rho, sd, K, tol=1e-5, min.sd=1e-3, min.k=10, max.iter=100, verbose=FALSE) {

    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
  
    n <- length(x)
    sinr <- sum(sin(x))
    cosr <- sum(cos(x))

    est.mu <- FALSE 
    if (missing(mu)) {  
        mu <- atan(sinr, cosr)
        est.mu <- TRUE
    }
    est.rho <- FALSE
    if (missing(sd)) {
        if (missing(rho)) {
            sd <- sqrt(-2*log(sqrt(sinr^2 + cosr^2)/n))
            if (is.na(sd) || sd < min.sd) sd <- min.sd
            est.rho <- TRUE
        } else {
            sd <- sqrt(-2*log(rho))
        }
    }
     
    xdiff <- 1+tol
    iter <- 0
    if (missing(K)) {
        range <- max(mu, x) - min(mu, x)
        K <- (range+6*sd)%/%(2*pi)+1
        K <- max(min.k, K)
    }
    
    while (xdiff > tol & iter <= max.iter) {
           iter <- iter + 1
           mu.old <- mu
           sd.old <- sd

           z <- .Fortran("mlewrpno",
                    as.double(x),
                    as.double(mu),
                    as.double(sd),
                    as.integer(n),
                    as.integer(K),
                    as.integer(est.mu),
                    as.integer(est.rho),
                    w=double(n),
                    wk=double(n),
                    wm=double(n),
                    PACKAGE="circular"
           )
           w <- z$w
           wk <- z$wk
           wm <- z$wm
           
           if (est.mu) {
               mu <- sum(x)/n
               if (any(wk!=0)) {
                   mu <- mu + 2*pi*mean(wk[wk!=0]/w[wk!=0])
               }
           }
           if (est.rho) {
               if  (any(wm!=0)) {
                    sd <- sqrt(sum(wm[wm!=0]/w[wm!=0])/n)
               } else {
                    sd <- min.sd
               }
           }

           if (verbose) {
               cat("mu: ", mu, "\n")
               cat("rho: ", exp(-sd^2/2), "\n")              
               cat("sd: ", sd, "\n")
           }
           xdiff <- max(abs(mu - mu.old), abs(sd - sd.old))
    }

    rho <- exp(-sd^2/2)

    if (units=="degrees") {
        mu <- mu/pi*180
        sd <- sd/pi*180
    }
    
    attr(mu, "circularp") <- xcircularp
    attr(mu, "class") <- "circular"
    
    result <- list()
    
    result$call <- match.call()
    result$mu <- mu
    result$rho <- rho
    result$sd <- sd
    result$est.mu <- est.mu
    result$est.rho <- est.rho
    result$convergence <- TRUE
    if (iter > max.iter) {
        result$convergence <- FALSE
    }
    class(result) <- "mle.wrappednormal"
    return(result)
}

#############################################################
#                                                           #
#	print.mle.wrappednormal function                    #
#	Author: Claudio Agostinelli                         #
#	E-mail: claudio@unive.it                            #
#	Date: November, 19, 2003                                #
#	Version: 0.1-2                                      #
#                                                           #
#	Copyright (C) 2003 Claudio Agostinelli              #
#                                                           #
#############################################################

print.mle.wrappednormal <- function(x, digits = max(3, getOption("digits") - 3), ...) {
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("mu: ")
    cat(format(x$mu, digits=digits), "\n")
    cat("\n")
    cat("rho: ")    
    cat(format(x$rho, digits=digits), "\n")
    cat("\n")
    cat("sd: ")       
    cat(format(x$sd, digits=digits), "\n")
    cat("\n")   
    if (!x$est.mu) cat("mu is known\n")
    if (!x$est.rho) {
        cat("rho and sd are known\n")
    }
    if (!x$convergence) cat("\nThe convergence is not achieved after the prescribed number of iterations \n")
    invisible(x)
}

#############################################################
#                                                           #
#   axis.circular function                                  #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: October, 07, 2003                                 #
#   Version: 0.3-2                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################
 
axis.circular <- function(at, labels,  units = c("radians", "degrees"), template=c("none", "geographics"), modulo = c("asis", "2pi", "pi"), zero=0, rotation = c("counter", "clock"), tick=TRUE, lty, lwd, cex, col, font, tcl=0.025, tcl.text=0.125, digits=2) {

  units <- match.arg(units)
  template <- match.arg(template)
  modulo <- match.arg(modulo)
  rotation <- match.arg(rotation)

  if (missing(at)) {
      if (units=="radians") {
          at <- c(0, pi/2, pi, 3/2*pi)
      } else {
          at <- c(0, 90, 180, 270)
      }   
  }

  at <- as.circular(at, units=units, template=template, modulo=modulo, zero=zero, rotation=rotation)
  atcircularp <- circularp(at)
  zero <- atcircularp$zero
  rotation <- atcircularp$rotation
  atasis <- at
  attr(atasis, "circularp") <- attr(atasis, "class") <- NULL
  at <- conversion.circular(at, units="radians")
  attr(at, "circularp") <- attr(at, "class") <- NULL

  attext <- round(at/pi, digits=digits)
  
  if (missing(cex)) cex <- par("cex.axis")
  if (missing(col)) col <- par("col.axis")
  if (missing(font)) font <- par("font.axis")
  if (missing(lty)) lty <- par("lty")
  if (missing(lwd)) lwd <- par("lwd")

  if (missing(labels) & all(at==c(0, pi/2, pi, 3/2*pi))) { 
      if (template=="geographics") {
          labels <- c("N", "E", "S", "W")         
      } else {
          if (units=="radians") {
              labels <- c("0", expression(frac(pi,2)), expression(pi), expression(frac(3*pi,2)))
          } else {
              labels <- c("0", "90", "180", "270")      
          }
      }   
  } else {
      if (missing(labels)) {
          if (units!="radians") {
              labels <- as.character(round(atasis, digits=digits))
          }
      }
  }

if (!missing(labels) && length(at)!=length(labels)) stop("'at' and 'labels' must have the same length")
  
  if (rotation=="clock") at <- -at 
  at <- at + zero
        
  r <- 1+tcl*c(-1/2,1/2)
  r.l <- 1-tcl.text 
  z <- cos(at)
  y <- sin(at)
  
  for (i in 1:length(at)) {
       if (tick) {
           lines.default(z[i]*r, y[i]*r, col=col, lty=lty, lwd=lwd)
       }
       if (missing(labels) & units=="radians") {
              labeltext <- substitute(at*pi, list(at=attext[i]))
       } else {
              labeltext <- labels[i]
       }
       text.default(z[i]*r.l, y[i]*r.l, labeltext, cex=cex, col=col)    
  }
     
  text(0, 0, "+", cex=1)
}

#############################################################
#                                                           #
#   ticks.circular function                                 #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: September, 22, 2003                               #
#   Version: 0.2-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################
 
ticks.circular <- function(x, template=c("none", "geographics"), zero, rotation, tcl=0.025, col, ...) {

  template <- match.arg(template)
  if (missing(col)) col <- par("col")

  x <- as.circular(x, template=template)
  xcircularp <- attr(x, "circularp")
  xzero <- xcircularp$zero
  xrotation <- xcircularp$rotation
  attr(x, "circularp") <- attr(x, "class") <- NULL

  if (missing(zero)) {
      if (template=="geographics") {
          zero <- pi/2
      } else {
          zero <- xzero
      }
  }
  
  if (missing(rotation)) {
      if (template=="geographics") {
          rotation <- "clock"
      } else {
          rotation <- xrotation
      }
  }
  
  if (rotation=="clock") x <- -x 
  x <- x + zero
        
  r <- 1+tcl*c(-1/2,1/2)
  z <- cos(x)
  y <- sin(x)

  for (i in 1:length(x)) {
       lines.default(z[i]*r, y[i]*r, col=col, ...)
  }
}

#############################################################
#                                                           #
#   plot.circular function                                  #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 18, 2003                                #
#   Version: 0.2-5                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################
 
plot.circular <- function(x, pch = 16, cex = 1, stack = FALSE, axes = TRUE, sep = 0.025, shrink = 1, bins, ticks = FALSE, tcl=0.025, tcl.text=0.125, col, tol = 0.04, uin, xlim=c(-1, 1), ylim=c(-1, 1), main=NULL, digits=2, ...) {

    if (is.matrix(x) | is.data.frame(x)) {
        nseries <- ncol(x)
    } else {
        nseries <- 1
    }
    xx <- as.data.frame(x)
  
    xcircularp <- attr(as.circular(xx[,1]), "circularp")
    type <- xcircularp$type
    units <- xcircularp$units
    template <- xcircularp$template
    modulo <- xcircularp$modulo
    zero <- xcircularp$zero
    rotation <- xcircularp$rotation

    xlim <- shrink * xlim
    ylim <- shrink * ylim
    midx <- 0.5 * (xlim[2] + xlim[1])
    xlim <- midx + (1 + tol) * 0.5 * c(-1, 1) * (xlim[2] - xlim[1])
    midy <- 0.5 * (ylim[2] + ylim[1])
    ylim <- midy + (1 + tol) * 0.5 * c(-1, 1) * (ylim[2] - ylim[1])
    oldpin <- par("pin")
    xuin <- oxuin <- oldpin[1]/diff(xlim)
    yuin <- oyuin <- oldpin[2]/diff(ylim)
    if (missing(uin)) {
        if (yuin > xuin) yuin <- xuin
        else xuin <- yuin
    } else {
        if (length(uin) == 1) uin <- uin * c(1, 1)
        if (any(c(xuin, yuin) < uin)) stop("uin is too large to fit plot in")
        xuin <- uin[1]; yuin <- uin[2]
    }    
    xlim <- midx + oxuin/xuin * c(-1, 1) * diff(xlim) * 0.5
    ylim <- midy + oyuin/yuin * c(-1, 1) * diff(ylim) * 0.5
    plot(cos(seq(0, 2 * pi, length = 1000)), sin(seq(0, 2 * pi, length = 1000)), axes = FALSE, xlab = "", ylab = "", main = main, type = "l", xlim=xlim, ylim=ylim, xaxs="i", yaxs="i")
    
    if (missing(bins)) {
         bins <- NROW(x)
    } else {
         bins <- round(bins)
         if (bins<=0) stop("bins must be non negative")
    }

    if (axes) {
        axis.circular(units = units, template=template, zero=zero, rotation=rotation, digits=digits, cex=cex, tcl=tcl, tcl.text=tcl.text)
    }

    if (missing(col)) {
    col <- seq(nseries)
    } else {
    if (length(col)!=nseries) {
        col <- rep(col, nseries)[1:nseries]
    }
    }
    pch <- rep(pch, nseries, length.out=nseries)
    
    if (!is.logical(ticks)) stop("ticks must be logical")

    arc <- (2 * pi)/bins
    pos.bins <- ((1:nseries)-1/2)*arc/nseries-arc/2
            
    if (ticks) {
        at <- (0:bins)/bins*2*pi
        if (rotation=="clock") at <- -at
        at <- at + zero

        ticks.circular(circular(x=at, type="angles", units="radians", modulo="asis", zero=zero, rotation=rotation), tcl=tcl)
    }

    for (iseries in 1:nseries) {
      
     x <- xx[,iseries]
         x <- as.circular(x)    
         x <- conversion.circular(x, units="radians")
     n <- length(x)

     if (!stack) {
             if (rotation=="clock") x <- -x
             x <- x + zero 
             z <- cos(x)
             y <- sin(x)
         r <- 1+(iseries-1)*sep*shrink
         points.default(z*r, y*r, cex=cex, pch=pch[iseries], col = col[iseries], ...)
     } else {
         bins.count <- c(1:bins)
         xmod <- x %% (2 * pi)
         for (i in 1:bins) {
          bins.count[i] <- sum(xmod <= i * arc & xmod > (i - 1) * arc)
         }
         mids <- seq(arc/2, 2 * pi - pi/bins, length = bins) + pos.bins[iseries]
             if (rotation=="clock") mids <- -mids
             mids <- mids + zero
 
         index <- cex*sep
         for (i in 1:bins) {
          if (bins.count[i] != 0) {
              for (j in 0:(bins.count[i] - 1)) {
               r <- 1 + j * index
               z <- r * cos(mids[i])
               y <- r * sin(mids[i])
               points.default(z, y, cex=cex, pch=pch[iseries], col=col[iseries], ...)
              }
           }
          }
     }
    }
return(invisible(list(zero=zero, rotation=rotation, next.points=nseries*sep)))
}

#############################################################
#                                                           #
#   points.circular function                                #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: September, 22, 2003                               #
#   Version: 0.1-2                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################
 
points.circular <- function(x, pch = 16, cex = 1, stack = FALSE, sep = 0.025, shrink=1, bins, col, next.points, plot.info, zero, rotation, ...) {

    if (is.matrix(x) | is.data.frame(x)) {
        nseries <- ncol(x)
    } else {
        nseries <- 1
    }
    xx <- as.data.frame(x)
  
    xcircularp <- attr(as.circular(xx[,1]), "circularp")
    type <- xcircularp$type
    modulo <- xcircularp$modulo
    if (missing(plot.info)) {
        if (missing(zero)) zero <- xcircularp$zero
        if (missing(rotation)) rotation <- xcircularp$rotation
        if (missing(next.points)) next.points <- 0
    } else {
        zero <- plot.info$zero
        rotation <- plot.info$rotation
        if (missing(next.points))
            next.points <- plot.info$next.points
    }
        
    x <- conversion.circular(x, units="radians")

    if (missing(bins)) {
    bins <- NROW(x)
    } else {
    bins <- round(bins)
    if (bins<=0) stop("bins must be non negative")
    }

    if (is.matrix(x) || is.data.frame(x)) {
        nseries <- ncol(x)
    } else {
        nseries <- 1
    }
 
    if (missing(col)) {
    col <- seq(nseries)
    } else {
    if (length(col)!=nseries) {
        col <- rep(col, nseries)[1:nseries]
    }
    }
    pch <- rep(pch, nseries, length.out=nseries)

    arc <- (2 * pi)/bins
    pos.bins <- ((1:nseries)-1/2)*arc/nseries-arc/2
            
    for (iseries in 1:nseries) {
     x <- xx[,iseries] 
         x <- as.circular(x)
         x <- conversion.circular(x, units="radians")
         n <- length(x)

     if (!stack) {
             if (rotation=="clock") x <- -x
             x <- x + zero
             z <- cos(x)
             y <- sin(x)
         r <- 1+next.points+(iseries-1)*sep*shrink
         points.default(z*r, y*r, cex=cex, pch=pch[iseries], col = col[iseries], ...)
     } else {
         bins.count <- c(1:bins)

         for (i in 1:bins) {
          bins.count[i] <- sum(x <= i * arc & x > (i - 1) * arc)
         }
         mids <- seq(arc/2, 2 * pi - pi/bins, length = bins) + pos.bins[iseries]
             if (rotation=="clock") mids <- -mids
             mids <- mids + zero
 
         index <- cex*sep
         for (i in 1:bins) {
          if (bins.count[i] != 0) {
              for (j in 0:(bins.count[i] - 1)) {
               r <- 1 + j * index
               z <- r * cos(mids[i])
               y <- r * sin(mids[i])
               points.default(z, y, cex=cex, pch=pch[iseries], col=col[iseries], ...)
              }
           }
          }
     }
    }
return(invisible(list(zero=zero, rotation=rotation, next.points=next.points+nseries*sep)))
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   plot.edf function                                       #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 29, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

plot.edf <- function(x, type = "s", xlim = c(0, 2 * pi), ylim = c(0, 1), ...) {
    x <- as.circular(x)
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    x <- x %% (2 * pi)
    x <- sort(x)
    n <- length(x)
    plot.default(c(0, x, 2 * pi), c(0, seq(1:n)/n, 1), type=type, xlim=xlim, ylim=ylim, ...)
}

#############################################################
#                                                           #
#   lines.edf function                                      #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: August, 01, 2003                                  #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

lines.edf <- function(x, type = "s", ...) {
    x <- as.circular(x)
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    x <- x %% (2 * pi)
    x <- sort(x)
    n <- length(x)
    lines.default(c(0, x, 2 * pi), c(0, seq(1:n)/n, 1), type=type, ...)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   pp.plot function                                        #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

pp.plot <- function(x, ref.line = TRUE, tol=1e-20,  xlab = "von Mises Distribution", ylab = "Empirical Distribution", ...) {
        x <- as.circular(x)
        xcircularp <- circularp(x)
        units <- xcircularp$units
        x <- conversion.circular(x, units="radians")

        res <- mle.vonmises(x)
        mu <- res$mu
    kappa <- res$kappa

    n <- length(x)
    x <- sort(x %% (2 * pi))
    z <- (1:n)/(n + 1)
    
    y <- pvonmises(q=x, mu=mu, kappa=kappa, tol=tol)
    
    plot.default(z, y, xlab=xlab, ylab=ylab, ...)
    if (ref.line)
        abline(0, 1)
        if (units=="degrees") mu <- mu/pi*180
        attr(mu, "circularp") <- xcircularp
        attr(mu, "class") <- "circular"
    invisible(list(mu=mu, kappa=kappa))
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

###############################################################
#                                                             #
#       R port: Claudio Agostinelli  <claudio@unive.it>       #
#                                                             #
#       Date: January, 14, 2003                               #
#       Version: 0.1-6                                        #
#                                                             #
###############################################################

rad <- function(x) {
    (x * pi)/180
}
###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################


#############################################################
#                                                           #
#   range.circular function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: December, 23, 2003                                #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-1                                           #
#############################################################

range.circular <- function(x, test = FALSE, na.rm=FALSE, finite=FALSE, ...) {
       x <- as.circular(x)
       units <- circularp(x)$units
       x <- conversion.circular(x, units="radians")
       if (finite) 
           x <- x[is.finite(x)]
        else if (na.rm) 
           x <- x[!is.na(x)]
        if (length(x)) 
            c(min(x), max(x))
        else c(NA, NA)
    x <- sort(x %% (2*pi))
    n <- length(x)
    spacings <- c(diff(x), x[1] - x[n] + 2*pi)
    range <- 2*pi - max(spacings)
        if (units=="degrees") 
            rangenew <- range/pi*180
        else
            rangenew <- range
        result <- rangenew
    if(test == TRUE) {
        stop <- floor(1/(1 - range/(2*pi)))
        index <- c(1:stop)
        sequence <- ((-1)^(index - 1)) * exp(log(gamma(n + 1)) - log(gamma(index + 1)) - log(gamma(n - index + 1))) * (1 - index * (1 - range/(2 * pi)))^(n - 1)
        p.value <- sum(sequence)
        result <- list(range=rangenew, p.value=p.value)
    }
    result
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rao.spacing.test function                               #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 26, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

rao.spacing.test <- function(x, alpha = 0) {
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="degrees")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    if (!any(c(0, 0.01, 0.025, 0.05, 0.1, 0.15)==alpha)) stop("'alpha' must be one of the following values: 0, 0.01, 0.025, 0.05, 0.1, 0.15")
    x <- sort(x %% 360)
    n <- length(x)
 if (n < 4) {
     warning("Sample size too small")
     U <- NA
 } else {
    spacings <- c(diff(x), x[1] - x[n] + 360)
    U <- 1/2 * sum(abs(spacings - 360/n))
 }
    result <- list(statistic=U, alpha=alpha, n=n)
    class(result) <- "rao.spacing.test"
    return(result)
}

#############################################################
#                                                           #
#   print.rao.spacing.test function                         #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.rao.spacing.test <- function(x, digits=4, ...) {
    U <- x$statistic
    alpha <- x$alpha
    n <- x$n
    data(rao.table)
    if (n <= 30)
    table.row <- n - 3
    else if (n <= 32)
         table.row <- 27
    else if (n <= 37)
         table.row <- 28
    else if (n <= 42)
         table.row <- 29
    else if (n <= 47)
         table.row <- 30
    else if (n <= 62)
         table.row <- 31
    else if (n <= 87)
         table.row <- 32
    else if (n <= 125)
         table.row <- 33
    else if (n <= 175)
         table.row <- 34
    else if (n <= 250)
         table.row <- 35
    else if (n <= 350)
         table.row <- 36
    else if (n <= 450)
         table.row <- 37
    else if (n <= 550)
         table.row <- 38
    else if (n <= 650)
         table.row <- 39
    else if (n <= 750)
         table.row <- 40
    else if (n <= 850)
         table.row <- 41
    else if (n <= 950)
         table.row <- 42
    else table.row <- 43
        
    cat("\n")
    cat("       Rao's Spacing Test of Uniformity", "\n", "\n")
    cat("Test Statistic =", round(U, digits=digits), "\n")
   
    if (alpha == 0) {
        if (U > rao.table[table.row, 1])
        cat("P-value < 0.001", "\n", "\n")
        else if (U > rao.table[table.row, 2])
             cat("0.001 < P-value < 0.01", "\n", "\n")
        else if (U > rao.table[table.row, 3])
             cat("0.01 < P-value < 0.05", "\n", "\n")
        else if (U > rao.table[table.row, 4])
             cat("0.05 < P-value < 0.10", "\n", "\n")
        else cat("P-value > 0.10", "\n", "\n")
    } else {
        table.col <- (1:4)[alpha == c(0.001, 0.01, 0.05, 0.1)]
        critical <- rao.table[table.row, table.col]
        cat("Level", alpha, "critical value =", critical, "\n")
        if (U > critical)
        cat("Reject null hypothesis of uniformity \n\n")
        else
            cat("Do not reject null hypothesis of uniformity \n\n")
  }
invisible(x)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rao.test function                                       #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 25, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

rao.test <- function(..., alpha = 0) {
        y <- list(...)
        x <- list()
        for (i in 1:length(y)) {
             if (is.data.frame(y[[i]])) {
                 x <- c(x, as.list(y[[i]]))
             } else if (is.matrix(y[[i]])) {
                        for (j in 1:ncol(y[[i]])) {
                             x <- c(x, list(y[[i]][,j]))
                        }
             } else if (is.list(y[[i]])) {
                        x <- c(x, y[[i]])
             } else {
                 x <- c(x, list(y[[i]]))
             }
        }
        if (length(x)<2) stop("There must be at least two samples")
        for (i in 1:length(x)) {
             x[[i]] <- as.circular(x[[i]])
             x[[i]] <- conversion.circular(x[[i]], units="radians")
             attr(x[[i]], "circularp") <- attr(x[[i]], "class") <- NULL
        }
        if (!any(c(0, 0.01, 0.025, 0.05, 0.1, 0.15)==alpha)) stop("'alpha' must be one of the following values: 0, 0.01, 0.025, 0.05, 0.1, 0.15")
    n <- unlist(lapply(x, length))
    k <- length(x)
    c.data <- lapply(x, cos)
    s.data <- lapply(x, sin)
    x <- unlist(lapply(c.data, mean))
    y <- unlist(lapply(s.data, mean))
    s.co <- unlist(lapply(c.data, var))
    s.ss <- unlist(lapply(s.data, var))
    s.cs <- c(1:k)
    for(i in 1:k) {
        s.cs[i] <- var(c.data[[i]], s.data[[i]])
    }
    s.polar <- 1/n * (s.ss/x^2 + (y^2 * s.co)/x^4 - (2 * y * s.cs)/x^3)
    tan <- y/x
    H.polar <- sum(tan^2/s.polar) - (sum(tan/s.polar))^2/sum(1/s.polar)
    U <- x^2 + y^2
    s.disp <- 4/n * (x^2 * s.co + y^2 * s.ss + 2 * x * y * s.cs)
    H.disp <- sum(U^2/s.disp) - (sum(U/s.disp))^2/sum(1/s.disp)
        
    result <- list(statistic=c(H.polar, H.disp), df=k-1, p.value=c((1 - pchisq(H.polar, k - 1)), (1 - pchisq(H.disp, k - 1))), alpha=alpha)
    class(result) <- "rao.test"
    return(result)
}

#############################################################
#                                                           #
#   print.rao.test function                                 #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.rao.test <- function(x, digits=4, ...) {
        statistic <- x$statistic
        p.value <- x$p.value
        alpha <- x$alpha
        df <- x$df
    cat("\n")
    cat("Rao's Tests for Homogeneity", "\n")
    if(alpha == 0) {
        cat("\n")
        cat("       Test for Equality of Polar Vectors:", "\n", "\n")
        cat("Test Statistic =", round(statistic[1], digits=digits), "\n")
        cat("Degrees of Freedom =", df, "\n")
        cat("P-value of test =", round(p.value[1], digits=digits), "\n", "\n")
        cat("       Test for Equality of Dispersions:", "\n", "\n")
        cat("Test Statistic =", round(statistic[2], digits=digits), "\n")
        cat("Degrees of Freedom =", df, "\n")
        cat("P-value of test =", round(p.value[2], digits=digits), "\n", "\n")
    } else {
        cat("\n")
        cat("       Test for Equality of Polar Vectors:", "\n", "\n")
        cat("Test Statistic =", round(statistic[1], digits=digits), "\n")
        cat("Degrees of Freedom =", df, "\n")
        cat("Level", alpha, "critical value =", round(qchisq(1 - alpha, df), digits=digits), "\n")
        if (statistic[1] > qchisq(1 - alpha, df)) {
            cat("Reject null hypothesis of equal polar vectors", "\n", "\n")
        } else { 
                    cat("Do not reject null hypothesis of equal polar vectors", "\n", "\n")
                }
        cat("       Test for Equality of Dispersions:", "\n", "\n")
        cat("Test Statistic =", round(statistic[2], digits=digits), "\n")
        cat("Degrees of Freedom =", df, "\n")
        cat("Level", alpha, "critical value =", round(qchisq(1 - alpha, df), digits=digits), "\n")
        if (statistic[2] > qchisq(1 - alpha, df)) {
            cat("Reject null hypothesis of equal dispersions", "\n", "\n")
        } else {
                    cat("Do not reject null hypothesis of equal dispersions", "\n", "\n")
                }
    }
        invisible(x)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rayleigh.test function                                  #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 25, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

rayleigh.test <- function(x, mu) {
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    n <- length(x)
    if (missing(mu)) {
        ss <- sum(sin(x))
        cc <- sum(cos(x))
        rbar <- (sqrt(ss^2 + cc^2))/n
        z <- (n * rbar^2)
        p.value <- exp( - z)
        if (n < 50)
        temp <- (1 + (2 * z - z^2)/(4 * n) - (24 * z - 132 * z^2 + 76 * z^3 - 9 * z^4)/(288 * n^2))
        else temp <- 1
        result <- list(statistic = rbar, p.value = p.value * temp, mu=NA)
    } else {
    r0.bar <- (sum(cos(x - mu)))/n
    z0 <- sqrt(2 * n) * r0.bar
    pz <- pnorm(z0)
    fz <- dnorm(z0)
    p.value <- 1 - pz + fz * ((3 * z0 - z0^3)/(16 * n) + (15 * z0 + 305 * z0^3 - 125 * z0^5 + 9 * z0^7)/(4608 * n^2))
    result <- list(statistic = r0.bar, p.value = p.value, mu=mu)
    }
    class(result) <- "rayleigh.test"
    return(result)
}

#############################################################
#                                                           #
#   print.rayleigh.test function                            #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.rayleigh.test <- function(x, digits=4, ...) {
    rbar <- x$statistic
    p.value <- x$p.value
    mu <- x$mu
    cat("\n", "      Rayleigh Test of Uniformity \n")
    if (is.na(mu)) {
        cat("       General Unimodal Alternative \n\n")
    } else {
        cat("       Alternative with Specified Mean Direction: ", mu, "\n\n")
    }
    cat("Test Statistic: ", round(rbar, digits=digits), "\n")
    cat("P-value: ", round(p.value, digits=digits), "\n\n")
    invisible(x)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rho.circular function                                   #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 24, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

rho.circular <- function(x, na.rm=FALSE) {
    if (na.rm) 
        x <- x[!is.na(x)]
    if (any(is.na(x))) return(NA)

    n <- length(x)
    x <- as.circular(x)
    x <- conversion.circular(x, units="radians")
    
    sinr <- sum(sin(x))
    cosr <- sum(cos(x))
    sqrt(sinr^2 + cosr^2)/n
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rose.diag function                                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: September, 22, 2003                               #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-2                                           #
#############################################################

rose.diag <- function(x, pch = 16, axes = TRUE, shrink = 1, bins, ticks = TRUE, tcl=0.025, col, tol = 0.04, uin, xlim=c(-1, 1), ylim=c(-1, 1), prop = 1, main=NULL, ...) {
  
    if (is.matrix(x) | is.data.frame(x)) {
        nseries <- ncol(x)
    } else {
        nseries <- 1
    }
    xx <- as.data.frame(x)
  
    xcircularp <- attr(as.circular(xx[,1]), "circularp")
    type <- xcircularp$type
    units <- xcircularp$units
    template <- xcircularp$template
    modulo <- xcircularp$modulo
    zero <- xcircularp$zero
    rotation <- xcircularp$rotation

    xlim <- shrink * xlim
    ylim <- shrink * ylim
    midx <- 0.5 * (xlim[2] + xlim[1])
    xlim <- midx + (1 + tol) * 0.5 * c(-1, 1) * (xlim[2] - xlim[1])
    midy <- 0.5 * (ylim[2] + ylim[1])
    ylim <- midy + (1 + tol) * 0.5 * c(-1, 1) * (ylim[2] - ylim[1])
    oldpin <- par("pin")
    xuin <- oxuin <- oldpin[1]/diff(xlim)
    yuin <- oyuin <- oldpin[2]/diff(ylim)
    if (missing(uin)) {
        if (yuin > xuin) yuin <- xuin
        else xuin <- yuin
    } else {
        if (length(uin) == 1) uin <- uin * c(1, 1)
        if (any(c(xuin, yuin) < uin)) stop("uin is too large to fit plot in")
        xuin <- uin[1]; yuin <- uin[2]
    }    
    xlim <- midx + oxuin/xuin * c(-1, 1) * diff(xlim) * 0.5
    ylim <- midy + oyuin/yuin * c(-1, 1) * diff(ylim) * 0.5
    plot(cos(seq(0, 2 * pi, length = 1000)), sin(seq(0, 2 * pi, length = 1000)), axes = FALSE, xlab = "", ylab = "", main = main, type = "l", xlim=xlim, ylim=ylim, xaxs="i", yaxs="i")
    
    if (missing(bins)) {
    bins <- NROW(x)
    } else {
    bins <- round(bins)
    if (bins<=0) stop("bins must be non negative")
    }
    
    if (axes) {
    axis.circular(units = units, template=template, modulo = modulo, zero=zero, rotation=rotation)
    }

    if (missing(col)) {
    col <- seq(nseries)
    } else {
    if (length(col)!=nseries) {
        col <- rep(col, nseries)[1:nseries]
    }
    }
    pch <- rep(pch, nseries, length.out=nseries)
    
    if (!is.logical(ticks)) stop("ticks must be logical")

    arc <- (2 * pi)/bins
    pos.bins <- ((1:nseries)-1/2)*arc/nseries-arc/2
            
    if (ticks) {
        at <- (0:bins)/bins*2*pi
        if (rotation=="clock") at <- -at
        at <- at + zero

        ticks.circular(circular(x=at, type="angles", units="radians", modulo="asis", zero=zero, rotation=rotation), tcl=tcl)
    }

    for (iseries in 1:nseries) {
      
    x <- xx[,iseries]
        x <- as.circular(x)    
        x <- conversion.circular(x, units="radians")
        if (rotation=="clock") x <- -x
        x <- x + zero 
        x <- x %% (2 * pi)
        n <- length(x)
    freq <- c(1:bins)
    arc <- (2 * pi)/bins
    for(i in 1:bins) {
        freq[i] <- sum(x <= i * arc & x > (i - 1) * arc)
    }
    rel.freq <- freq/n
    radius <- sqrt(rel.freq) * prop
    sector <- seq(0, 2 * pi - (2 * pi)/bins, length = bins)
    mids <- seq(arc/2, 2 * pi - pi/bins, length = bins)
    for(i in 1:bins) {
        if(rel.freq[i] != 0) {
            lines.default(c(0, radius[i] * cos(sector[i])), c(0, radius[i] * sin(sector[i])), col=col[iseries], ...)
            lines.default(c(0, radius[i] * cos(sector[i] + (2 * pi)/bins)), c(0, radius[i] * sin(sector[i] + (2 * pi)/bins)), col=col[iseries], ...)
            lines.default(c(radius[i] * cos(sector[i]), radius[i] * cos(sector[i] + (2 * pi)/bins)), c(radius[i] * sin(sector[i]), radius[i] * sin(sector[i] + (2 * pi)/bins)), col=col[iseries], ...)
        }
    }
   }
return(invisible(list(zero=zero, rotation=rotation, next.points=0)))    
}
###############################################################
#       rstable function                                      #
#       Date: January, 22, 2002                               #
#       Version: 0.1                                          #
#                                                             #
###############################################################
#                                                             #
#   This  R code is based on C functions gsl_ran_levy and     #
#     gsl_ran_levy_skew from GNU Scientifi Library            #
#      copyrighted under GNU general license by               #
#     James Theiler, Brian Gough and Keith Briggs.            #
#                                                             #
###############################################################     
#   Here the original comments in the code:
#
#   The stable Levy probability distributions have the form
#
#   p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
#
#   with 0 < alpha <= 2. 
#
#   For alpha = 1, we get the Cauchy distribution
#   For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.
#
#   Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
#   Simulation". The original reference given there is,
#
#   J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
#   simulating stable random variates". Journal of the American
#   Statistical Association, JASA 71 340-344 (1976).
#
#   The following routine for the skew-symmetric case was provided by
#   Keith Briggs.
#
#   The stable Levy probability distributions have the form
#
#   2*pi* p(x) dx
#
#     = \int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for alpha!=1
#     = \int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|)))   for alpha==1
#
#   with 0<alpha<=2, -1<=beta<=1, sigma>0.
#
#   For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
#
#   For alpha = 1, beta=0, we get the Lorentz distribution
#   For alpha = 2, beta=0, we get the Gaussian distribution
#
#   See A. Weron and R. Weron: Computer simulation of Lvy alpha-stable 
#   variables and processes, preprint Technical University of Wroclaw.
#   http://www.im.pwr.wroc.pl/~hugo/Publications.html
#
###############################################################################


rstable <- function(n, scale = 1, index = stop("no index arg"), skewness = 0) {

alpha <- index
beta <- skewness

if (alpha > 2 | alpha <= 0) {stop("rstable is not define for index outside the interval 0 < index <= 2\n")}

if (beta > 1 | beta < -1) {stop("rstable is not define for skewness outside the interval -1 <= skewness <= 1\n")}

if (beta==0) {

## cauchy case
  if (alpha == 1) { 
      return(scale*rcauchy(n, location = 0, scale = 1))
  }

## gaussian case
  if (alpha == 2) { 
      return(rnorm(n, mean = 0, sd = sqrt(2)*scale)) 
  }

## general case

  rngstab <- vector(length=0)
  for (i in 1:n) {
       u <- 0 
       while (u == 0 | u == 1) {
              u <-  pi * (runif(1, min=0, max=1) - 0.5)
       }

       v <- 0
       while (v == 0) {
              v <- rexp(1,rate=1)   
       }

       t <-  sin (alpha * u) / (cos (u)^(1 / alpha))
       s <- (cos ((1 - alpha) * u) / v)^((1 - alpha) / alpha)
  rngstab <- c(rngstab, t*s)
  }

  return (scale * rngstab);
} else {

   rngstab <- vector(length=0)
   for (i in 1:n) {

       u <- 0 
       while (u == 0 | u == 1) {
              u <-  pi * (runif(1, min=0, max=1) - 0.5)
       }

       v <- 0
       while (v == 0) {
              v <- rexp(1,rate=1)   
       }

       if (alpha == 1) {
           X <-  (((pi/2) + beta * u) * tan (u) - beta * log ((pi/2) * v * cos (u) / ((pi/2) + beta * u))) / (pi/2)
           rngstab <- c(rngstab, (scale * (X + beta * log (scale) / (pi/2))))
       } else {
           t <- beta * tan ((pi/2) * alpha)
           B <- atan (t) / alpha
           S <-  (1 + t * t)^(1/(2 * alpha))

           X <-  S * sin (alpha * (u + B)) / (cos (u)^(1 / alpha)) * (cos (u - alpha * (u + B)) / v)^((1 - alpha) / alpha)
          rngstab <- c(rngstab, (scale * X))
      }
  }
return(rngstab)
}

}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   summary.circular                                        #
#   Authors: Claudio Agostinelli, David Andel               #
#   Email: claudio@unive.it, andel@ifi.unizh.ch             #
#   Date: August, 01, 2003                                  #
#   Copyright (C) 2003 Claudio Agostinelli, David Andel     #
#                                                           #
#   Version 0.3-1                                           #
#############################################################

summary.circular <- function(object, ...) {
  if (is.matrix(object)) {
    return(summary.matrix(object, ...))
  }
  else {
    nas <- is.na(object)
    object <- object[!nas]
    n <- length(object)
    result <- c(n, mean.circular(object), rho.circular(object))
    names(result) <- c("n", "Mean", "Rho")
    if(any(nas))
      c(result, "NA's" = sum(nas))
    else result
  }
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rtriangular function                                    #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

rtriangular <- function(n, rho, units=c("radians", "degrees"), ...) {
        if (rho < 0 | rho > 4/pi^2) stop("'rho' must be between 0 and 4/pi^2")
        units <- match.arg(units)
    u <- matrix(c(runif(n)), ncol = 1)
    get.theta <- function(u, rho)
    {
        if(u < 0.5) {
            a <- pi * rho
            b <-  - (4 + pi^2 * rho)
            c <- 8 * pi * u
            theta1 <- ( - b + sqrt(b^2 - 4 * a * c))/(2 * a)
            theta2 <- ( - b - sqrt(b^2 - 4 * a * c))/(2 * a)
            min(theta1, theta2)
        }
        else {
            a <- pi * rho
            b <- 4 - 3 * pi^2 * rho
            c <- (2 * pi^3 * rho) - (8 * pi * u)
            theta1 <- ( - b + sqrt(b^2 - 4 * a * c))/(2 * a)
            theta2 <- ( - b - sqrt(b^2 - 4 * a * c))/(2 * a)
            max(theta1, theta2)
        }
    }
    theta <- apply(u, 1, get.theta, rho)
    theta[theta > pi] <- theta[theta > pi] - 2 * pi
        if (units=="degrees") theta <- theta*180/pi
    return(circular(theta, units=units, ...))
}

#############################################################
#                                                           #
#   dtriangular function                                    #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

dtriangular <- function(x, rho) {
    if (rho < 0 | rho > 4/pi^2) stop("'rho' must be between 0 and 4/pi^2")
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    attr(x, "circularp") <-  NULL
    d <- (4 - pi^2 * rho + 2 * pi * rho * abs(pi - x))/(8 * pi)
    d <- unclass(d)
    return(d)

}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   trigonometric.moment function                           #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

trigonometric.moment <- function(x, p = 1, center = FALSE) {
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    attr(x, "circularp") <-  NULL

    sinr <- sum(sin(x))
    cosr <- sum(cos(x))
    circmean <- atan(sinr, cosr)
    sin.p <- sum(sin(p * (x - circmean * center)))/n
    cos.p <- sum(cos(p * (x - circmean * center)))/n
    mu.p <- atan(sin.p, cos.p)
    rho.p <- sqrt(sin.p^2 + cos.p^2)
    if (units=="degrees") mu.p <- mu.p/pi*180
    attr(mu.p, "circularp") <- xcircularp
    attr(mu.p, "class") <- "circular"
    return(list(mu=mu.p, rho=rho.p, cos=cos.p, sin=sin.p, p=p, n=n, call=match.call()))
}

var <- function(x, ...) UseMethod("var")

var.default <- function(x, y = NULL, na.rm = FALSE, use, ...) base::var(x=x, y=y, na.rm=na.rm, use=use)

#var.matrix <- function(x, ...) {
#    apply(x, 2, var, ...)
#}

var.data.frame <- function(x, ...) {
    sapply(x, var, ...)
}


###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   var.circular function                                   #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: November, 19, 2003                                #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

var.circular <- function (x, na.rm = FALSE, only.var=TRUE, ...)  {
    x <- as.circular(x)
    x <- conversion.circular(x, units="radians")

    if (is.matrix(x)) {
        apply(x, 2, var.circular, na.rm=na.rm, only.var=only.var)
    } else {

       if (na.rm) 
           x <- x[!is.na(x)]
       if (any(is.na(x))) return(NA)

       n <- length(x)
       c <- sum(cos(x))
       s <- sum(sin(x))
       r <- sqrt(c^2 + s^2)
       rbar <- r/n
       circvar <- 1 - rbar
       if (only.var) {
           return(circvar)
       } else {
           return(c(n,r,rbar,circvar))
       }
   }
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rvonmises function                                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 21, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

rvonmises <- function(n, mu, kappa, units=c("radians", "degrees"), ...) {
    units <- match.arg(units)
    if (units=="degrees") {
        mu <- mu/180*pi
    }

    vm <- 1:n
    a <- 1 + (1 + 4 * (kappa^2))^0.5
    b <- (a - (2 * a)^0.5)/(2 * kappa)
    r <- (1 + b^2)/(2 * b)
    obs <- 1
    while (obs <= n) {
       U1 <- runif(1, 0, 1)
       z <- cos(pi * U1)
       f <- (1 + r * z)/(r + z)
       c <- kappa * (r - f)
       U2 <- runif(1, 0, 1)
       if (c * (2 - c) - U2 > 0) {
           U3 <- runif(1, 0, 1)
           vm[obs] <- sign(U3 - 0.5) * acos(f) + mu
           vm[obs] <- vm[obs] %% (2 * pi)
           obs <- obs + 1
       } else {
           if (log(c/U2) + 1 - c >= 0) {
           U3 <- runif(1, 0, 1)
           vm[obs] <- sign(U3 - 0.5) * acos(f) + mu
           vm[obs] <- vm[obs] %% (2 * pi)
           obs <- obs + 1
           }
       }
    }
    if (units=="degrees") vm <- vm/pi*180
    vm <- circular(vm, units=units, ...)
    return(vm)
}

#############################################################
#                                                           #
#   dvonmises function                                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 21, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

dvonmises <- function (x, mu, kappa) {
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    attr(x, "class") <- attr(x, "circularp") <-  NULL
    
    if (units=="degrees") {
        mu <- mu/180*pi
    }  
  
    return(1/(2 * pi * besselI(x = kappa, nu = 0, expon.scaled = TRUE)) * (exp(cos(x - mu) -1))^kappa)
}

#############################################################
#                                                           #
#   pvonmises function                                      #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-1                                           #
#############################################################

pvonmises <- function(q, mu, kappa, tol = 1e-020) {

     q <- as.circular(q)
     qcircularp <- circularp(q)
     units <- qcircularp$units
     q <- conversion.circular(q, units="radians")
     n <- length(q)
     attr(q, "class") <- attr(q, "circularp") <-  NULL
    
     if (units=="degrees") {
         mu <- mu/180*pi
     }  
     attr(mu, "class") <- attr(mu, "circularp") <-  NULL
 
     q <- q %% (2 * pi)
     mu <- mu %% (2 * pi)
     pvm.mu0 <- function(q, kappa, tol) {
        flag <- TRUE
        p <- 1
        sum <- 0
        while (flag) {
           term <- (besselI(x=kappa, nu=p, expon.scaled = FALSE) * sin(p * q))/p
           sum <- sum + term
           p <- p + 1
           if (abs(term) < tol)
           flag <- FALSE
    }
    return(q/(2 * pi) + sum/(pi * besselI(x=kappa, nu=0, expon.scaled = FALSE)))
     }
     result <- rep(NA, n)
     if (mu == 0) {
         for (i in 1:n) {
          result[i] <- pvm.mu0(q[i], kappa, tol)
         }
     } else {
         for (i in 1:n) {
           
         if (q[i] <= mu) {
             upper <- (q[i] - mu) %% (2 * pi)
             if (upper == 0)
             upper <- 2 * pi
             lower <- ( - mu) %% (2 * pi)
             result[i] <- pvm.mu0(upper, kappa, tol) - pvm.mu0(lower, kappa, tol)
             } else {
             upper <- q[i] - mu
             lower <- mu %% (2 * pi)
             result[i] <- pvm.mu0(upper, kappa, tol) + pvm.mu0(lower, kappa, tol)
            }
         }     
     }
    return(result)
}

#############################################################
#                                                           #
#   dmixedvonmises function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 21, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

dmixedvonmises <- function(x, mu1, mu2, kappa1, kappa2, p) {
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    attr(x, "class") <- attr(x, "circularp") <-  NULL
    
    if (units=="degrees") {
        mu1 <- mu1/180*pi
        mu2 <- mu2/180*pi
    }  
  
    return(p/(2 * pi * besselI(x=kappa1, nu=0, expon.scaled = TRUE)) * (exp(cos(x - mu1) - 1))^kappa1 + (1 - p)/(2 * pi * besselI(x=kappa2, nu=0, expon.scaled = TRUE)) * (exp(cos(x - mu2) - 1))^kappa2)
}

#############################################################
#                                                           #
#   rmixedvonmises function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 21, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

rmixedvonmises <- function(n, mu1, mu2, kappa1, kappa2, p, units=c("radians", "degrees"), ...) {
    units <- match.arg(units)
    result <- rep(NA, n)
    test <- runif(n)
    n1 <- sum(test < p)
    n2 <- n - n1
    res1 <- rvonmises(n1, mu1, kappa1, units=units, ...)
    res2 <- rvonmises(n2, mu2, kappa2, units=units, ...)
    result[test < p] <- res1
    result[test >= p] <- res2
    return(result)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   watson.test function                                    #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 27, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

watson.test <- function(x, alpha = 0, dist = c("uniform", "vonmises")) {
    dist <- match.arg(dist)
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    if (!any(c(0, 0.01, 0.025, 0.05, 0.1)==alpha)) stop("'alpha' must be one of the following values: 0, 0.01, 0.025, 0.05, 0.10")
    n <- length(x)

    if (dist == "uniform") {
    u <- sort(x)/(2 * pi)
    u.bar <- mean.default(u)
    i <- 1:n
    sum.terms <- (u - u.bar - (2 * i - 1)/(2 * n) + 0.5)^2
    u2 <- sum(sum.terms) + 1/(12 * n)
    u2 <- (u2 - 0.1/n + 0.1/(n^2)) * (1 + 0.8/n)
        result <- list(statistic=u2, alpha=alpha, n=n, dist=dist, row=NA)
    } else {
        res <- mle.vonmises(x, bias=FALSE)
    mu.hat <- res$mu
    kappa.hat <- res$kappa
    x <- (x - mu.hat) %% (2 * pi)
    x <- matrix(x, ncol = 1)
    z <- apply(x, 1, pvonmises, mu=0, kappa=kappa.hat)
    z <- sort(z)
    z.bar <- mean.default(z)
    i <- 1:n
    sum.terms <- (z - (2 * i - 1)/(2 * n))^2
    Value <- sum(sum.terms) - n * (z.bar - 0.5)^2 + 1/(12 * n)                
    if (kappa.hat < 0.25)
        row <- 1
    else if (kappa.hat < 0.75)
            row <- 2
        else if (kappa.hat < 1.25)
            row <- 3
        else if (kappa.hat < 1.75)
            row <- 4
        else if (kappa.hat < 3)
            row <- 5
        else if (kappa.hat < 5)
            row <- 6
        else row <- 7   
        result <- list(statistic=Value, alpha=alpha, n=n, dist=dist, row=row)
    }
    class(result) <-"watson.test"
    return(result)
}

#############################################################
#                                                           #
#   print.watson.test function                              #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.watson.test <- function(x, digits=4, ...) {
    dist <- x$dist
    n <- x$n
    alpha <- x$alpha

    if (dist == "uniform") {
        u2 <- x$statistic
    cat("\n", "      Watson's Test for Circular Uniformity", "\n", "\n")
    crits <- c(99, 0.267, 0.221, 0.187, 0.152)
    if (n < 8) {
        warning("Total Sample Size < 8:  Results may not be valid", "\n", "\n")
    }
        cat("Test Statistic:", round(u2, digits=digits), "\n")
    if (alpha == 0) {
        if (u2 > 0.267)
        cat("P-value < 0.01", "\n", "\n")
        else if (u2 > 0.221)
             cat("0.01 < P-value < 0.025", "\n", "\n")
        else if (u2 > 0.187)
             cat("0.025 < P-value < 0.05", "\n", "\n")
        else if (u2 > 0.152)
             cat("0.05 < P-value < 0.10", "\n", "\n")
        else cat("P-value > 0.10", "\n", "\n")
    } else {
        index <- (1:5)[alpha == c(0, 0.01, 0.025, 0.05, 0.1)]
        Critical <- crits[index]
        if (u2 > Critical)
        Reject <- "Reject Null Hypothesis"
        else Reject <- "Do Not Reject Null Hypothesis"
        cat("Level", alpha, "Critical Value:", round(Critical, digits=digits), "\n")
        cat(Reject, "\n\n")
        }
    
    } else if (dist=="vonmises") {
               Value <- x$statistic
               row <- x$row
           cat("\n", "      Watson's Test for the von Mises Distribution \n\n")
               u2.crits <- cbind(c(0, 0.5, 1, 1.5, 2, 4, 100), c(0.052, 0.056, 0.066, 0.077, 0.084, 0.093, 0.096), c(0.061, 0.066, 0.079, 0.092, 0.101, 0.113, 0.117), c(0.081, 0.09, 0.11, 0.128, 0.142, 0.158, 0.164))
 
           if (alpha != 0) {
           if (alpha == 0.1)
               col <- 2
           else if (alpha == 0.05)
               col <- 3
           else if (alpha == 0.01)
            col <- 4
           Critical <- u2.crits[row, col]
           if (Value > Critical)
               Reject <- "Reject Null Hypothesis"
           else Reject <- "Do Not Reject Null Hypothesis"
           cat("Test Statistic:", round(Value, digits=digits), "\n")
           cat("Level", alpha, "Critical Value:", round(Critical, digits=digits), "\n")
           cat(Reject, "\n\n")
           } else {
           cat("Test Statistic:", round(Value, digits=digits), "\n")
           if (Value < u2.crits[row, 2])
               cat("P-value > 0.10", "\n", "\n")
           else if ((Value >= u2.crits[row, 2]) && (Value < u2.crits[row, 3]))
                cat("0.05 < P-value > 0.10", "\n", "\n")
           else if ((Value >= u2.crits[row, 3]) && (Value < u2.crits[row, 4]))
                cat("0.01 < P-value > 0.05", "\n", "\n")
           else cat("P-value < 0.01", "\n", "\n")
           }

           }
    invisible(x)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   watson.two.test function                                #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: July, 27, 2003                                    #
#   Version: 0.1                                            #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

watson.two.test <- function(x, y, alpha = 0) {
    x <- as.circular(x)
    x <- conversion.circular(x, units="radians")
    attr(x, "circularp") <- attr(x, "class") <- NULL
    y <- as.circular(y)
    y <- conversion.circular(y, units="radians")
    attr(y, "circularp") <- attr(y, "class") <- NULL

    n1 <- length(x)
    n2 <- length(y)
    n <- n1 + n2
    x <- cbind(sort(x %% (2 * pi)), rep(1, n1))
    y <- cbind(sort(y %% (2 * pi)), rep(2, n2))
    xx <- rbind(x, y)
    rank <- order(xx[, 1])
    xx <- cbind(xx[rank,  ], 1:n)
    a <- 1:n
    b <- 1:n
    for (i in 1:n) {
     a[i] <- sum(xx[1:i, 2] == 1)
     b[i] <- sum(xx[1:i, 2] == 2)
    }
    d <- b/n2 - a/n1
    dbar <- mean.default(d)
    u2 <- (n1 * n2)/n^2 * sum((d - dbar)^2)
    result <- list(statistic=u2, alpha=alpha, nx=n1, ny=n2)
    class(result) <- "watson.two.test"
    return(result)
}

#############################################################
#                                                           #
#   print.watson.two.test function                          #
#   Author: Claudio Agostinelli                             #
#   E-mail: claudio@unive.it                                #
#   Date: November, 19, 2003                                #
#   Version: 0.1-1                                          #
#                                                           #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#############################################################

print.watson.two.test <- function(x, digits=4, ...) {
    u2 <- x$statistic
    n1 <- x$nx
    n2 <- x$ny
    alpha <- x$alpha
    n <- n1 + n2
    cat("\n      Watson's Two-Sample Test of Homogeneity \n\n")
    if (n < 18)
    warning("Total Sample Size < 18:  Consult tabulated critical values \n\n")
    crits <- c(99, 0.385, 0.268, 0.187, 0.152)
    cat("Test Statistic:", round(u2, digits=digits), "\n")
    
    if (alpha == 0) {
    if (u2 > 0.385)
        cat("P-value < 0.001", "\n", "\n")
    else if (u2 > 0.268)
         cat("0.001 < P-value < 0.01", "\n", "\n")
    else if (u2 > 0.187)
         cat("0.01 < P-value < 0.05", "\n", "\n")
    else if (u2 > 0.152)
         cat("0.05 < P-value < 0.10", "\n", "\n")
    else cat("P-value > 0.10", "\n", "\n")
    } else {
    index <- (1:5)[alpha == c(0, 0.001, 0.01, 0.05, 0.1)]
    Critical <- crits[index]
    if (u2 > Critical)
        Reject <- "Reject Null Hypothesis"
    else Reject <- "Do Not Reject Null Hypothesis"
    cat("Level", alpha, "Critical Value:", round(Critical, digits=digits), "\n")
    cat(Reject, "\n\n")
    }
invisible(x)
}

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rwrappedcauchy function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

rwrappedcauchy <- function(n, mu = 0, rho = exp(-1), units=c("radians", "degrees"), ...) {
    units <- match.arg(units)
    if (units=="degrees") {
        mu <- mu/180*pi
    }
    
    if (rho == 0)
    result <- runif(n, 0, 2 * pi)
    else if (rho == 1)
         result <- rep(mu, n)
    else {
       scale <-  - log(rho)
       result <- rcauchy(n, mu, scale) %% (2 * pi)
    }
    if (units=="degrees") result <- result/pi*180
    result <- circular(result, units=units, ...)
    return(result)
}

#############################################################
#                                                           #
#   dwrappedcauchy function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: July, 24, 2003                                    #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1                                             #
#############################################################

dwrappedcauchy <- function(x, mu=0, rho=exp(-1)) {
    if (rho < 0 | rho > 1)
        stop("rho must be between 0 and 1")
   
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    attr(x, "circularp") <-  NULL
    
    if (units=="degrees") {
        mu <- mu/180*pi
    }  
    d <- (1 - rho^2)/((2 * pi) * (1 + rho^2 - 2 * rho * cos(x - mu)))
    d <- unclass(d)
    return(d)

  }

###############################################################
#                                                             #
#       Original Splus: Ulric Lund                            #
#       E-mail: ulund@calpoly.edu                             #
#                                                             #
###############################################################

#############################################################
#                                                           #
#   rwrappednormal function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: December, 17, 2003                                #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-1                                           #
#############################################################

rwrappednormal <- function(n, mu=0, rho, sd=1, units=c("radians", "degrees"), ...) {
    units <- match.arg(units)
    if (units=="degrees") {
        mu <- mu/180*pi
    }
    if (missing(rho)) {
        rho <- exp(-sd^2/2)
    }
    if (rho < 0 | rho > 1)
        stop("rho must be between 0 and 1")        
    if (rho == 0)
        result <- runif(n, 0, 2 * pi)
    else if (rho == 1)
        result <- rep(mu, n)
    else {
        sd <- sqrt(-2 * log(rho))
        result <- rnorm(n, mu, sd) %% (2 * pi)
    }
    if (units=="degrees") result <- result/pi*180
    result <- circular(result, units=units, ...)
    return(result)
}

#############################################################
#                                                           #
#   dwrappednormal function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: December, 17, 2003                                #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.1-2                                           #
#############################################################

dwrappednormal <- function(x, mu=0, rho, sd=1, K, min.k=10) {
  
    x <- as.circular(x)
    xcircularp <- circularp(x)
    units <- xcircularp$units
    x <- conversion.circular(x, units="radians")
    n <- length(x)
    im <- length(mu)
    attr(x, "class") <- attr(x, "circularp") <-  NULL
    
    if (units=="degrees") {
        mu <- mu/180*pi
    }  
    if (missing(rho)) {
        rho <- exp(-sd^2/2)
    }
    if (rho < 0 | rho > 1)
        stop("rho must be between 0 and 1")
    var <- -2 * log(rho)
    sd <- sqrt(var)
    
    if (missing(K)) {
        range <- abs(mu-x)
        K <- (range+6*sqrt(var))%/%(2*pi)+1
        K <- max(min.k, K)
    }

    z <- .Fortran("dwrpnorm",
                    as.double(x),
                    as.double(mu),
                    as.double(sd),
                    as.integer(n),
                    as.integer(im),
                    as.integer(K),
                    d=mat.or.vec(im, n),
                    PACKAGE="circular"
    )
    d <- z$d/sqrt(var * 2 * pi)

## we have to be careful to return the right dimension when x is a matrix and mu is a vector 
#    if (!is.null(dim(x))) {
#         d <- matrix(d, dim(x))
#    }
    return(d)
}
#############################################################
#                                                           #
#   rwrappedstable function                                 #
#   Author: Claudio Agostinelli                             #
#   Email: claudio@unive.it                                 #
#   Date: September, 22, 2003                               #
#   Copyright (C) 2003 Claudio Agostinelli                  #
#                                                           #
#   Version 0.2                                             #
#############################################################

rwrappedstable <- function(n,  scale=1, index, skewness, units=c("radians", "degrees"), ...) {
    units <- match.arg(units)
#    if (units=="degrees") {
#        scale <- scale/180*pi
#    }
    result <- rstable(n=n, scale=scale, index=index, skewness=skewness) %% (2 * pi)
    if (units=="degrees") result <- result/pi*180
    result <- circular(result, units=units, ...)
    return(result)
}

.First.lib <- function(lib, pkg) {
  library.dynam("circular", pkg, lib)
  cat("This is version 0.1 of circular package under tests \n")
  cat("Please report any bugs or comments to <Claudio Agostinelli> claudio@unive.it \n")
  cat("The package redefine how function 'density' and 'var' works \n")
  cat("In particular, for 'var' function (try 'methods(var)') \n notice that 'var.default' is an alias for the original 'var' function \n and that a method for data.frame is available.\n")
}
