Касательные к сложной поверхности

Пусть дан массив координат точек в пространстве, который задает сложную поверхность (например, геоид), а так же координаты точки в пространстве, лежащей за пределами этой поверхности. Необходимо найти кривую, образованную касательными прямыми, проведенными из данной точки к поверхности. Поверхность я задаю в Octave следующим образом:

clc; clear all; close all;
E = wgs84Ellipsoid('km');
a=E.SemimajorAxis;
b=E.SemiminorAxis;
e=E.Eccentricity;

phi = -pi/2:pi/500:pi/2;
lambda = 0:2*pi/250:2*pi;
N = zeros(length(phi),length(lambda));
h=waitbar(0, 'Calculating of geoid''s heights + coords...');
count=0;
for p = 1:length(phi)
  for l = 1:length(lambda)
    count=count+1;
    N(p,l) = geoid_height(phi(p), lambda(l), 'egm2008-1')/1000;
    n=a^2/sqrt(a^2*cos(phi(p))^2+b^2*sin(phi(p))^2);
    xc(p,l) = (n+N(p,l))*cos(phi(p))*cos(lambda(l));
    yc(p,l) = (n+N(p,l))*cos(phi(p))*sin(lambda(l));
    zc(p,l) = (b^2/a^2*n+N(p,l))*sin(phi(p)); 
    waitbar(count/(length(phi)*length(lambda)))
  end
end
close(h)
figure; hold on; grid on;
surf(xc,yc,zc,'facealpha',0.5,'edgecolor','none')
axis equal

Функция geoid_height - аналог матлабовской функции geoidheight (в Octave такой нет). Исходник кода брала тут: https://github.com/chutsu/proto/blob/master/octave/third_party/geographiclib_toolbox_1.49/geoid_height.m Она задает высоты геоида над эллипсоидом-референсом.

Координаты точки можно задать так:

Phi = deg2rad(150);
Lambda = deg2rad(60);
h = 700;
H = sqrt(a^2*cos(Phi)^2+b^2*sin(Phi)^2)+h;


xcam = H*cos(Phi)*cos(Lambda);
ycam = H*cos(Phi)*sin(Lambda);
zcam = H*sin(Phi);
plot3(xcam,ycam,zcam,'r*','linewidth',3);

Очевидно, что здесь нужен какой-то численный метод, однако в голову совсем ничего не приходит.