Orthogonal projections on ellipsoids in TikZ

The math problem and the algorithmic approach

As @Thruston suggests math is needed to solve the problem. Anyway this leads to a not trivial (quartic) equation which is hard to solve in the analytic way (let have a look at either similar question or Point-to-ellipse and point-to-ellipsoid distance equation analysis). So the idea is to solve that an equation numerically. On https://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/ I found a geometric and stable algorithm which find the point (orthogonal projection) on the ellipse by minimizing the distance from the original point.

The algorithm

The following steps and the image will suggest the idea.

1. Connect O and P in oder to get A_start (this allows to run the algorithm "on the right side" of the ellipse).
2. Draw a circle (blue) and get the midpoint of the two intersections with the blue circle and the ellipse.
3. Use the midpoint to draw a new smaller circle (purple) and iterate the process (i.e. red, orange, pink,...)

The code

The code needs the packages tikz and tkz-euclide and in particular usetikzlibrary{intersections} for the intersection points. I use tkz-euclide because I feel good with the commands. Anyway you can get the same result in pure tikz.

begin{tikzpicture}

% INITIAL DATA %
% the arbitrary point P
tkzDefPoint(3,2){P}
% the center of the ellipse
tkzDefPoint(0,0){O}
% use rotate=angle to set the desired orientation
path[draw,name path=theellipse,rotate=20] (O) ellipse (2cm and 1cm);
tkzLabelPoints[above right](P)
tkzLabelPoints[below left](O)

% STARTING POINT OF ALGORITHM %
path[name path=OP] (O)--(P);
path[name intersections={of=OP and theellipse,by={Aone}}];
% comment/erase if need next three code lines
tkzLabelPoint[above left](Aone){$A_{textrm{start}}$}
tkzDrawCircle[help lines](P,Aone)
tkzDrawPoints(Aone)

% ALGORITHM TO FIND THE ORTHOGONAL PROJECTION %
% set up a different number of steps if needed
% (algorithm converges relatively fast)
foreach i in {1,...,3}
{
% define a circle with center P through Aone
% (Astart for the first step)
tkzDefCircle[radius](P,Aone)
tkzGetLength{dPAone}
path[name path=circle] (P) circle (dPAone pt);

% find intersections of circle with ellipse (Aone, Atwo)
path[name intersections={of=circle and theellipse,by={Atwo,Aone}}];

% find a "proper" midpoint of Aone -- Atwo on the ellipse
tkzDefMidPoint(Aone,Atwo)tkzGetPoint{Aone}
path[name path=PAone] (P)--(Aone);
path[name intersections={of=PAone and theellipse,by={Aone}}];
}


% GET AND PRINT OUT THE DISTANCE
tkzDrawPoints(O,P,Aone)
tkzDrawSegment[red](P,Aone)
end{tikzpicture}

Just for comparison, you can do this very simply in Metapost using the solve macro and a suitable helper function.

documentclass[border=5mm]{standalone}
usepackage{luamplib}
begin{document}
mplibtextextlabel{enable}
begin{mplibcode}
beginfig(1);

    path e; pair p; numeric k;
    e = fullcircle xscaled 233 yscaled 144 rotated 10;
    p = 160 dir 142;

    vardef acute(expr t) =
        direction t of e dotprod (p - point t of e) > 0
    enddef;

    k = solve acute(0, 4);

    drawarrow p -- point k of e withcolor red;
    draw e; 
    dotlabel.top(btex $p$ etex, p);

endfig;
end{mplibcode}
end{document}

This is wrapped up in luamplib so you can compile it with lualatex.

Notes

• solve is explained on pages 176-177 of the Metafont book.

• The idea is that you define macro foo such that foo(x) is either true or false. Then you call solve foo(a, b) where a and b are values such that foo(a) is true and foo(b) is false. solve uses a binary search to find the edge value between a and b.

• In this case I have defined a macro called acute that uses the dotprod operator to tell us whether the tangent at point t of the ellipse makes an acute angle with the line from p to point t of the ellipse.

• solve finds the point at which the angle is no longer acute, which is therefore the point at which the line to p is orthogonal to the tangent, and is therefore the closest to p.

• Some skill and judgment is required to pick the correct initial values for different positions of p.

As you can see my explanation is rather longer than the code...

I suggest TikZ + gradient descent

documentclass[tikz]{standalone}
usepackage{tikz,tkz-euclide}

begin{document}

newcommand{boundellipse}[3]% center, xdim, ydim
    {(#1) ellipse (#2 and #3)}

makeatletter
xdefsx{-0.875} % shift x
xdefsy{0} % shift y
xdefra{1} % radius a
xdefrb{3} % radius b
xdefro{25} % rotation
pgfpointxy{0}{4}
xdefPx{thepgf@x}xdefPy{thepgf@y}

% let ang ("angle") be a free variable and run gradient descent
defang{234} % choose your favorite initial value
foreachiterationcounter in{1,...,20}{
    begin{tikzpicture}
        draw(-5,-3)rectangle(1,5);
        draw[shift={(-0.875,0)},rotate=25] boundellipse{0,0}{1}{3};
        node at (0,4)[circle,fill,inner sep=1.5pt]{};
        % evaluate Ellipse(ang)
        pgfpointxy{sx + rb*cos(ang)*sin(ro) + ra*sin(ang)*cos(ro)}
                    {sy - rb*cos(ang)*cos(ro) + ra*sin(ang)*sin(ro)}
        xdefQx{thepgf@x}xdefQy{thepgf@y}
        draw(Qx,Qy)circle(.1);
        % evaluate diff vector to target point
        xdefDx{thedimexprPx-Qx}
        xdefDy{thedimexprPy-Qy}
        draw[red,->](Qx,Qy)--+(Dx,Dy);
        % evaluate tangent line = d Ellipse(ang) / dang 
        pgfpointxy{- rb*sin(ang)*sin(ro) + ra*cos(ang)*cos(ro)}
                    {+ rb*sin(ang)*cos(ro) + ra*cos(ang)*sin(ro)}
        xdefTx{thepgf@x}
        xdefTy{thepgf@y}
        draw[blue,->](Qx,Qy)--+(Tx,Ty);
        % inner product
        pgfmathsetmacroInn{Dx*Tx + Dy*Ty}
        % rescale inner product
        pgfmathsetmacroinn{Inn / sqrt(Tx*Tx+Ty*Ty)}
        message{^^J thinbold: inn ^^J}
        % update angle
        pgfmathsetmacroang{ang + inn/10} % /10 is the step length
        xdefang{ang}
    end{tikzpicture}
}

end{document}