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))