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

[Russell Senior]

>>>>> "Alberto" == Alberto Vieira Ferreira Monteiro <albmont at centroin.com.br> writes:

Alberto> I guess this sample is required for some practical
Alberto> application, say a simulation for something done over the
Alberto> Earth. Then, I also guess that the sample does not have to be
Alberto> _absolutely_ exact, but a reasonable approximation can do
Alberto> it. And the ellipsoid is a rotation ellipsoid.

Alberto> This is my suggestion:

Alberto> (1) Divide the Ellipsoid by latitudes in _n_ horizontal
Alberto> slices in such a way that each slice can be considered
Alberto> "almost" spherical. Of course here lies the problem:
Alberto> depending on the purpose of the simulation, _n_ would be so
Alberto> big as to make it impractical

Alberto> (2) Compute the area of each slice (there are formulas for
Alberto> that, whose error is not very big - again, we rely on the
Alberto> purpose of the simulation)

Alberto> (3) Chose a random slice based on weight = area

Alberto> (4) Chose the random latitude by a uniform from the minimum
Alberto> to the maximum latitude (a much better approximation would
Alberto> give higher weight to the latitude closer to the equator)

Alberto> (5) lon = 2 pi runif(1) # :-)

Alberto> Now the question is: do you know the formulas to compute the
Alberto> area in (2)?  I know these formulas exist, I learned them in
Alberto> the last century, but I can't remember them and I don't know
Alberto> how to find them.

This is essentially the approach I (and my local helpers) have taken.
With garden-variety calculus, we have derived a function t(z) that
maps equal area over the oblate ellipsoid.  We select t from a uniform
distribution and then find the corresponding z dimension.  The
expression is quite hairy, but apparently analytically correct (I
haven't found any errors yet), if possibly unstable numerically.

The derivative with respect to z of area on the ellipsoid at a plane
of latitude intersecting at a height z above the equator, where the
ellipsoid has been scaled such that a is the polar radius and the
equatorial radius is 1, is: 

  dA/dz = 2 * pi * f

where

  f = sqrt((z^2 * (1 - a^2))/a^4 + 1)

You find the function t(z) such that dA/dt is constant.

Then you select from a uniform distribution and then find the value of
z that corresponds.


-- 
Russell Senior         ``I have nine fingers; you have ten.''
seniorr at aracnet.com