# function plots a tree
# written by Liam Revell 2011

plotTree<-function(tree){
	# make binary
	tree<-multi2di(tree)
	# reorder
	tree<-reorder(tree)
	# count nodes and tips
	n<-length(tree$tip); m<-tree$Nnode
	# Y coordinates for nodes
	Y<-matrix(NA,m+n,1)
	# first, assign y coordinates to all the tip nodes
	Y[1:n]<-1:n
	# now reorder tree "pruningwise"
	tree<-reorder(tree,"pruningwise")
	# get Y coordinates of the nodes
	nodes<-unique(tree$edge[,1])
	for(i in 1:m){
		desc<-tree$edge[which(tree$edge[,1]==nodes[i]),2]
		Y[nodes[i]]<-(min(Y[desc])+max(Y[desc]))/2
	}
	# reorder "cladewise"
	tree<-reorder(tree)
	# compute node heights
	root<-length(tree$tip)+1
	node.height<-matrix(NA,nrow(tree$edge),2)
	for(i in 1:nrow(tree$edge)){
		if(tree$edge[i,1]==root){
			node.height[i,1]<-0.0
			node.height[i,2]<-tree$edge.length[i]
		} else {
			node.height[i,1]<-node.height[match(tree$edge[i,1],tree$edge[,2]),2]
			node.height[i,2]<-node.height[i,1]+tree$edge.length[i]
		}
	}
	# open plot
	plot.new(); plot.window(xlim=c(0,max(node.height)*1.05),ylim=c(0,max(Y)))
	for(i in 1:nrow(tree$edge)){
		points(node.height[i,],c(Y[tree$edge[i,2]],Y[tree$edge[i,2]]))
 		lines(node.height[i,],c(Y[tree$edge[i,2]],Y[tree$edge[i,2]]))
	}
	for(i in 1:m) lines(node.height[which(tree$edge[,1]==nodes[i]),1],Y[tree$edge[which(tree$edge[,1]==nodes[i]),2]])
	for(i in 1:n) text(node.height[which(tree$edge[,2]==i),2],Y[i],tree$tip.label[i],pos=4)
}