vcv.phylo() with ancestral nodes

An R-sig-phylo user just posted a request for a function to compute the so-called phylogenetic variance-covariance matrix. That is, the N×N matrix (for N species) containing in position i,j the height above the root of the MRCA of species i and j (as done in the ape function vcv.phylo) - but including ancestral nodes. As it turns out, I have a code fragment to do this in my function anc.trend() (link here). I pasted the fragment into a function vcvPhylo() as follows:

vcvPhylo<-function(phy,anc.nodes=T){
  if(anc.nodes){
    D<-dist.nodes(phy)
    ntips<-length(phy$tip.label)
    Cii<-D[ntips+1,]
    C<-D; C[,]<-0
    for(i in 1:nrow(D)) for(j in 1:ncol(D)) 
      C[i,j]<-(Cii[i]+Cii[j]-D[i,j])/2
    dimnames(C)[[1]][1:length(phy$tip)]<-dimnames(C)[[2]][1:length(phy$tip)]<-phy$tip.label
    C<-C[c(1:ntips,(ntips+2):nrow(C)),c(1:ntips,(ntips+2):ncol(C))]
    return(C)
  } else return(vcv.phylo(phy))
}

and it seems to do the trick.