问点到椭圆弧最短距离的算法

EN

double x0,y0,rx,ry,a0,a1;   // elliptic arc center,semi-axises,start/end angles CW
void ellarc_closest_point(double &x_out,double &y_out,double x_in,double y_in)
    {
    int e,i;
    double ll,l,aa,a,da,x,y,b0,b1;
    while (a0>=a1) a0-=pi2;                 // just make sure a0<a1
    b0=a0; b1=a1; da=(b1-b0)/25.0;          // 25 sample points in first iteration
    ll=-1; aa=a0;                           // no best solution yet
    for (i=0;i<3;i++)                       // recursions more means more accurate result
        {
        // sample arc a=<b0,b1> with step da
        for (e=1,a=b0;e;a+=da)
            {
            if (a>=b1) { a=b1; e=0; }
            // elliptic arc sampled point
            x=x0+rx*cos(a);
            y=y0-ry*sin(a);                 // mine y axis is in reverse order therefore -
            // distance^2 to x_in,y_in
            x-=x_in; x*=x;
            y-=y_in; y*=y; l=x+y;
            // remember best solution
            if ((ll<0.0)||(ll>l)) { aa=a; ll=l; }
            }
        // use just area near found solution aa
        b0=aa-da; if (b0<a0) b0=a0;
        b1=aa+da; if (b1>a1) b1=a1;
        // 10 points per area stop if too small area already
        da=0.1*(b1-b0); if (da<1e-6) break;
        }
    x_out=x0+rx*cos(aa);
    y_out=y0-ry*sin(aa);    // mine y axis is in reverse order therefore -
    }