如何找到该代码的良好近似值？

Question

我想回答这个问题但我觉得有点错误...

这是我的代码：

unitsize(1cm);
import geometry;

pair[] EllipseTangentNodes(ellipse el, pair M)
{
  pair[] op;
  op.push(intersectionpoints(tangents(el,(point) M)[0],el)[0]);
  op.push(intersectionpoints(tangents(el,(point) M)[1],el)[0]);
  return op;
}
arc ArcEllipse(ellipse el, pair M, pair N, bool direction=CCW)
{
  return arc(el, (point) M, (point) N, direction);
}

transform t=shift((-0.875,0));
ellipse e=t*ellipse((0,0), 1, 3, 25);
ellipse e1=shift((3,3))*ellipse((0,0), 2, 1, 35);
pair C=(0,4);

pair[] inter=EllipseTangentNodes(e,C);
pair A=inter[0],B=inter[1];
pair[] inter1=EllipseTangentNodes(e1,C);
pair A1=inter1[0],B1=inter1[1];

arc arcellipse=ArcEllipse(e,A,B);
arc arcellipse1=ArcEllipse(e1,B1,A1);

pair ahi(arc p, pair M, int n=1000)
{
  real k=arclength(p)/n;
  pair Q=relpoint(p,0);
  int j=0;
  while ( j <= n)
  {
    if (abs(M-Q) > abs(M-relpoint(p,j*k)))
      {
        Q=relpoint(p,j*k);
      }
    ++j;
  }
  return Q;
}
pair aho(arc p, pair M, int n=1000)
{
  real k=arclength(p)/n;
  real[] dis;
  for (int i=0; i < n; ++i) { dis[i]=abs(M-relpoint(p,i*k));}
  real[] d=sort(dis);
  return relpoint(p,search(d,min(dis))*k);
}
real Aho(arc p, pair M, int n=1000)
{
  real k=arclength(p)/n;
  real[] dis;
  for (int i=0; i < n; ++i) { dis[i]=abs(M-relpoint(p,i*k));}
  return min(dis);
}

pair G=relpoint(arcellipse,0.8);
pair F=relpoint(arcellipse,0.87);

pair T=ahi(ArcEllipse(e,G,F),C);
pair Tt=aho(ArcEllipse(e,G,F),C);

write(T); // (-1.32499014517229,2.41066630618231)
write((point) T @ arcellipse); // true
write(Tt); // (-1.21624693207659,2.30875337124993)
write((point) Tt @ arcellipse); // true

dot("G",G,dir(20));
dot("F",F,dir(55));
dot("T",T,dir(-115),red);

draw(C--T,darkgreen);
draw(e);
dot("C",C,dir(90),red);
dot("A",A);
dot("B",B,dir(145));
draw(C--B,blue);
draw(C--A,green);
shipout(bbox(2mm,invisible));

我觉得这部分好像失败了！

write(T); // (-1.32499014517229,2.41066630618231)
write((point) T @ arcellipse); // true
write(Tt); // (-1.21624693207659,2.30875337124993)
write((point) Tt @ arcellipse); // true

我怎样才能得到一个很好的近似值或者您可以得到更好的代码？

Answer 1

一种简单的方法，使用以下Asymptote graph()函数进行近似：

import solids;
size(200,0);

currentprojection=orthographic(camera=(70,32,-24));

guide3 gell=rotate(42,Y)*rotate(20,X)*path3(scale(5,3)*circle((0,0),1),XYplane);
triple A=(3,-7,2);

pair fdist(real t){return (t,abs(relpoint(gell,t)-A));}

int n=1024;
guide gfd= graph(fdist, 0,1,n);
real mind=min(gfd).y;
pair t=intersectionpoints(gfd,(-1,mind)--(2,mind))[0];
triple B=relpoint(gell,t.x);

draw(A--B,orange+thick());
draw(surface(gell),lightgray+opacity(.3),meshpen=nullpen,render(merge=true));
draw(gell,deepblue+thick());

revolution b=sphere(A,mind);
draw(surface(b),paleblue+opacity(.5));

dot(A--B,deepgreen);

label("$A$",A,Z+X);
label("$B$",B,Z+X);