.packageName <- "plotrix"
# axis.break places a break marker at the position "breakpos" 
# in user coordinates on the axis nominated - see axis().

axis.break<-function(axis=1,breakpos=NULL,bgcol="white",breakcol="black",
 style="slash",brw=0.02) {
 
 # get the coordinates of the outside of the plot
 figxy<-par("usr")
 # flag if either axis is logarithmic
 xaxl<-par("xlog")
 yaxl<-par("ylog")
 # calculate the x and y offsets for the break
 xw<-(figxy[2]-figxy[1])*brw
 yw<-(figxy[4]-figxy[3])*brw
 # if no break position was given, put it just off the plot origin
 if(is.null(breakpos))
  breakpos<-ifelse(axis%%2,figxy[1]+xw*2,figxy[3]+yw*2)
 if(xaxl && (axis == 1 || axis == 3)) breakpos<-log10(breakpos)
 if(yaxl && (axis == 2 || axis == 4)) breakpos<-log10(breakpos)
 # set up the "blank" rectangle (left, bottom, right, top)
 switch(axis,
  br<-c(breakpos-xw/2,figxy[3]-yw/2,breakpos+xw/2,figxy[3]+yw/2),
  br<-c(figxy[1]-xw/2,breakpos-yw/2,figxy[1]+xw/2,breakpos+yw/2),
  br<-c(breakpos-xw/2,figxy[4]-yw/2,breakpos+xw/2,figxy[4]+yw/2),
  br<-c(figxy[2]-xw/2,breakpos-yw/2,figxy[2]+xw/2,breakpos+yw/2),
  stop("Improper axis specification."))
 # get the current setting of xpd
 old.xpd<-par("xpd")
 # don't cut the break off at the edge of the plot
 par(xpd=TRUE)
 # correct for logarithmic axes
 if(xaxl) br[c(1,3)]<-10^br[c(1,3)]
 if(yaxl) br[c(2,4)]<-10^br[c(2,4)]
 # draw the "blank" rectangle
 rect(br[1],br[2],br[3],br[4],col=bgcol,border=bgcol)
 if(style == "slash") {
  # calculate the slash ends
  if(axis == 1 || axis == 3) {
   xbegin<-c(breakpos-xw,breakpos)
   xend<-c(breakpos,breakpos+xw)
   ybegin<-c(br[2],br[2])
   yend<-c(br[4],br[4])
   if(xaxl) {
    xbegin<-10^xbegin
    xend<-10^xend
   }
  }
  else {
   xbegin<-c(br[1],br[1])
   xend<-c(br[3],br[3])
   ybegin<-c(breakpos-yw,breakpos)
   yend<-c(breakpos,breakpos+yw)
   if(yaxl) {
    ybegin<-10^ybegin
    yend<-10^yend
   }
  }
 }
 else {
  # calculate the zigzag ends
  if(axis == 1 || axis == 3) {
   xbegin<-c(breakpos-xw/2,breakpos-xw/4,breakpos+xw/4)
   xend<-c(breakpos-xw/4,breakpos+xw/4,breakpos+xw/2)
   ybegin<-c(ifelse(yaxl,10^figxy[3+(axis==3)],figxy[3+(axis==3)]),br[4],br[2])
   yend<-c(br[4],br[2],ifelse(yaxl,10^figxy[3+(axis==3)],figxy[3+(axis==3)]))
   if(xaxl) {
    xbegin<-10^xbegin
    xend<-10^xend
   }
  }
  else {
   xbegin<-c(ifelse(xaxl,10^figxy[1+(axis==4)],figxy[1+(axis==4)]),br[1],br[3])
   xend<-c(br[1],br[3],ifelse(xaxl,10^figxy[1+(axis==4)],figxy[1+(axis==4)]))
   ybegin<-c(breakpos-yw/2,breakpos-yw/4,breakpos+yw/4)
   yend<-c(breakpos-yw/4,breakpos+yw/4,breakpos+yw/2)
   if(yaxl) {
    ybegin<-10^ybegin
    yend<-10^yend
   }
  }
 }
 # draw the segments
 segments(xbegin,ybegin,xend,yend,col=breakcol,lty=1)
 # restore xpd
 par(xpd=old.xpd)
}
boxed.labels<-function(x,y,labels,col="white",border=TRUE,xpad=0.6,ypad=0.6,...) {
 widths<-strwidth(labels)
 heights<-strheight(labels)
 rect(x-widths*xpad,y-heights*ypad,x+widths*xpad,y+heights*ypad,
  col=col,border=border)
 text(x,y,labels,...)
}
color.scale<-function(x,redrange,greenrange,bluerange,scale.up=FALSE){
 ncolors<-length(x)
 if(length(redrange) > 1) {
  reds<-rescale(x,redrange)
  if(min(reds) < 0 || max(reds) > 255) reds<-rescale(reds,c(0,255))
 }
 else reds<-rep(redrange,ncolors)
 if(length(greenrange) > 1) {
  greens<-rescale(x,greenrange)
  if(min(greens) < 0 || max(greens) > 255) greens<-rescale(greens,c(0,255))
 }
 else greens<-rep(greenrange,ncolors)
 if(length(bluerange) > 1) {
  blues<-rescale(x,bluerange)
  if(min(blues) < 0 || max(blues) > 255) blues<-rescale(blues,c(0,255))
 }
 else blues<-rep(bluerange,ncolors)
 colormatrix<-cbind(reds,greens,blues)
 colvec<-apply(colormatrix,1,rgb.to.hex,scale.up)
 return(colvec)
}
# display a pie chart at an arbitrary location on an existing plot

floating.pie<-function (xpos,ypos,x,edges=200,radius=1,col=NULL,...) {
 if (!is.numeric(x) || any(is.na(x) | x<=0))
  stop("floating.pie: `sectors' values must be positive.")
 x <- c(0, cumsum(x)/sum(x))
 dx <- diff(x)
 nx <- length(dx)
 if (is.null(col)) col<-rainbow(nx)
 else if(length(col) < nx) col<-rep(col,nx)
 # get the center values in radians
 bc<-2*pi*(x[1:nx]+dx/2)
 for (i in 1:nx) {
  n <- max(2, floor(edges * dx[i]))
  t2p <- 2 * pi * seq(x[i], x[i + 1], length = n)
  xc <- c(cos(t2p)+xpos,xpos) * radius
  yc <- c(sin(t2p)+ypos,ypos) * radius
  polygon(xc, yc, col = col[i],...)
  t2p <- 2 * pi * mean(x[i + 0:1])
  xc <- cos(t2p) * radius
  yc <- sin(t2p) * radius
  lines(c(1, 1.05) * xc, c(1, 1.05) * yc)
 }
 return(bc)
}

# place text labels at the specified distance from x,y on the radial lines
# specified by angles.

pie.labels<-function(x,y,angles,labels,radius=1,col="white",border=TRUE,...) {
 if(nargs()<4)
  stop("Usage: pie.labels(x,y,angles,labels,radius=1,col=\"white\",border=TRUE,...)")
 xc<-cos(angles)*radius+x
 yc<-sin(angles)*radius+y
 boxed.labels(xc,yc,labels,col=col,border=border,...)
}
get.gantt.info<-function(format="%Y/%m/%d") {
 cat("Enter the label, start and finish time for each task.\n")
 cat("Default format for time is year/month/day e.g. 2005/2/22\n")
 cat("Enter a blank label to end.\n")
 nextlabel<-"dummy"
 tasklabels<-NA
 taskstarts<-NA
 taskends<-NA
 priorities<-NA
 while(nchar(nextlabel)) {
  nextlabel<-readline("Task label - ")
  if(nchar(nextlabel)) {
   if(is.na(tasklabels[1])) tasklabels<-nextlabel
   else tasklabels<-c(tasklabels,nextlabel)
   nextstart<-as.POSIXct(strptime(readline("Task begins - "),format=format))
   if(is.na(taskstarts[1])) taskstarts<-nextstart
   else taskstarts<-c(taskstarts,nextstart)
   nextend<-nextstart-1
   while(nextend < nextstart) {
    nextend<-as.POSIXct(strptime(readline("Task ends - "),format=format))
    if(nextend < nextstart) cat("Task cannot end before it starts!\n")
    else {
     if(is.na(taskends[1])) taskends<-nextend
     else taskends<-c(taskends,nextend)
    }
   }
   nextpriority<-0
   while(nextpriority < 1 || nextpriority > 10)
    nextpriority<-as.numeric(readline("Task priority (1-10) - "))
   if(is.na(priorities[1])) priorities<-nextpriority
   else priorities<-c(priorities,nextpriority)
  }
 }
 return(list(labels=tasklabels,starts=taskstarts,ends=taskends,
  priorities=priorities))
}

axis.POSIXct.tickpos<-function (side,x,format,...) {
 x <- as.POSIXct(x)
 range <- par("usr")[if(side%%2) 1:2 else 3:4]
 d <- range[2] - range[1]
 z <- c(range, x[is.finite(x)])
 if (d < 1.1 * 60) {
  sc <- 1
  if(missing(format))
  format <- "%S"
 }
 else if (d < 1.1 * 60 * 60) {
  sc <- 60
  if(missing(format)) format <- "%M:%S"
 }
 else if (d < 1.1 * 60 * 60 * 24) {
  sc <- 60 * 24
  if(missing(format)) format <- "%H:%M"
 }
 else if (d < 2 * 60 * 60 * 24) {
  sc <- 60 * 24
  if(missing(format)) format <- "%a %H:%M"
 }
 else if (d < 7 * 60 * 60 * 24) {
  sc <- 60 * 60 * 24
  if(missing(format)) format <- "%a"
 }
 else {
  sc <- 60 * 60 * 24
 }
 if (d < 60 * 60 * 24 * 50) {
  zz <- pretty(z/sc)
  z <- zz * sc
  class(z) <- c("POSIXt", "POSIXct")
  if(missing(format)) format <- "%b %d"
 }
 else if (d < 1.1 * 60 * 60 * 24 * 365) {
  class(z) <- c("POSIXt", "POSIXct")
  zz <- as.POSIXlt(z)
  zz$mday <- 1
  zz$isdst <- zz$hour <- zz$min <- zz$sec <- 0
  zz$mon <- pretty(zz$mon)
  m <- length(zz$mon)
  m <- rep.int(zz$year[1], m)
  zz$year <- c(m, m + 1)
  z <- as.POSIXct(zz)
  if (missing(format)) format <- "%b"
 }
 else {
  class(z) <- c("POSIXt", "POSIXct")
  zz <- as.POSIXlt(z)
  zz$mday <- 1
  zz$isdst <- zz$mon <- zz$hour <- zz$min <- zz$sec <- 0
  zz$year <- pretty(zz$year)
  z <- as.POSIXct(zz)
  if(missing(format)) format <- "%Y"
 }
 z <- z[z >= range[1] & z <= range[2]]
 return(z)
}

gantt.chart<-function(x=NULL,format="%Y/%m/%d",xlim=NULL,
 taskcolors=NULL,main="",ylab="") {
 oldpar<-par(no.readonly=TRUE)
 if(is.null(x)) x<-get.x(format=format)
 ntasks<-length(x$labels)
 plot.new()
 charheight<-strheight("M",units="inches")
 maxwidth<-max(strwidth(x$labels,units="inches"))*1.5
 if(is.null(xlim)) xlim=range(c(x$starts,x$ends))
 if(is.null(taskcolors))
  taskcolors<-color.gradient(c(255,0),c(0,0),c(0,255),max(x$priorities))
 par(mai=c(0,maxwidth,charheight*5,0.1))
 par(omi=c(0.1,0.1,0.1,0.1))
 plot(x$starts,1:ntasks,xlim=xlim,ylim=c(0.5,ntasks+0.5),
     main="",xlab="",ylab=ylab,axes=FALSE,type="n")
 box()
 if(nchar(main)) mtext(main,3,2)
 axis.POSIXct(3,xlim)
 topdown<-seq(ntasks,1)
 axis(2,at=topdown,labels=x$labels,las=2)
 xrange<-par("usr")[1:2]
 abline(v=axTicks(3),col="darkgray",lty=3)
# abline(v=axis.POSIXct.tickpos(3,xlim),col="darkgray",lty=3)
 half.height <- 0.25
 for(i in 1:ntasks) {
  rect(x$starts[i],topdown[i]-half.height,
   x$ends[i],topdown[i]+half.height,
   col=taskcolors[x$priorities[i]],
   border=FALSE)
 }
 par(oldpar)
 invisible(x)
}
rgb.to.hex<-function(rgb,scale.up=TRUE) {
 if(length(rgb) != 3) stop("rgb must be an rgb triplet")
 if(any(rgb < 0) || any(rgb > 255)) stop("all rgb must be between 0 and 255")
 # if it looks like a 0-1 value AND scale.up is TRUE, get the 0-255 equivalent
 if(all(rgb <= 1) && scale.up) rgb<-rgb*255 
 hexdigit<-c(0:9,letters[1:6])
 return(paste("#",hexdigit[rgb[1]%/%16+1],hexdigit[rgb[1]%%16+1],
  hexdigit[rgb[2]%/%16+1],hexdigit[rgb[2]%%16+1],
  hexdigit[rgb[3]%/%16+1],hexdigit[rgb[3]%%16+1],
  sep="",collapse=""))
}

color.gradient<-function(reds,greens,blues,nslices=50,scale.up=FALSE) {
 maxncol<-max(c(length(reds),length(greens),length(blues)))
 if(maxncol < 2) {
  cat("color.gradient: Must specify at least two values for one color\n")
  return(NULL)
 }
 if(length(reds) < nslices) {
  reds<-approx(reds,n=nslices)$y
  # take care of any values < 0 or > 255
  if(min(reds) < 0 || max(reds) > 255) reds<-rescale(reds,c(0,255))
 }
 else {
  # chop off extra values so they don't mess up cbind()
  if(length(reds) > nslices) reds<-reds[1:nslices]
 }
 if(length(greens) < nslices) {
  greens<-approx(greens,n=nslices)$y
  if(min(greens) < 0 || max(greens) > 255) greens<-rescale(greens,c(0,255))
 }
 else if(length(greens) > nslices) greens<-greens[1:nslices]
 if(length(blues) < nslices) {
  blues<-approx(blues,n=nslices)$y
  if(min(blues) < 0 || max(blues) > 255) blues<-rescale(blues,c(0,255))
 }
 else if(length(blues) > nslices) blues<-blues[1:nslices]
 colormatrix<-cbind(reds,greens,blues)
 colvec<-apply(colormatrix,1,rgb.to.hex,scale.up)
 return(colvec)
}

gradient.rect<-function(xleft,ybottom,xright,ytop,reds,greens,blues,
 nslices=50,gradient="x") {
 # assume that the user will never want black gradients, so scale up
 colvec<-color.gradient(reds,greens,blues,nslices,scale.up=TRUE)
 if(!is.null(colvec)) {
  if(gradient == "x") {
   if(length(xleft) == 1) {
    xinc<-(xright-xleft)/(nslices-1)
    xlefts<-seq(xleft,xright-xinc,length=nslices)
    xrights<-xlefts+xinc
   }
   else {
    xlefts<-xleft
    xrights<-xright
   }
   rect(xlefts,ybottom,xrights,ytop,col=colvec,lty=0)
  }
  else {
   if(length(ybottom) == 1) {
    yinc<-(ytop-ybottom)/(nslices-1)
    ybottoms<-seq(ybottom,ytop-yinc,length=nslices)
    ytops<-ybottoms+yinc
   }
   else {
    ybottoms<-ybottom
    ytops<-ytop
   }
   rect(xleft,ybottoms,xright,ytops,col=colvec,lty=0)
  }
 }
 return(colvec)
}
# linearly transforms a vector of numbers to a new range

rescale<-function(x,newrange) {
 if(missing(x) | missing(newrange)) {
  usage.string<-paste("Usage: rescale(x,newrange)\n",
   "\twhere x is a numeric object and newrange is the new min and max\n",
   sep="",collapse="")
  stop(usage.string)
 }
 if(is.numeric(x) && is.numeric(newrange)) {
  xrange<-range(x)
  if(xrange[1] == xrange[2]) stop("rescale: can't rescale a constant vector!")
  mfac<-(newrange[2]-newrange[1])/(xrange[2]-xrange[1])
  return(newrange[1]+(x-xrange[1])*mfac)
 }
 else {
  warning("Only numeric objects can be rescaled")
  return(x)
 }
}

# plots data as radial lines or a polygon on a 24 hour "clockface" going 
# clockwise. radial.pos and radial.range should be in hours.
# Remember to convert hour/minute values to hour/decimal values.
# example: clock24.plot(rnorm(16)+3,seq(5.5,20.5))

clock24.plot<-function(lengths,clock.pos,rp.type="r",...) {
 npos<-length(lengths)
 # if no positions are given, spread the lines out over the circle 
 if(missing(clock.pos)) clock.pos<-seq(0,24-24/(npos+1),length=npos)
 radial.range<-range(clock.pos)
 radial.range[1]<-radial.range[1]-(radial.range[2]-radial.range[1])/(npos-1)
 newrange<-c(pi*(2.5-radial.range[1]/12),pi*(0.5+(24-radial.range[2])/12))
 # rescale to a range of pi/2 to 2.5*pi
 # starting at "midnight" and going clockwise
 clock.radial.pos<-rescale(c(clock.pos,radial.range),newrange)[1:npos]
 clock.labels<-as.character(seq(100,2400,by=100))
 clock.label.pos<-seq(29*pi/12,pi/2,by=-pi/12)
 radial.plot(lengths,clock.radial.pos,newrange,clock.labels,clock.label.pos,
  rp.type=rp.type,...)
}

# plots data as radial lines or a polygon starting at the right and going
# counterclockwise.
# angles should be given in 0-360 values, use radial.plot for radians
# example: polar.plot(rnorm(20)+3,seq(90,280,by=10))

polar.plot<-function(lengths,polar.pos,labels,label.pos,rp.type="r",...) {
 npos<-length(lengths)
 # if no positions are given, add the average distance between positions so that
 # the first and last line don't overlap
 if(missing(polar.pos)) polar.pos<-seq(0,360-360/(npos+1),length=npos)
 if(missing(labels)) {
  label.pos<-seq(0,340,by=20)
  labels<-as.character(label.pos)
  label.range<-c(0,pi*340/180)
 }
 if(missing(label.pos)) label.pos<-polar.pos
 polar.range<-range(polar.pos)
 newrange<-c(pi*polar.range[1]/180,pi*(2-(360-polar.range[2])/180))
 # rescale to radians
 radial.pos<-rescale(c(polar.pos,polar.range),newrange)[1:npos]
 nlabels<-length(labels)
 label.pos<-rescale(c(label.pos,0,360),c(0,2*pi))[1:nlabels]
 radial.plot(lengths,radial.pos,newrange,labels,label.pos,rp.type=rp.type,...)
}

# plots radial lines of length 'lengths' or a polygon with corresponding
# vertices specified by 'radial.pos' in radians.
# starts at the 'east' position and goes counterclockwise
# label.prop is the proportion of max(lengths) that gives the
# radial position of the labels

radial.plot<-function(lengths,radial.pos,radial.range,labels,label.pos,
 rp.type="r",label.prop=1.1,main="",xlab="",ylab="",...) {
 maxlength<-label.prop*max(lengths)
 # calculate missing range as above
 if(missing(radial.range)) {
  radial.range<-range(radial.pos)
  radial.range[1]<-radial.range[1]-(radial.range[2]-radial.range[1])/
  (length(radial.pos)-1)
 }
 plot(c(-maxlength,maxlength),c(-maxlength,maxlength),type="n",axes=FALSE,
  main=main,xlab=xlab,ylab=ylab,...)
 # get the vector of x positions
 xpos<-cos(radial.pos)*lengths
 # get the vector of y positions
 ypos<-sin(radial.pos)*lengths
 # plot radial lines if rp.type == "r"    
 if(rp.type == "r") segments(0,0,xpos,ypos,...)
 else polygon(xpos,ypos,...)
 if(missing(labels)) {
  if(length(radial.pos) <= 20) {
   labels<-as.character(round(radial.pos,2))
   label.pos<-radial.pos
  }
  else {
   label.pos<-seq(0,1.95*pi,length=19)
   labels<-as.character(round(label.pos,2))
  }
 }
 if(missing(label.pos)) {
  xpos<-cos(radial.pos)*maxlength
  ypos<-sin(radial.pos)*maxlength
 }
 else {
  xpos<-cos(label.pos)*maxlength
  ypos<-sin(label.pos)*maxlength
 }
 text(xpos,ypos,adj=0.5,labels)
}
spread.labels<-function(x,y,labels=NULL,offset,col="white",border=FALSE,...) {
 if(missing(x))
  stop("Usage: spread.labels(x,y,labels=NULL,offset,col=\"white\",...)")
 if(diff(range(x)) < diff(range(y))) {
  sort.index<-sort.list(y)
  x<-x[sort.index]
  y<-y[sort.index]
  ny<-length(y)
  offsets<-rep(c(offset,-offset),ny/2+1)[1:ny]
  newy<-seq(y[1],y[ny],length=ny)
  segments(x+offsets,newy,x,y)
  boxed.labels(x+offsets,newy,labels,col=col,border=border,...)
 }
 else {
  sort.index<-sort.list(x)
  x<-x[sort.index]
  y<-y[sort.index]
  nx<-length(x)
  offsets<-rep(c(offset,-offset),nx/2+1)[1:nx]
  newx<-seq(x[1],x[nx],length=nx)
  segments(newx,y+offsets,x,y)
  boxed.labels(newx,y+offsets,labels,col=col,border=border,...)
 }
}
# staxlab produces staggered axis tick labels
# note that barplot() tends to mess things up by plotting an X axis 
# even when axes=F

staxlab<-function(side=1,at,labels,nlines=2,top.line=0.5,line.spacing=0.8) {
 if(missing(labels)) stop("Usage: staxlab(side=1,at,labels,nlines=2)")
 nlabels<-length(labels)
 if(missing(at)) at<-1:nlabels
 linepos<-rep(top.line,nlines)
 for(i in 2:nlines) linepos[i]<-linepos[i-1]+line.spacing
 linepos<-rep(linepos,ceiling(nlabels/nlines))[1:nlabels]
 axis(side=side,at=at,labels=rep("",nlabels))
 mtext(text=labels,side=side,line=linepos,at=at)
}
# thigmophobe returns the direction (as 1|2|3|4 - see pos= in the text function) 
# _away_ from the nearest point where x and y are vectors of 2D coordinates

thigmophobe<-function(x,y) {
 # get the current upper and lower limits of the plot
 plot.span<-par("usr")
 x.span<-plot.span[2] - plot.span[1]
 y.span<-plot.span[4] - plot.span[3]
 # if either axis is logarithmic, transform the values into logarithms
 if(par("xlog")) x<-log(x)
 if(par("ylog")) y<-log(y)
 # scale the values to the plot span
 # this avoids the numerically larger
 # axis dominating the distance measure
 x<-x/x.span
 y<-y/y.span
 # get the distance matrix as a full matrix
 xy.dist<-as.matrix(dist(cbind(x,y)))
 lenx<-length(x)
 nearest.index<-rep(0,lenx)
 for(index in 1:lenx)
  nearest.index[index]<-as.numeric(names(which.min(xy.dist[-index,index])))
 # get the x and y differences for each point to the nearest point
 xdiff<-x - x[nearest.index]
 ydiff<-y - y[nearest.index]
 # first set the east/west direction
 dir.ew<-ifelse(xdiff > 0,4,2)
 # now do the north/south
 dir.ns<-ifelse(ydiff > 0,3,1)
 dir.away<-ifelse(abs(xdiff)>abs(ydiff),dir.ew,dir.ns)
 # set any congruent points to N/S labels or they'll overprint
 for(i in 1:lenx) {
  if(!xdiff[i] & !ydiff[i])
   dir.away[c(i,nearest.index[i])]<-c(1,3)
 }
 return(dir.away)
}

# thigmophobe.labels positions labels at points so that they
# are most distant from the nearest other point, where the
# points are described as x and y coordinates.

thigmophobe.labels<-function(x,y,labels=1:length(x),...) {
 
 if(!missing(x) && !missing(y)) {
  text.pos<-thigmophobe(x,y)
  text(x,y,labels,pos=text.pos,...)
 }
 else
  cat("Usage: thigmophobe.labels(x,y,labels=1:length(x),...)\n")
}
