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
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
