#--Read in table of data (from Sanderson et al. 2006):
# This refers to a sample of 20 clusters of galaxies with Chandra X-ray data
file <- "http://www.sr.bham.ac.uk/~ajrs/papers/sanderson06/table1.txt"

#--By default, read.table() will treat non-numeric values as factors:
# Sometimes this is annoying, in which case use "as.is=T"; you can specify the
# type of all the input columns explicitly -- see ?read.table
A <- read.table(file,  header=TRUE, sep="|")

#--Rename factor levels:
levels(A$cctype) <- c("Cool core clusters", "Non-cool core clusters")

#--Source errobar script:
source("http://www.sr.bham.ac.uk/~ajrs/R/scripts/errorbar.R")

#--Change output device to pdf file:
pdf(file="cluster_z_vs_kT.pdf", height=5.5, width=8.25, pointsize=16)

#---Remove extra top margin:
par(mar=c(3, 3, 0.5, 1))  # Trim margin around plot [bottom, left, top, right]
par(tcl=0.35)  # Switch tick marks to insides of axes
par(mgp=c(1.5, 0.2, 0))  # Set margin lines; default c(3, 1, 0) [title,labels,line]
par(xaxs="r", yaxs="r")  # Extend axis limits by 4% ("i" does no extension)
par(lwd=2)

#--Define different colours & point styles:
cols <- c("red", "blue")
pchs <- c(20, 15)	# 20=Small filled circle; 15=filled square

#--Add in columns with colour & point style defined as per cool-core type:
A <- transform(A, col=cols[factor(cctype)], 
                  pch=pchs[factor(cctype)])

#--Convert colour from a factor to a character vector:
A$col <- as.character(A$col)

print(class(A$col))

#--Now plot mean temperature vs. redshift, with colours & point styles 
#   appropriate for CC status:
plot(Tx ~ z, data=A, xlab="Redshift", ylab="Temperature (keV)", col=col,
     pch=pch, log="xy")

#--Add errorbars (using function in sourced errorbar.R script):
with(A, errorbar(z, 0, 0, Tx, -Tx.le, Tx.ue, bnd=F, col=col))	# bnd=F means input data are errors not confidence bounds

#--Add a plot legend:
legend(x="topleft", c(levels(A$cctype)), col=cols, text.col=cols, 
       inset=0.05, pch=pchs, bty="n")

#--Ensure graphics device finishes cleanly:
dev.off()