# # 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 gray 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="gray") 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="gray") 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="gray") 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="gray") # 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 gray 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="gray", col="gray") xx1A <- xx1 yy1A <- yy1 xx3A <- xx3 yy3A <- yy3 # create X1X3 gray 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="gray", col="gray") xx1B <- xx1 yy1B <- yy1 xx3B <- xx3 yy3B <- yy3 # create X1X2 gray 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="gray", col="gray") xx1C <- xx1 yy1C <- yy1 xx3C <- xx3 yy3C <- yy3 # plot boundaries to gray 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")