How to create a QQ plot of azimuths to test rotational symmetry of a spherical point dataset? I am trying to do some Q-Q plots and the Kuiper test for rotational symmetry about the mean direction.  These are unit vector, spherical data.  What I am struggling to understand is the rotation from my own directions, where φ is the azimuth direction (0 to 360) and θ is the dip up/down direction (0 to 180), to the rotated/modified φ' and θ'. I'll limit this posting to the Q-Q plot question.
Here is a small example set of 10 azimuth/φ angles: 
146.6
198.0
182.5
148.5
183.0
130.5
344.9
355.5
82.3
166.2
So, first the graphical method: Q-Q plot. Images, including #1 below, are from Fisher, Lewis and Embleton 1987 (reference below). I have the φ values, but the plot requires the rotated/modified φ' values.

Now below in Image #2 is the method for doing the rotation, equation 3.8.  I'm thinking I can set the ψ0 to zero, so I could use the "simpler" equation 3.9.

But equations 3.8 or 3.9 require equations 3.5 and 3.6, which are below in Images #3 and #4 respectively.  I'm able to calculate the Resultant length (R) and x-,y- and z-hat, so I have those, and which allows me to get φ-hat and φ-hat.  So, I'm not sure what to do with equation 3.9 or 3.9, or how to rotate my numbers.  Any suggestions?  Apologies for the long post of images.


Fisher, N.I., T. Lewis and B.J.J. Embleton. 1987. Statistical Analysis of Spherical Data. Cambridge, Cambridge University Press.
 A: You just want to study the azimuths of a set of spherical points $P_i$ relative to their spherical mean $\bar P$.  The most straightforward solution solves the spherical triangles $(N,\bar P, P_i)$ where $N$ is the North Pole.
Let the co-latitudes of the points $P_i$ and $\bar P$ (angles from North) be $a$ and $b$ respectively.  Let $\gamma$ be the angle between them: it's just the difference between the longitudes of the same two points.  The azimuth, with due East being zero and orienting angles counterclockwise, is determined by
$$\arctan_2(\sin(b)\cos(a) - \cos(b)\sin(a)\cos(\gamma),\  \sin(a)\sin(\gamma))$$
where $\arctan_2(y,x)$ is the angle of a point $(x,y)$ in the plane.  (This is supposed to be a numerically stable version of the formula, but I haven't tested it extensively.)
In this example, $100$ points were generated according to a (symmetric) Fisher-von Mises distribution distributed throughout the southern and western hemispheres, along with another $50$ points focused in the south and east.  The resulting distribution is not symmetric.

The mean point is shown as a red triangle.
Relative to the mean point, there is a cluster of points to its right (East) and upward (North), creating a swath of azimuths in the QQ plot between $0$ and $1$ (expressed in radians).  The diffuse cluster to its west creates a broader swath of azimuths between $3$ and $5$.  The QQ plot is clearly not uniform (for otherwise it would lie close to the dashed diagonal line), reflecting the bimodality of the spherical point distribution.

The R code that produced this example can be used to generate azimuthal QQ plots for any data.  It assumes the spherical coordinates are provided as rows in an array; the relevant rows are indexed by "phi" and "theta".
#
# Spherical triangle, two sides and included angle given.
# Returns the angle `alpha` opposite `a`, in radians between 0 and 2*pi.
#
SAS <- function(a, gamma, b) {
  atan2(sin(b)*cos(a) - cos(b)*sin(a)*cos(gamma), sin(a)*sin(gamma)) %% (2*pi)
}
#
# Cartesian coordinate conversion (for generating points).
#
xyz.to.spherical <- function(xyz) {
  xyz <- matrix(xyz, nrow=3)
  x <- xyz[1,]; y <- xyz[2,]; z <- xyz[3,]
  r2 <- x^2 + y^2
  rho <- sqrt(r2 + z^2)
  theta <- pi/2 - atan2(z, sqrt(r2))
  phi <- atan2(y, x)
  theta[x==0 && y==0] <- sign(z) * pi/2
  return (rbind(rho, theta, phi))
}
#
# Generate random points on the sphere.
#
library(MASS)
set.seed(17)
n.1 <- 100
n.2 <- 50
mu.1 <- c(0,-1,-1/4) * 2        # Center of first distribution
mu.2 <- c(1,1,-1/2) * 5         # Center of the second distribution
Sigma <- outer(1:3, 1:3, "==")  # Identity covariance matrix
xyz.1 <- t(mvrnorm(n.1, mu.1, Sigma)) # Each column is a point
xyz.2 <- t(mvrnorm(n.2, mu.2, Sigma))
xyz <- cbind(xyz.1, xyz.2)      # The Cartesian coordinates
rtf <- xyz.to.spherical(xyz)    # The spherical coordinates (also in columns)
#
# Compute the spherical mean and the azimuths relative to that mean.
#
mean.rtf <- xyz.to.spherical(rowMeans(xyz))
a <- SAS(rtf["theta",], rtf["phi",]-mean.rtf["phi",], mean.rtf["theta",])
#
# Plot the data and a QQ plot of the azimuths.
#
par(mfrow=c(1,2))
plot(c(-pi, pi), c(-1,1), type="n", 
     xlab="Phi", ylab="Cos(theta)", main="Sample points and their mean")
abline(h=0, col="Gray") # The Equator
abline(v=0, col="Gray") # The Prime Meridian
points(rtf["phi",], cos(rtf["theta", ]), col="#00000080")
points(mean.rtf["phi",], cos(mean.rtf["theta",]), bg="Red", pch=24, cex=1.25)

plot(c(0,1), c(0,2*pi), type="n",
     xlab="Quantile", ylab="Azimuth", main="Azimuthal QQ Plot")
abline(c(0, 2*pi), lty=3, lwd=2, col="Gray")
points(seq(0, 1, along.with=a), sort(a))

