From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Vectors on ellipsoid surface
Date: Sun, 23 May 2010 03:50:04 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 67
Message-ID: <hta8lc$e4f$>
References: <ht5ldf$7tf$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1274586604 14479 (23 May 2010 03:50:04 GMT)
NNTP-Posting-Date: Sun, 23 May 2010 03:50:04 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:638415

"Els " <> wrote in message <ht5ldf$7tf$>...
> I plotted an ellipsoid [x,y,z]=ellipsoid   (standard in Matlab)
> Now I have a point C(x,y,z) which lays within this surface.
> I am looking for the projection of C on the surface of the ellipsoid,called C', where the normal vector perpendicular to the surface is in line with the vector between C and C'.
> Thanks in advance. 

  The following is an analytic solution of your problem.  It is not an iterative method but gets its solutions directly using matlab's 'roots' function.  It assumes an ellipsoid with center at the origin and axes aligned with coordinates axes, satifying the equation

 x^2/a^2 + y^2/b^2 + z^2/c^2 = 1

The point (p,q,r) can be either inside or outside this ellipsoid.  Points (x,y,z) on the ellipsoid are found such that (p,q,r) lies on the line through (x,y,z) orthogonal to the surface of the ellipsoid there.  In general there will be six solutions to this problem though frequently some of these will be complex-valued.  There are always at least two real solutions.

  It should be noted that if two of the three semi-axes a, b, and c, are equal, giving an oblate or prolate spheroid, then there are in general only five solutions, though again some of these may be complex too.  It should be noted that in this case if (p,q,r) lies right on the axis of revolution, infinite rings of solutions about this axis become possible.  In the code below in this case the algorithm still produces six solutions, but two of them will be equal - that is, equal to within round off error.

  The algorithm below computes the appropriate coefficients of a sextic polynomial in the variable t

 A*t^6 + B*t^5 + C*t^4 + D*t^3 + E*t^2 + F*t + G = 0	

and then computes a set of (x,y,z) points on the surface for each root t according to the equations

 x = a^2*p/(a^2+t)
 y = b^2*q/(b^2+t)
 z = c^2*r/(c^2+t)

The significance of the t variable is that it is the common ratio

 (p-x)/(x/p^2) = (q-y)/(y/q^2) = (r-z)/(z/r^2) = t

between components of the line between (p,q,r) and (x,y,z) and the components of a line through (x,y,z) orthogonal to the surface there.

  The code assumes that a, b, c and p, q, r, have been defined.

- - - - - - - - - - - - - - - - - -
 % The program

 % Compute some convenient intermediate values
 p2 = p^2; q2 = q^2; r2 = r^2;
 a2 = a^2; b2 = b^2; c2 = c^2;
 sbc = b2+c2; sca = c2+a2; sab = a2+b2;
 mbc = b2*c2; mca = c2*a2; mab = a2*b2;
 s1 = a2+b2+c2; s2 = mbc+mca+mab; s3 = a2*b2*c2;

 % Now for the polynomial's actual coefficients
 A = 1;
 B = 2*s1;
 C = s1^2+2*s2-p2*a2-q2*b2-r2*c2;
 D = 2*(s1*s2+s3-p2*a2*sbc-q2*b2*sca-r2*c2*sab);
 E = s2^2+2*s3*s1-p^2*a^2*(sbc^2+2*mbc) ...
 F = 2*s3*(s2-p2*sbc-q2*sca-r2*sab);
 G = s3*(s3-p2*mbc-q2*mca-r2*mab);

 % Now find the polynomial's six roots
 t = roots([A,B,C,D,E,F,G]);

 % Derive the corresponding sets of coordinates.
 x = a^2*p./(a^2+t);
 y = b^2*q./(b^2+t);
 z = c^2*r./(c^2+t);
 % The resulting three six-element column vectors give the solutions.
- - - - - - - - - - - - - - - - - -

  It should be pointed out that, depending on the position of the (p,q,r) point relative to the ellipsoid, there can be a comparatively large effect from round off errors.   For example, if the ellipsoid is nearly a sphere and the point (p,q,r) is nearly at the origin, all points on the surface come close to being solutions, and this would have the effect of magnifying any round off errors in the solutions actually arrived at.  This difficulty is inherent in the nature of the stated problem.

Roger Stafford