[R] random uniform sample of points on an ellipsoid (e.g. WG

[R Heberto Ghezzo, Dr]

OK, I stand corrected, my previous post generated a point on an ellipsoid as a proyection of an uniform random  point in a sphere.
now, to generate a random point in a sphere I generate an uniform on one axis and an uniform angle in the second, since the sphere has the property that the perifery of all slices of equal thickness parallel to an equator have equal surface.
So, if I devide an oblate ellipsoid, like a geoid, in slices parallel to the equator, I can easily compute the surface area of the slices  
  lower <- 0
  upper <- pi/2
  a <- 1
  b <- 1.001
  e <- sqrt(1-a^2/b^2)
#
#   area revolution
#
  fun3 <- function(t,az,bz){
    2 * pi * bz * sin(t) * sqrt(az^2 * sin(t)^2 + bz^2 * cos(t)^2)
  }
  are3 <- integrate(fun3,lower,upper,az=a, bz=b)
  are3
  s <- 2 * pi *(a*a + b*b*atanh(e)/e)
  print(s/2)
so the integration works ok on the semi-ellipsoid, Now I devide it in 100 slices
#
  v <- rep(0,100)
  for(i in 1:100){
    l <- (i-1)*pi/200
    u <- i*pi/200
    v[i] <- integrate(fun3,lower=l, upper=u, az=a, bz = b)$value
  }
  print(sum(v))
#[1] 6.291565
  s/2
#[1] 6.29157
  w <- cumsum(v)
  w <- w/w[100]
#
now I choose a slice at random proportional to its area and a point in the latitude axis
the longitude is random 0-2pi.
#
 x <- runif(1)
 xx <- max(which(w < x))
 lat <- (xx+runif(1)) * pi / 200
 lon <- 2 * runif(1) * pi
#
 cat(lat*180/pi,lon*180/pi)
#
This point should be 'approximately' at random on the surface of the ellipsoid (revolution oblate). The error in the distribution of the point versus perfectly random is 'only' in the determination of latitude where I interpolate within a slice with a, straight line ' xx + runif(1) ' instead of following the curve. Taking 1000 slices instead of 100 obviously reduces this error to almost nothing, but not zero.
Am I correct this time?
Heberto Ghezzo
McGill University
Montreal - Canada