[R] points on a sphere

[Mike Lonergan]

Hi Richard,

Here's a non-iterative way of approximately laying out points on the surface
of a sphere. It basically lays out a spiral on the surface (like peeling an
apple) then drops the points along it. I don't know if it will be any actual
use to you, but I quite like it (and would be interested to know what breaks
it).

Cheers,

Mike.


----------------------------------------------------------------
Mike Lonergan

The Observatory, Buchanan Gardens 	ph: (+44) 1334 461803
University of St Andrews		fax:(+44) 1334 461800
Fife, KY16 9LZ, Scotland		email: mel at mcs.st-and.ac.uk
--------------------------------------------------------------------



##########

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

##########

     > -----Original Message-----
     > From: owner-r-help at stat.math.ethz.ch
     > [mailto:owner-r-help at stat.math.ethz.ch]On Behalf Of Richard Rowe
     > Sent: 26 October 2002 00:46
     > To: r-help at stat.math.ethz.ch
     > Subject: [R] points on a sphere
     >
     >
     > Not an R question directly, but has anyone got a method for
     > placing a
     > moderately large number of (near) equi-spaced points on a
     > sphere?  I have a
     > nasty feeling platonic solids are needed for exact solutions and I'm
     > thinking of samplings involving around 200 - 1000
     > regularly-spaced points,
     >
     > Thanks,
     >
     > Richard Rowe
     >
     > Richard Rowe
     > Senior Lecturer
     > Department of Zoology and Tropical Ecology, James Cook University
     > Townsville, Queensland 4811, Australia
     > fax (61)7 47 25 1570
     > phone (61)7 47 81 4851
     > e-mail: Richard.Rowe at jcu.edu.au
     > http://www.jcu.edu.au/school/tbiol/zoology/homepage.html
     >
     > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
     > .-.-.-.-.-.-.-.-.-
     > r-help mailing list -- Read
     > http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
     > Send "info", "help", or "[un]subscribe"
     > (in the "body", not the subject !)  To:
     > r-help-request at stat.math.ethz.ch
     > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
     > ._._._._._._._._._
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: points on unit sphere.txt
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20021028/c5b5a4b0/pointsonunitsphere.txt