points.on.unit.sphere<-function(n)
{ # attempts to regularly distribute approximately n points on 
  # the surface of a unit sphere non-iteratively
  # by laying them out along a spiral with a fixed (angular) pitch, c,
  # x,y,z are the cartesian coordinates of the points,
  # theta is their longitude, phi their lattitude (in radians)
  # by Mike Lonergan, mel@mcs.st-and.ac.uk

  c<-sqrt(n*pi)/2
  theta<-c(0,2*pi)
  for(i in 3:floor(n/2))
      theta[i]<-theta[i-1]+pi/(c*cos(theta[i-1]/(2*c)-pi/2))
  #   theta[i]<-sqrt(2*c+theta[i-1]^2)
  if (2*floor(n/2)==n)
     theta<-c(theta,2*pi*c-rev(theta))
  else
     theta<-c(theta, pi*c, 2*pi*c-rev(theta))
     


  pts<-data.frame(theta=theta)

  pts$phi<-theta/(2*c)-pi/2
  pts$x<-cos(theta)*cos(pts$phi)
  pts$y<-sin(theta)*cos(pts$phi)
  pts$z<-sin(pts$phi)

  pts
}


nearest<-function(data)
{ # takes a dataframe with columns x, y, z and returns the (straightline) nearest
  # neighbour distances between the points in its rows
  # inefficient, but adequate for checking points.on.unit.sphere.

 res<-NA
 for (i in 1:dim(data)[1])
    res[i]<-sqrt(min((data$x[-i]-data$x[i])^2 + (data$y[-i]-data$y[i])^2 +
          (data$z[-i]-data$z[i])^2))
 
res
}

points.on.unit.sphere(1000)->pous
nearest(pous)->npous

par(mfrow=c(2,2))
hist(npous,main="nearest neighbour distances")
plot(pous$x+sign(pous$z),pous$y,main="plan view")
plot(pous$x+sign(pous$y),pous$z,main="side.view")
plot(npous,ylab="nearest neighbour distances",xlab="theta")
length(which(npous<0.06))