根据两点以及半径示圆心

function [xo1,yo1,xo2,yo2]=getmid(x0,y0,x1,y1, R)
    oa = (1 + (x1-x0)^2/(y1-y0)^2);
    ob = -(2 * x1 + (x1-x0)*(x1^2 + 2*y1*y0 - y0^2 - x0^2 -y1^2)/(y1-y0)^2);
    oc = x1^2 + 1/4 * (x1^2+2*y1*y0 - y0^2-x0^2-y1^2)^2/(y1-y0)^2 -R^2;
    xo1 = (-ob - sqrt(ob^2 - 4* oa *oc))/(2* oa);
    xo2 = (-ob + sqrt(ob^2 - 4* oa *oc))/(2* oa);
    yo1 = ((x1^2 +y1^2 - x0^2 -y0^2) - 2*(x1-x0)*xo1)/(2 *(y1-y0));
    yo2 = ((x1^2 +y1^2 - x0^2 -y0^2) - 2*(x1-x0)*xo2)/(2 *(y1-y0));