Найти кратчайший маршрут между двумя маяками

restart;
M1:=[5*cos(phi),5*sin(phi),0]: N1:=[5,0,0]: S:=[0,0,18]:
# Длина образующей полного конуса 
l:=sqrt(349):
# Отношения, в которых точки N и M делят образующую 
lambda1:=2/(l-2): lambda2:=1/(l-1):
# Функции расстояний в R^3 и в R2 
dist1:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2+(X[3]-Y[3])^2):
dist2:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2):
N:=simplify((N1+~lambda1*~S)/~(1+lambda1));
M:=simplify((M1+~lambda2*~S)/~(1+lambda2));
# Расстояние между M и N в R^3
d:=simplify(dist1(M,N));
# Значение угла phi в радианах
phi:=solve(d=9);
# Тот же угол phi в градусах
evalf(180*%/Pi);
# Угол между N1S и M1S на развёртке
alpha:=5*phi/l;
# Тот же угол в градусах
evalf(180*%/Pi);
# Координаты точек N и M на развёртке
N:=[l-2,0]: M:=[(l-1)*cos(alpha),(l-1)*sin(alpha)]:
# Точное (символьное) значение расстояния по геодезической конуса и его приближённое значени с 20 значащими цифрами 
simplify(dist2(N,M));
`Длина пути по геодезической усечённого конуса`=evalf[20](%);
 
# Длина N2S
l1:=sqrt(4^2+(18-3.6)^2);
 
# Движение через верхнее основание
 
# Углы s=N2O1P и t=N2O1Q в R^3
P:=[4*cos(s),4*sin(s),3.6]: Q:=[4*cos(t),4*sin(t),3.6]:
# Те же углы на развёртке
s1:=4*s/l1: t1:=4*t/l1:
# Точки P и Q на развёртке
P1:=[l1*cos(s1),l1*sin(s1)]: Q1:=[l1*cos(t1),l1*sin(t1)]:
# Целевая функция (общее расстояние)
F:=dist2(N,P1)+dist1(P,Q)+dist2(Q1,M):
# График F
plot3d(F, s=0..phi,t=s..phi,  axes=normal); 
L:=Optimization:-Minimize(F, {s>=0,s<=phi,t>=s,t<=phi},initialpoint={s=0.2,t=2});
`Длина пути через верхнее основание`=L[1];
# Углы s и t в градусах
s=eval(s*180/Pi,L[2])^o; t=eval(t*180/Pi,L[2])^o; 
 
# Движение через нижнее основание
 
# Углы s=N1O1P и t=N1O1Q в R^3
P:=[5*cos(s),5*sin(s),0]: Q:=[5*cos(t),5*sin(t),0]:
# Те же углы на развёртке
s1:=5*s/l: t1:=5*t/l:
# Точки P и Q на развёртке
P1:=[l*cos(s1),l*sin(s1)]: Q1:=[l*cos(t1),l*sin(t1)]:
# Целевая функция (общее расстояние)
F:=dist2(N,P1)+dist1(P,Q)+dist2(Q1,M):
# График F
plot3d(F, s=0..phi,t=s..phi,  axes=normal); 
Optimization:-Minimize(F, {s>=0,s<=phi,t>=s,t<=phi},initialpoint={s=0.1,t=2.4});
`Длина пути через нижнее основание`=%[1];

Длина геодезической только по конусу Lgc=
12.076998118505042854014738050457339440878321503798
Угол пересечения основания в градусах xº=
62.3948681324413550591376897686796979862332580073
Длины первого, второго участков по конусу и длина среднего участка по основанию
Lg1= 3.0250215836862859361474627548811158088984959378176
Lg2= 1.9322234860239374422901608246917857346366856705269
Lgm= 7.0892966334710640405113776854836550026257596954593
Общая длина геодезической
Lg= 12.046541703181287418949001265056556546160941303804

Длина геодезической только по конусу Lgc=
12.076998118505042854014738050457339440878321503798
Угол пересечения основания в градусах xº=
69.538690986113603376507943120961011300183363510552
Длины первого, второго участков по конусу и длина среднего участка по основанию
Lg1= 1.0715758486571688677955649015848165387910289317378
Lg2= 2.1528082016445345901870919921570426670392763862362
Lgm= 9.3690846549046458326413543175601694232712698846677
Общая длина геодезической
Lg= 12.593468705206349290624011211302028629101575202642

r1=4;
r2=5;
h=3.6;
dl1=1;
dl2=2;
d=9;
sgn=1;          /*Через какое основание считать геодезическую (верхнее "1", нижнее "-1")*/
sa=1/sqrt((h/(r2-r1))^2+1);
r=if(sgn==1,r1,r2);
l=r/sa;
l1=r2/sa-dl1;
l2=r2/sa-dl2;
gm=sa*asin(sqrt((d^2-(l2-l1)^2)/l1/l2)/2/sa);
print"Длина геодезической только по конусу Lgc=";
print Lgc=sqrt(l1^2+l2^2-2*l1*l2*cos(2*gm));
x=gm/2+pi/4;        /*Начальное значение угла пересечения x для метода Ньютона*/
z=10;           /*Число итераций для метода Ньютона*/
y=x-acos(l/l1*cos x)/2-acos(l/l2*cos x)/2+sgn*(gm-sa*x);
py=1-sin x/sqrt((l1/l)^2-(cos x)^2)/2-sin x/sqrt((l2/l)^2-(cos x)^2)/2-sgn*sa;
x=x-y/py;
if(z--,gotor-3,0);  /*Конец итераций*/
print"Угол пересечения основания в градусах xº=";
print x*180/pi;
print"Длины первого, второго участков по конусу и длина среднего участка по основанию";
print"Lg1=",Lg1=sgn*(sqrt(l1^2-(l*cos x)^2)-l*sin x);
print"Lg2=",Lg2=sgn*(sqrt(l2^2-(l*cos x)^2)-l*sin x);
print"Lgm=",Lgm=2*sin x*r;
print"Общая длина геодезической через основание";
print"Lg=",Lg=Lg1+Lgm+Lg2;