plot(pressure)
text(150, 600, 
     "Pressure (mm Hg)\nversus\nTemperature (Celsius)")





#
#  Comment:
# 
#  Examples of the use of standard high-level plotting functions.
# 
#  In each case, extra output is also added using low-level 
#  plotting functions.
#


par(mfrow=c(3, 2))

# Scatterplot
x <- c(0.5, 2, 4, 8, 12, 16)
y1 <- c(1, 1.3, 1.9, 3.4, 3.9, 4.8)
y2 <- c(4, .8, .5, .45, .4, .3)
par(las=1, mar=c(4, 4, 2, 4))
plot.new()
plot.window(range(x), c(0, 6))
lines(x, y1)
lines(x, y2)
points(x, y1, pch=16, cex=2)
points(x, y2, pch=21, bg="white", cex=2)
par(col="grey50", fg="grey50", col.axis="grey50")
axis(1, at=seq(0, 16, 4))
axis(2, at=seq(0, 6, 2))
axis(4, at=seq(0, 6, 2))
box(bty="u")
mtext("Travel Time (s)", side=1, line=2, cex=0.8)
mtext("Responses per Travel", side=2, line=2, las=0, cex=0.8)
mtext("Responses per Second", side=4, line=2, las=0, cex=0.8)
text(4, 5, "Bird 131")
par(mar=c(5.1, 4.1, 4.1, 2.1), col="black", fg="black", col.axis="black")

# Histogram
# Random data
Y <- rnorm(50)
# Make sure no Y exceed [-3.5, 3.5]
Y[Y < -3.5 | Y > 3.5] <- NA
x <- seq(-3.5, 3.5, .1)
dn <- dnorm(x)
par(mar=c(4.5, 4.1, 3.1, 0))
hist(Y, breaks=seq(-3.5, 3.5), ylim=c(0, 0.5), 
     col="grey80", freq=FALSE)
lines(x, dnorm(x), lwd=2)
par(mar=c(5.1, 4.1, 4.1, 2.1))

# Barplot
# Modified from example(barplot)
par(mar=c(2, 3.1, 2, 2.1))
midpts <- barplot(VADeaths, col=grey(0.5 + 1:5/12), 
                  names=rep("", 4))
mtext(sub(" ", "\n", colnames(VADeaths)),
      at=midpts, side=1, line=0.5, cex=0.5)
text(rep(midpts, each=5), apply(VADeaths, 2, cumsum) - VADeaths/2,
     VADeaths, col=rep(c("white", "black"), times=2:3, cex=0.8))
par(mar=c(5.1, 4.1, 4.1, 2.1))

# Boxplot
# Modified example(boxplot) - itself from suggestion by Roger Bivand
par(mar=c(3, 4.1, 2, 0))
     boxplot(len ~ dose, data = ToothGrowth,
             boxwex = 0.25, at = 1:3 - 0.2,
             subset= supp == "VC", col="grey90",
             xlab="",
             ylab="tooth length", ylim=c(0,35))
     mtext("Vitamin C dose (mg)", side=1, line=2.5, cex=0.8)
     boxplot(len ~ dose, data = ToothGrowth, add = TRUE,
             boxwex = 0.25, at = 1:3 + 0.2,
             subset= supp == "OJ", col="grey70")
     legend(1.5, 9, c("Ascorbic acid", "Orange juice"), bty="n",
            fill = c("grey90", "grey70"))
par(mar=c(5.1, 4.1, 4.1, 2.1))

# Persp
# Almost exactly example(persp)
    x <- seq(-10, 10, length= 30)
     y <- x
     f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
     z <- outer(x, y, f)
     z[is.na(z)] <- 1
# 0.5 to include z axis label
par(mar=c(0, 0.5, 0, 0), lwd=0.1)
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "grey80")
par(mar=c(5.1, 4.1, 4.1, 2.1), lwd=1)

# Piechart
# Example 4 from help(pie)
par(mar=c(0, 2, 1, 2), xpd=FALSE, cex=0.5)
     pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12)
     names(pie.sales) <- c("Blueberry", "Cherry",
         "Apple", "Boston Cream", "Other", "Vanilla")
     pie(pie.sales, col = gray(seq(0.4,1.0,length=6)))




#
# Comment:
#
# A sophisticated example of adding further output to a basic plot.
# 
# Most of the functions defined are just for calculating values
# relevant to the data analysis.  
# 
# The function plotPars() is the one of interest for seeing how
# the drawing of the plot is done.
#


params <- function(N, breaks, p=seq(0.001, 1, length=100)) {
  list(N=N, T=1/breaks, p=p, q=1-p)
}

pdfcomp <- function(comp, params) {
  n <- params$T
  p <- params$p
  q <- params$q
  y <- round(comp/n)
  choose(n, comp)*p^comp*q^(n-comp) / (1 - q^n)
}

# Expected num sherds (for a vessel) [=completeness]
expcomp <- function(params) {
  params$T*params$p/(1-params$q^params$T)
}

# Variance of num sherds (for a vessel)
varcomp <- function(params) {
  n <- params$T
  p <- params$p
  q <- params$q
  # From Johnson & Kotz
  (n*p*q / (1 - q^n)) - (n^2*p^2*q^n / (1 - q^n)^2)
  # n^2 times Thomas Yee's formula
  # n^2*((p*(1 + p*(n - 1)) / (n*(1 - q^n))) - (p^2 / (1 - q^n)^2))
}

# Expected value of completeness (for a sample of vessels)
expmeancomp <- function(params) {
  expcomp(params)
}

# Variance of completeness (for a sample of vessels)
# Use the expected number of vessels in sample as denominator
varmeancomp <- function(params) {
  varcomp(params)/(numvess(params))
}

numvess <- function(params) {
  params$N*(1-params$q^params$T)
}

ecomp <- function(p, T, comp) {
  q <- 1 - p
  T*p/(1 - q^T) - comp
}

estN <- function(comp, broke, n) {
  T <- 1/broke
  n / (1 - (1 - uniroot(ecomp, c(0.00001, 1), T=T, comp=comp)$root)^T)
}

nvessscale <- function(params, xlim, ylim, new=TRUE) {
  if (new)
    par(new=TRUE)
  plot(0:1, c(1, params$N), type="n", axes=!new, ann=FALSE,
       xlim=xlim, ylim=ylim)
}

compscale <- function(params, xlim, ylim, new=TRUE) {
  if (new)
    par(new=TRUE)
  plot(0:1, c(1, params$T), type="n", axes=!new, ann=FALSE,
       xlim=xlim, ylim=ylim)
}

lowerCI <- function(p, N, breaks, lb) {
  params <- params(N, breaks, p)
  expmeancomp(params) - 2*sqrt(varmeancomp(params)) - lb
}

upperCI <- function(p, N, breaks, lb) {
  params <- params(N, breaks, p)
  expmeancomp(params) + 2*sqrt(varmeancomp(params)) - lb
}

critP <- function(comp, params) {
  c(uniroot(lowerCI, c(0.00001, 1), N=params$N,
            breaks=1/params$T, lb=max(comp))$root,
    if (upperCI(0.00001, params$N, 1/params$T, min(comp)) > 0) 0
    else uniroot(upperCI, c(0.00001, 1), N=params$N,
                 breaks=1/params$T, lb=min(comp))$root)
}

anncomp <- function(params, comp, xlim, ylim, cylim) {
  cp <- critP(comp, params)
  nv <- numvess(params(params$N, 1/params$T, cp))
  nvessscale(params, xlim, ylim)
  polygon(c(cp[2], cp[2], 0, 0, cp[1], cp[1]),
          c(0, nv[2], nv[2], nv[1], nv[1], 0),
          col="grey90", border=NA)
  text(0, nv[1], paste(round(nv[1]),
                       " (", round(100*nv[1]/params$N), "%)", sep=""),
       adj=c(0, 0), col="grey")
  text(0, nv[2], paste(round(nv[2]), 
                       " (", round(100*nv[2]/params$N), "%)", sep=""),
       adj=c(0, 1), col="grey")
  compscale(params, xlim, cylim)
  segments(1, min(comp), cp[2], comp, col="grey")
  segments(1, max(comp), cp[1], comp, col="grey")
  text(1, comp, paste(comp, collapse="-"), adj=c(1, 0), col="grey")
}

plotPars <- function(params, comp, xlim=NULL, ylim=NULL) {
  mean <- expmeancomp(params)
  var <- 2*sqrt(varmeancomp(params))
  lb <- mean - var
  ub <- mean + var
  par(mar=c(5, 4, 4, 4))
  if (is.null(ylim))
    cylim <- ylim
  else
    cylim <- c(1 + ((ylim[1] - 1)/(params$N - 1))*(params$T - 1),
               1 + ((ylim[2] - 1)/(params$N - 1))*(params$T - 1))
  nvessscale(params, xlim, ylim, new=FALSE)
  compscale(params, xlim, cylim)
  polygon(c(params$p, rev(params$p)), c(lb, rev(ub)),
          col="grey90", border=NA)
  anncomp(params, comp, xlim, ylim, cylim)
  nvessscale(params, xlim, ylim)
  mtext("Number of Vessels", side=2, line=3)
  mtext("Sampling Fraction", side=1, line=3)
  lines(params$p, numvess(params))
  par(new=TRUE)
  compscale(params, xlim, cylim)
  mtext("Completeness", side=4, line=3)
  axis(4)
  lines(params$p, mean, lty="dashed")
  lines(params$p, lb, lty="dotted")
  lines(params$p, ub, lty="dotted")
  mtext(paste("N = ", round(params$N),
              "     brokenness = ", round(1/params$T, 3), sep=""),
        side=3, line=2)
}

par(cex=0.8, mar=c(3, 3, 3, 3))
p6 <- params(estN(1.2, 0.5, 200), 0.5)
plotPars(p6, 1.2)
nvessscale(p6, NULL, NULL)
pcrit <- 1 - (1 - 200/estN(1.2, 0.5, 200))^(1/p6$T)
lines(c(0, pcrit), c(200, 200))
lines(c(pcrit, pcrit), c(200, 0))



#
# Comment:
#
# A bit of mucking around is required to get the second (whole-world)
# map positioned correctly;  this provides an example of calling a 
# plotting function to perform calculations but do no drawing (see the
# second call to the map() function).
#
# Makes use of the "maps" and "mapproj" packages to draw the maps.
#


library(maps)
par(mar=rep(0, 4))
map("nz", fill=TRUE, col="grey80")
points(174.75, -36.87, pch=16, cex=2)
arrows(172, -36.87, 174, -36.87, lwd=3)
text(172, -36.87, "Auckland  ", adj=1, cex=2)
# mini world map as guide
maplocs <- map(projection="sp_mercator", wrap=TRUE, lwd=0.1, 
               col="grey", ylim=c(-60, 75),
               interior=FALSE, orientation=c(90, 180, 0), add=TRUE,
               plot=FALSE)
xrange <- range(maplocs$x, na.rm=TRUE)
yrange <- range(maplocs$y, na.rm=TRUE)
aspect <- abs(diff(yrange))/abs(diff(xrange))
# customised to 6.5 by 4.5 figure size
par(fig=c(0.99 - 0.5, 0.99, 0.01, 0.01 + 0.5*aspect*4.5/6.5), 
    mar=rep(0, 4), new=TRUE)
plot.new()
plot.window(xlim=xrange,
            ylim=yrange)
map(projection="sp_mercator", wrap=TRUE, lwd=0.1, ylim=c(-60, 75),
    interior=FALSE, orientation=c(90, 180, 0), add=TRUE)
symbols(-.13, -0.8, circles=1, inches=0.1, add=TRUE)



#
# Comment:
#
# Code by Arden Miller (Department of Statistics, The University of Auckland).
# 
# Lots of coordinate transformations being done "by hand".
# This code is not really reusable;  just a demonstration that very 
# pretty results are possible if you're sufficiently keen.
#


par(mfrow=c(2, 1), pty="s", mar=rep(1, 4)) 
# Create plotting region and plot outer circle
plot(c(-1.1, 1.2), c(-1.1, 1.2),
     type="n", xlab="", ylab="", 
     xaxt="n", yaxt="n", cex.lab=2.5)
angs <- seq(0, 2*pi, length=500)
XX <- sin(angs)
YY <- cos(angs)
lines(XX, YY, type="l")

# Set constants
phi1 <- pi*2/9
k1 <- sin(phi1)
k2 <- cos(phi1)

# Create grey regions
obsphi <- pi/12
lambdas <- seq(-pi, pi, length=500)
xx <- cos(pi/2 - obsphi)*sin(lambdas)
yy <- k2*sin(pi/2 - obsphi)-k1 * cos(pi/2 - obsphi)*cos(lambdas)
polygon(xx, yy, col="grey")
lines(xx, yy, lwd=2)
theta1sA <- seq(-obsphi, obsphi, length=500)
theta2sA <- acos(cos(obsphi)/cos(theta1sA))
theta1sB <- seq(obsphi, -obsphi, length=500)
theta2sB <-  -acos(cos(obsphi)/cos(theta1sB))
theta1s <- c(theta1sA, theta1sB)
theta2s <- c(theta2sA, theta2sB)
xx <- cos(theta1s)*sin(theta2s+pi/4)
yy <- k2*sin(theta1s)-k1*cos(theta1s)*cos(theta2s+pi/4)
polygon(xx, yy, col="grey")
lines(xx, yy, lwd=2)
xx <- cos(theta1s)*sin(theta2s-pi/4)
yy <- k2*sin(theta1s)-k1*cos(theta1s)*cos(theta2s-pi/4)
polygon(xx, yy, col="grey")
lines(xx, yy, lwd=2)

# Plot longitudes
vals <- seq(0, 7, 1)*pi/8
for(lambda in vals){
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, pi+phi, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=.5)
}

# Grey out polar cap
phi <- 5.6*pi/12
lambdas <- seq(-pi, pi, length=500)
xx <- cos(phi)*sin(lambdas)
yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
polygon(xx, yy, col="grey")

# Plot Latitudes
vals2 <- seq(-2.8, 5.6, 1.4)*pi/12
for(phi in vals2){
  if (k1*sin(phi) > k2 * cos(phi)) 
    crit <- pi 
  else 
    crit <- acos((-k1*sin(phi))/(k2*cos(phi)))
  lambdas <- seq(-crit, crit, length=500)
  xx <- cos(phi)*sin(lambdas)
  yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
  lines(xx, yy, lwd=.5)
}


# Plots axes and label
lines(c(0.00, 0.00), c(k2*sin(pi/2), 1.11), lwd=4)
lines(c(0.00, 0.00), c(-1, -1.12), lwd=4)
a2x <- sin(-pi/4)
a2y <- cos(-pi/4)*(-k1)
lines(c(a2x, 1.5*a2x), c(a2y, 1.5*a2y), lwd=4)
k <- sqrt(a2x^2+a2y^2)
lines(c(-a2x/k, 1.2*(-a2x/k)), c(-a2y/k, 1.2*(-a2y/k)), lwd=4)
a3x <- sin(pi/4)
a3y <- cos(pi/4)*(-k1)
lines(c(a3x, 1.5*a3x), c(a3y, 1.5*a3y), lwd=4)
k <- sqrt(a3x^2+a3y^2)
lines(c(-a3x/k, 1.2*(-a3x/k)), c(-a3y/k, 1.2*(-a3y/k)), lwd=4)
text(0.1, 1.12, expression(bold(X[1])))
text(-1.07, -.85, expression(bold(X[2])))
text(1.11, -.85, expression(bold(X[3])))

# set plot region and draw outer circle
plot(c(-1.1, 1.2),  c(-1.1, 1.2),
     type="n", xlab="", ylab="", 
     xaxt="n", yaxt="n", cex.lab=2.5)
angs <- seq(0, 2*pi, length=500)
XX <- sin(angs)
YY <- cos(angs)
lines(XX, YY, type="l")

# set constants
phi1 <- pi*2/9
k1 <- sin(phi1)
k2 <- cos(phi1)
obsphi <- pi/24

# create X2X3 grey region and plot boundary
crit <- acos((-k1*sin(obsphi))/(k2 * cos(obsphi)))
lambdas <- seq(-crit, crit, length=500)
xx1 <- cos(obsphi)*sin(lambdas)
yy1 <- k2*sin(obsphi)-k1 * cos(obsphi)*cos(lambdas)
obsphi <-  -pi/24
crit <- acos((-k1*sin(obsphi))/(k2 * cos(obsphi)))
lambdas <- seq(crit, -crit, length=500)
xx3 <- cos(obsphi)*sin(lambdas)
yy3 <- k2*sin(obsphi)-k1 * cos(obsphi)*cos(lambdas)
ang1 <-  atan(xx1[500]/yy1[500])
ang2 <- pi+atan(xx3[1]/yy3[1])
angs <- seq(ang1, ang2, length=50)
xx2 <- sin(angs)
yy2 <- cos(angs)
ang4 <-  atan(xx1[1]/yy1[1])
ang3 <-  -pi+ atan(xx3[500]/yy3[500])
angs <- seq(ang3, ang4, length=50)
xx4 <- sin(angs)
yy4 <- cos(angs)
xxA <- c(xx1, xx2, xx3, xx4)
yyA <- c(yy1, yy2, yy3, yy4)
polygon(xxA, yyA, border="grey", col="grey")
xx1A <- xx1
yy1A <- yy1
xx3A <- xx3
yy3A <- yy3

# create X1X3 grey region and plot boundary
obsphi <- pi/24
crit <- pi/2-obsphi
theta1sA <- c(seq(-crit, crit/2, length=200), seq(crit/2, crit, length=500))
theta2sA <- asin(cos(crit)/cos(theta1sA))
theta1sB <- seq(crit, crit/2, length=500)
theta2sB <-  pi-asin(cos(crit)/cos(theta1sB))
theta1s <- c(theta1sA, theta1sB)
theta2s <- c(theta2sA, theta2sB)
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s+pi/4)
xx1 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]+pi/4)
yy1 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]+pi/4)
theta2s <-  -theta2s
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s+pi/4)
xx3 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]+pi/4)
yy3 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]+pi/4)
rev <- seq(length(xx3), 1, -1)
xx3 <- xx3[rev]
yy3 <- yy3[rev]
ang1 <-  pi+atan(xx1[length(xx1)]/yy1[length(yy1)])
ang2 <-  pi+atan(xx3[1]/yy3[1])
angs <- seq(ang1, ang2, length=50)
xx2 <- sin(angs)
yy2 <- cos(angs)
ang4 <-  pi+atan(xx1[1]/yy1[1])
ang3 <-  pi+atan(xx3[length(xx3)]/yy3[length(yy3)])
angs <- seq(ang3, ang4, length=50)
xx4 <- sin(angs)
yy4 <- cos(angs)
xxB <- c(xx1, -xx2, xx3, xx4)
yyB <- c(yy1, -yy2, yy3, yy4)
polygon(xxB, yyB, border="grey", col="grey")
xx1B <- xx1
yy1B <- yy1
xx3B <- xx3
yy3B <- yy3

# create X1X2 grey region and plot boundary
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s-pi/4)
xx1 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]-pi/4)
yy1 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]-pi/4)
theta2s <-  -theta2s
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s-pi/4)
xx3 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]-pi/4)
yy3 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]-pi/4)
rev <- seq(length(xx3), 1, -1)
xx3 <- xx3[rev]
yy3 <- yy3[rev]
ang1 <-  pi+atan(xx1[length(xx1)]/yy1[length(yy1)])
ang2 <-  pi+atan(xx3[1]/yy3[1])
angs <- seq(ang1, ang2, length=50)
xx2 <- sin(angs)
yy2 <- cos(angs)
ang4 <-  pi+atan(xx1[1]/yy1[1])
ang3 <-  pi+atan(xx3[length(xx3)]/yy3[length(yy3)])
angs <- seq(ang3, ang4, length=50)
xx4 <- sin(angs)
yy4 <- cos(angs)
xx <- c(xx1, -xx2, xx3, xx4)
yy <- c(yy1, -yy2, yy3, yy4)
polygon(xx, yy, border="grey", col="grey")
xx1C <- xx1
yy1C <- yy1
xx3C <- xx3
yy3C <- yy3


# plot boundaries to grey regions
lines(xx1C[2:45], yy1C[2:45], lwd=2)
lines(xx1C[69:583], yy1C[69:583], lwd=2)
lines(xx1C[660:1080], yy1C[660:1080], lwd=2)
lines(xx3C[13:455], yy3C[13:455], lwd=2)
lines(xx3C[538:1055], yy3C[538:1055], lwd=2)
lines(xx3C[1079:1135], yy3C[1079:1135], lwd=2)
lines(xx1A[6:113], yy1A[6:113], lwd=2)
lines(xx1A[153:346], yy1A[153:346], lwd=2)
lines(xx1A[389:484], yy1A[389:484], lwd=2)
lines(xx3A[1:93], yy3A[1:93], lwd=2)
lines(xx3A[140:362], yy3A[140:362], lwd=2)
lines(xx3A[408:497], yy3A[408:497], lwd=2)
lines(xx1B[2:45], yy1B[2:45], lwd=2)
lines(xx1B[69:583], yy1B[69:583], lwd=2)
lines(xx1B[660:1080], yy1B[660:1080], lwd=2)
lines(xx3B[13:455], yy3B[13:455], lwd=2)
lines(xx3B[538:1055], yy3B[538:1055], lwd=2)
lines(xx3B[1079:1135], yy3B[1079:1135], lwd=2)

# Plot longitudes
vals <- seq(-7, 8, 1)*pi/8
for(lambda in vals){
  sl <- sin(lambda)
  cl <- cos(lambda)
  phi <- atan(((0-1)*k2*cl)/(k1))
  angs <- seq(phi, 5.6*pi/12, length=500)
  xx <- cos(angs)*sl
  yy <- k2*sin(angs)-k1*cos(angs)*cl
  lines(xx, yy, lwd=.5)
}


# Plot Latitudes
# vals2 <- seq(-2.8, 5.6, 1.4)*pi/12
vals2 <- c(-1.5, 0, 1.5, 3.0, 4.5, 5.6)*pi/12
for(phi in vals2){
  if (k1*sin(phi) > k2 * cos(phi)) 
    crit <- pi 
  else 
    crit <- acos((-k1*sin(phi))/(k2*cos(phi)))
  lambdas <- seq(-crit, crit, length=500)
  xx <- cos(phi)*sin(lambdas)
  yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
  lines(xx, yy, lwd=.5)
}


# create lines for X1X2- and X1X3-planes
lambda <- pi/4
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, pi+phi, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=2)
lambda <- 3*pi/4
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, pi+phi, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=2)

# create line for X2X3-plane
phi <- 0
crit <- acos((-k1*sin(phi))/(k2 * cos(phi)))
lambdas <- seq(-crit, crit, length=500)
xx <- cos(phi)*sin(lambdas)
yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
lines(xx, yy, lwd=2)

# create axes
lines(c(0.00, 0.00), c(k2*sin(pi/2), 1.11), lwd=4)
lines(c(0.00, 0.00), c(-1, -1.12), lwd=4)
a2x <- sin(-pi/4)
a2y <- cos(-pi/4)*(-k1)
lines(c(a2x, 1.5*a2x), c(a2y, 1.5*a2y), lwd=4)
a3x <- sin(pi/4)
a3y <- cos(pi/4)*(-k1)
lines(c(a3x, 1.5*a3x), c(a3y, 1.5*a3y), lwd=4)
k <- sqrt(a3x^2+a3y^2)
lines(c(-a3x/k, 1.2*(-a3x/k)), c(-a3y/k, 1.2*(-a3y/k)), lwd=4)
k <- sqrt(a2x^2+a2y^2)
lines(c(-a2x/k, 1.2*(-a2x/k)), c(-a2y/k, 1.2*(-a2y/k)), lwd=4)


# add text
text(-1.07, -.85, expression(bold(X[2])))
text(1.11, -.85, expression(bold(X[3])))
text(0.1, 1.12, expression(bold(X[1])))

lines(XX, YY, type="l")




#
# Comment:
#
# An example of a one-off image drawn using the grid system.
#
# The code is somewhat modular and general, with functions
# for producing different shapes, but the sizes and
# locations used in this particular image assume a 2:1 aspect ratio.
#
# The gradient-fill background (dark at the top to lighter at the 
# bottom) is achieved by filling multiple overlapping polygons with
# slowly changing shades of grey.
#

library(gridBase)
pushViewport(viewport(xscale=c(0, 1), yscale=c(0.5, 1),
             clip=TRUE))
             
res <- 50
for (i in 1:res)
  grid.rect(y=1 - (i-1)/res, just="top",
            gp=gpar(col=NULL, fill=grey(0.5*i/res)))

moon <- function(x, y, size) {
  angle <- seq(-90, 90, length=50)/180*pi
  x1 <- x + size*cos(angle)
  y1 <- y + size*sin(angle)
  mod <- 0.8
  x2 <- x + mod*(x1 - x)
  grid.polygon(c(x1, rev(x2)), c(y1, rev(y1)),
               default.unit="native",
               gp=gpar(col=NULL, fill="white"))
}

moon(.1, .9, .03)

star <- function(x, y, size) {
  x1 <- c(x,           x + size*.1, x + size*.5, x + size*.1,
          x,           x - size*.1, x - size*.5, x - size*.1) + .05
  y1 <- c(y - size,    y - size*.1, y,           y + size*.1,
          y + size*.7, y + size*.1, y,           y - size*.1) + .05
  grid.polygon(x1, y1, 
               default.unit="native",
               gp=gpar(col=NULL, fill="white"))
}

star(.5, .7, .02)
star(.8, .9, .02)
star(.72, .74, .02)
star(.62, .88, .02)

grid.circle(runif(20, .2, 1), runif(20, .6, 1), r=.002,
            default.unit="native",
            gp=gpar(col=NULL, fill="white"))

hill <- function(height=0.1, col="black") {
  n <- 100
  x <- seq(0, 1, length=n)
  y1 <- sin(runif(1) + x*2*pi)
  y2 <- sin(runif(1) + x*4*pi)
  y3 <- sin(runif(1) + x*8*pi)
  y <- 0.6 + height*((y1 + y2 + y3)/3)
  grid.polygon(c(x, rev(x)), c(y, rep(0, n)),
               default.unit="native",
               gp=gpar(col=NULL, fill=col))
}

hill()

rdir <- function(n) {
  sample(seq(-45, 45, length=10), n)/180*pi
}

grid.text("Once upon a time ...",
          x=.15, y=.51, just="bottom",
          default.unit="native",
          gp=gpar(col="white", fontface="italic", fontsize=10))

popViewport()

grid.rect()



##new plot

par(mfrow=c(2, 2))
z <- 2 * volcano        # Exaggerate the relief
x <- 10 * (1:nrow(z))   # 10 meter spacing (S to N)
y <- 10 * (1:ncol(z))   # 10 meter spacing (E to W)
# Don't draw the grid lines :  border = NA
par(mar=rep(0, 4))
persp(x, y, z, theta = 135, phi = 30, col = "light grey", scale = FALSE,
      ltheta = -120, shade = 0.75, border = NA, box = FALSE)
mtext("persp()", side=3, line=-2)
par(mar=c(3, 3, 2, 0.5))
# Note that  example(trees)  shows more sensible plots!
N <- nrow(trees)
attach(trees)
# Girth is diameter in inches
symbols(Height, Volume, circles=Girth/24, inches=FALSE,
        main="", xlab="", ylab="", bg=grey(Girth/max(Girth)))
mtext("symbols()", side=3, line=0.5)
par(mar=rep(0.5, 4))
contour(x, y, z, asp=1, labcex=0.35, axes=FALSE)
rect(0, 0, 870, 620)
mtext("contour()", side=3, line=-1.5)
image(x, y, z, asp=1, col=grey(0.5 + 1:12/24), xlab="", ylab="", axes=FALSE)
rect(min(x)-5, min(y)-5, max(x)+5, max(y)+5)
mtext("image()", side=3, line=-1.5)



par(mar=rep(0, 4), lwd=0.1)
z <- 2 * volcano        
x <- 10 * (1:nrow(z))   
y <- 10 * (1:ncol(z))   
trans <- persp(x, y, z, theta = 135, phi = 30, 
               scale = FALSE, ltheta = -120, 
               # shade=0.5, border=NA, 
               box = FALSE)
box(col="grey", lwd=1)

trans3d <- function(x,y,z,pmat) {
  tmat <- cbind(x,y,z,1)%*% pmat
  tmat[,1:2] / tmat[,4]
}

summit <- trans3d(x[20], y[31], max(z), trans)
points(summit[1], summit[2], pch=16)
summitlabel <- trans3d(x[20], y[31], max(z) + 50, trans)
text(summitlabel[1], summitlabel[2], "Summit")

drawRoad <- function(x, y, z, trans) {
  road <- trans3d(x, y, z, trans)
  lines(road[,1], road[,2], lwd=5)
  lines(road[,1], road[,2], lwd=3, col="grey")
}
with(volcano.summitRoad,
     drawRoad(srx, sry, srz, trans))
with(volcano.upDownRoad,
     {
       clipudx <- udx
       clipudx[udx < 230 & udy < 300 | 
               udx < 150 & udy > 300] <- NA
       drawRoad(clipudx, udy, udz, trans)
     })
with(volcano.accessRoad,
     drawRoad(arx, ary, arz, trans))




grid.rect(gp=gpar(col="grey"))
pushViewport(viewport(w=0.8, h=0.8))
offset <- unit(3, "mm")
grid.text("CIE-LUV\nchroma = 40\nluminance = 80",
          gp=gpar(fontfamily="mono"))
t <- seq(0, 2*pi, length=7)[-7]
x <- 0.5 + .4*cos(t)
y <- 0.5 + .4*sin(t)
rad <- .13
cols <- hcl(t/pi*180, 40, 80)
grid.circle(x, y, r=unit(rad, "npc"),
            gp=gpar(fill=cols))
labels <- c("r", "g", "b")
rgbcols <- col2rgb(cols)
hbars <- seq(50, 250, 50)
for (i in 1:6) {
  pushViewport(viewport(x[i], y[i], width=.2, height=.2,
                        layout=grid.layout(2, 3,
                          heights=unit(c(1, 1), c("null", "lines")))))
  for (j in 1:3) {
    pushViewport(viewport(layout.pos.col=j,
                          layout.pos.row=1))
    grid.text(labels[j], y=unit(-.5, "lines"),
              gp=gpar(fontfamily="mono"))
    grid.lines(x=c(.1, .9), y=0)
    grid.rect(y=0, width=.6, height=rgbcols[j, i]/300,
              just="bottom",
              gp=gpar(fill=NA))
    popViewport()
  }
  popViewport()
}
popViewport()
