# R functions (see <http://www.r-project.org/>)

# readCgats v3
readCgats <- function(file, dec=".") {
    tmp  <- scan(file, n=2000, what="character", comment.char="#", quiet = TRUE)
    tmp2 <- charmatch(c("BEGIN_DATA_FORMAT", "NUMBER_OF_FIELDS", "NUMBER_OF_SETS","END_DATA_FORMAT"), tmp)
    bdf  <- tmp2[1]+1
    edf  <- as.integer(tmp[tmp2[2]+1])+bdf-1
    if(!is.finite(edf)) {
        edf <- tmp2[4]-1}
    dataformat <- tmp[bdf:edf]
    nos <- as.integer(tmp[tmp2[3]+1])
    
    tmp <- scan(file,n=500, what="character", sep="\n", blank.lines.skip = FALSE, quiet = TRUE)
    skiplines <- grep("BEGIN_DATA[[:space:]]*$", tmp)
    
    cgats <- read.table(file,dec = dec, col.names=dataformat, skip=skiplines, nrows=nos, header = FALSE)
}

# writeCgats v1
writeCgats <- function(dframe, file, id=FALSE, dec=".") {

    if(id==TRUE) {
        na <- c("SAMPLE_ID",names(dframe))
        nc <- ncol(dframe)+1
    } else {
        na <- names(dframe)
        nc <- ncol(dframe)
    }
        header <- list("R",
            c("CREATED", paste("\"", format(Sys.time(), usetz = TRUE), "\"", sep="")),
            c("NUMBER_OF_FIELDS",   nc),
            "BEGIN_DATA_FORMAT",
            na,
            "END_DATA_FORMAT",
            c("NUMBER_OF_SETS", nrow(dframe)),
            "BEGIN_DATA")
    
    write(header[[1]],file, sep="\t")
        for (i in 2:length(header))
            write(header[[i]],file, sep="\t", ncolumns=1000,append=TRUE)
        write.table(format(dframe, digits=5), file = file, append = TRUE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = dec, row.names = id,
            col.names = FALSE, qmethod = "escape")
        write("END_DATA", file, append=TRUE)
}

# planck v1
# spektrale spezifische Ausstrahlung M des schwarzen Strahlers in W/m2/m
# lambda in nm
# T in Kelvin
planck <- function(lambda,T) {
    
    lambda<-lambda*1e-9
    
    c <- 299792458  # Lichtgeschwindingkeit m/s
    h <- 6.62606896e-34 # Plancksches Wirkungsquantum Js
    k <- 1.38065042e-23 # Boltzmannkonstante J/K

    
    M <- 2*pi*h*c^2/(lambda)^5*1/(exp(h*c/(lambda*k*T))-1)
    M
}

# end R functions

# R console:

p <- data.frame(t(planck(38:78*10,1000)))
for(x in seq(from=1100, to=10000,by=100)) {
    q <- data.frame(t(planck(38:78*10,x)))
    p <- rbind(p,q)
}

names(p) <- paste("nm",38:78*10, sep="")
rownames(p) <- paste(seq(from=1000, to=10000, by=100), "K", sep="")

p<-data.frame(RGB_R=0, RGB_G=0, RGB_B=0, p)

writeCgats(p, "planck_1-10kK.txt",id=TRUE)

# bash / argyllCMS (see <http://argyllcms.com/>):

txt2ti3 -2 -d planck_1-10kK.txt planck_1-10kK
spec2cie planck_1-10kK.ti3 planck_1-10kK_s2c.ti3

# R console:

p2<-readCgats("planck_1-10kK_s2c.ti3")

p2$x<-p2$XYZ_X/(p2$XYZ_X+p2$XYZ_Y+p2$XYZ_Z)
p2$y<-p2$XYZ_Y/(p2$XYZ_X+p2$XYZ_Y+p2$XYZ_Z)

cie2deg<-read.table("cie2deg.table.txt",header=TRUE) 
# CIE 1931 Std. observer extracted from CIE Excel file "Selected Colorimetric Tables"
# <http://www.cie.co.at/index.php/index.php?i_ca_id=298>
# <http://files.cie.co.at/204.xls>
# <./cie2deg.table.txt>

cie2deg$x<-cie2deg$xbar/(cie2deg$xbar+cie2deg$ybar+cie2deg$zbar)
cie2deg$y<-cie2deg$ybar/(cie2deg$xbar+cie2deg$ybar+cie2deg$zbar)

plot(cie2deg$x, cie2deg$y, asp=1, xlab="x", ylab="y", pch="x", col="blue", cex=0.5)
abline(h=0:8/10, v=0:8/10, col="gray")
lines(cie2deg$x, cie2deg$y)
lines(cie2deg[c(1,81),"x"], cie2deg[c(1,81),"y"])
text(cie2deg[c(1,20:50,81),"x"], cie2deg[c(1,20:50,81),"y"], labels=cie2deg[c(1,20:50,81),"lambda"], cex=0.5, pos=2, col="blue")
lines(p2$x,p2$y,col="red")
points(p2[c(11,21,31,41,51,71,91),"x"],p2[c(11,21,31,41,51,71,91),"y"], pch="x",cex=0.5)
text(x=p2[c(11,21,31,41,51,71,91),"x"],y=p2[c(11,21,31,41,51,71,91),"y"], labels=p2[c(11,21,31,41,51,71,91),"SAMPLE_NAME"], pos=2, cex=0.5)

# plot saved as blackbody_xy.pdf