ダメ出し：三角ダイアグラムを描く（Drawing ternary diagrams）

http://staff.aist.go.jp/a.noda/programs/ternary/ternary.html

ベクトル化されている関数はその機能を使おう。短いプログラムは分かりやすいとは限らないが，長いプログラムが分かりづらいのは確か。

####################
#
# Ternary.plot
#
# Drawing ternary plot of compositional data
# (modified from triangle.plot in ade4)
#
# Noda Atsushi
# Last updated: 2008/09/17 16:27:46.

# 0.1.2: axes. 軸（三角形の辺）を描くかどうかを選択できるようにした．
# 0.1.1: プロットされる点の型 (pch) と色(col) と大きさ (cex) を選択できるようにした．
# 0.1: just beggining.

ternary.plot <- function (dat, label=rownames(dat), axes=TRUE, clabel=0,
    cpoint=1, cline=0, lwd=1, addline=TRUE, label.vertices=TRUE,
    cline.col="black", pch=par("pch"), col=par("col"), cex=par("cex")) {
 
    par.orig <- par(mar=par("mar"))
    on.exit(par(par.orig))

    nam <- colnames(dat)
    dat <- dat/rowSums(dat)
   
    xy <- cbind((dat[,1]+2*dat[,3])/2, sqrt(3)*dat[,1]/2)

    plot(0, 0, type="n", xlim=c(0, 1), ylim=c(-0.1, 1), xlab="",
         ylab="", xaxt="n", yaxt="n", asp=1, frame.plot=FALSE)
   
    if (addline) {
        int <- 10
        i <- 1:(int-1)/int
        i.2 <- i/2
        sqrt3.2.i <- sqrt(3)*i.2
        segments(1-i.2,  sqrt3.2.i,  1-i, 0, col="lightgrey")
        segments(i.2,  sqrt3.2.i,  i, 0, col="lightgrey")
        segments(1-i.2,  sqrt3.2.i,  i.2, sqrt3.2.i, col="lightgrey")
    }

    if (axes)
        polygon(c(0.5, 0, 1), c(sqrt3.2, 0, 0), lwd=1)
    if (label.vertices)
        text(c(0.5, 0, 1), c(sqrt3.2, 0, 0), labels=nam, cex=1.5, pos=c(3, 1, 1))
    if (cpoint > 0)
        points(xy, pch=pch, cex=cex*cpoint, col=col)
    if (cline > 0)
        lines(xy, lwd=1, col=cline.col)
    if (clabel > 0)
        text(xy[, 1], xy[, 2], label, clabel)

    return(invisible(xy))
}