###############################################################
## Modified March 5, 2002
# Based on rose.diag from CircStats Library (R-port)


vera2.diag <- function(x,s,breaks=c(9,6,3,1.5),color=c(2,3,4,5),bins, main = "", prop = 1, cex = 1, pch = 16, dotsep = 40, shrink = 1) {

        x <- x * (pi/180) # convert to radian
	#draw outer circle
	eqscplot(cos(seq(0, 2 * pi, length = 1000)), sin(seq(0, 2 * pi, length = 1000)),axes = FALSE, xlab = "", ylab = "",
	main = main, type = "l", xlim = shrink * c(-1, 1), ylim = shrink * c(-1, 1),lty=1)
#draw ticks and lables
	lines(c(0, 0), c(0.95, 1))
	text(0.005, 0.9, "N", cex = 1.5)
	lines(c(0, 0), c(-0.95, -1))
	text(0.005, -0.9, "S", cex = 1.5)
	lines(c(-1, -0.95), c(0, 0))
	text(-0.9, 0, "W", cex = 1.5)
	lines(c(0.95, 1), c(0, 0))
	text(0.9, 0, "E", cex = 1.5)
y<-breaks 
#assign to xlow only the values of x which are less then y
j<-1
i<-1
g<-1
xlow<-c(1:length(s))
for(j in 1:length(y)+1){
cy<-color[j]
if(j==1){
  xlow<-x}
else{

for(g in 1:length(s))
if (!is.na(s[g])& !is.na(y[j]) & s[g]< y[j]){
xlow[g]<- (x[g])
}else{
xlow[g]<-NA}

}

        n <- length(x) #number of data points
	freq <- c(1:bins) # 1,2,3,4....bins
	arc <- (2 * pi)/bins 

        for(i in 1:bins) {
            freq[i] <- sum(xlow == i*arc & xlow>(i-1)*arc)
            
        } # finds out the number of  values in each bin.
        print(freq)
        rel.freq <- freq/n # percentage of occurences in each segment
	radius <- sqrt(rel.freq) * prop  # proportional distance from center
	sector <- seq((2 * pi - (2 * pi)/bins)+pi/2+arc/2, 0+pi/2+arc/2, length = bins) #
	mids <- seq(arc/2, 2 * pi - pi/bins, length = bins) # middle of sectors (angle in radians)
	degmids <- mids / (pi/180)
	index <- cex/dotsep  #calculates radial spacing between dots (ie. 0.025)
	for(i in 1:bins) {

		if(rel.freq[i] != 0 & !is.na(rel.freq[i])) {
                startsec<-sector[i]+arc/4
		endsec<-sector[i]-arc/4 + ((2 * pi)/bins)
		x1<-0
		y1<-0

		x2<-radius[i] * cos(startsec)
		y2<-radius[i] * sin(startsec)

		x3<-radius[i] * cos(endsec)	
		y3<-radius[i] * sin(endsec)

		polygon(c(x1,x2,x3),c(y1,y2,y3),col=(cy))
			
              }
	}
}
      }