function []=testify(a1,b1,a,b,r)
%----這里的a1,b1是第一個psd中心點Op1在軸瓦坐標系XOY下的坐標----
%----a,b表示激光點在軸瓦坐標系XOY下的坐標，r是軸瓦的半徑，L是psd中心點到軸瓦邊緣的距離----
%----alphi指的是桿與鉛垂線的夾角，范圍是[-90,90]----
%------alphi的正負確定方法：測量時桿向左邊傾斜則為正值，右邊傾斜則為負值-------
%----其中：ki指的是直線mi的斜率，（xi,yi）是直線mi與圓的交點的坐標，Opi是psd的中心點，Ai是直線mi與圓的交點----

%------對第一個點給出a1,b1,alphi_1的值，求出L--------
clc;
alphi_1=input('alphi_1=')  %----給出alphi_1的值----
if (alphi_1~=0)&&(alphi_1~=90)&&(alphi_1~=-90)%----除開0度，90度，-90度時的計算流程----
  if alphi_1>0               %----根據給出的alphi_1判斷經過Op1點的直線m1的斜率----
      k1=tan(pi/2-alphi_1*pi/180);
  else
      k1=tan(pi/2+abs(alphi_1*pi/180));
  end
  %f1=sym('y-b1=k1*(x-a1)');%----第一條直線的方程----由于numeric函數不可用，所以直接把交點的計算坐標給出----
  %f=sym('x^2+y^2=r^2');       %----圓的方程----
  %[x,y]=solve(f,f1)         %----求出第一條直線與圓的交點,交點有兩個，需要相應的判斷----
  %pause
  x1=[ 1/2/(1+k1^2)*(2*k1^2*a1-2*b1*k1+2*(r^2-b1^2-k1^2*a1^2+2*b1*k1*a1+k1^2*r^2)^(1/2)),1/2/(1+k1^2)*(2*k1^2*a1-2*b1*k1-2*(r^2-b1^2-k1^2*a1^2+2*b1*k1*a1+k1^2*r^2)^(1/2))];
  y1=[ b1+1/2*k1/(1+k1^2)*(2*k1^2*a1-2*b1*k1+2*(r^2-b1^2-k1^2*a1^2+2*b1*k1*a1+k1^2*r^2)^(1/2))-k1*a1,b1+1/2*k1/(1+k1^2)*(2*k1^2*a1-2*b1*k1-2*(r^2-b1^2-k1^2*a1^2+2*b1*k1*a1+k1^2*r^2)^(1/2))-k1*a1];
                                    %----直線m1與圓的交點----
  double(x1);
  double(y1);
if y1(1)<y1(2)                  %----因為只取圓的下半部分，所以只取縱坐標值小的那個點計算----
    L=sqrt((x1(1)-a1)^2+(y1(1)-b1)^2);
else
    L=sqrt((x1(2)-a1)^2+(y1(2)-b1)^2);
end
end
switch alphi_1               %----三個特殊的角度時的L值計算----
    case 0                   %----夾角0度時的L值----
    L=b1+sqrt(r^2-a1^2);
    case 90                  %----夾角90度時的L值----
    L=a1+sqrt(r^2-b1^2);
    case -90                 %----夾角-90度時的L值----
    L=sqrt(r^2-b1^2)-a1;
end
L

%------對第二個點給出alphi_2,求出相應的a2,b2------
alphi_2=input('alphi_2=')    %----給出alphi_2的值----
if (alphi_2~=0)&&(alphi_2~=90)&&(alphi_2~=-90)   %----除開0度，90度，-90度時的計算流程----
  if alphi_2>0
      k2=tan(pi/2-alphi_2*pi/180);
  else
    k2=tan(pi/2+abs(alphi_2*pi/180));
  end
  %f2=sym('y-b1=k2*(x-a1)');             %----經過Op2點的直線方程m2----
  %[x,y]=solve(f,f2)                     %----算出m2與圓的交點A2有兩個，需要相應的判斷----
  %pause
  x2=[1/2/(1+k2^2)*(2*k2^2*a1-2*b1*k2+2*(-b1^2-k2^2*a1^2+2*b1*k2*a1+r^2+k2^2*r^2)^(1/2)),1/2/(1+k2^2)*(2*k2^2*a1-2*b1*k2-2*(-b1^2-k2^2*a1^2+2*b1*k2*a1+r^2+k2^2*r^2)^(1/2))];
  y2=[b1+1/2*k2/(1+k2^2)*(2*k2^2*a1-2*b1*k2+2*(-b1^2-k2^2*a1^2+2*b1*k2*a1+r^2+k2^2*r^2)^(1/2))-k2*a1,b1+1/2*k2/(1+k2^2)*(2*k2^2*a1-2*b1*k2-2*(-b1^2-k2^2*a1^2+2*b1*k2*a1+r^2+k2^2*r^2)^(1/2))-k2*a1];
  double(x2);
  double(y2);
  if y2(1)<y2(2)                                    %----取縱坐標值小的那個點計算----
      Op1A2=sqrt((x2(1)-a1)^2+(y2(1)-b1)^2);         %----Op1點到A2點的距離----
  else
      Op1A2=sqrt((x2(2)-a1)^2+(y2(2)-b1)^2);
  end
  Op1Op2=abs(Op1A2-L);                         %----兩個psd中心的距離----
  if (y2(1)<y2(2))&&(Op1A2>L)     %----如果Op1A2大于L，說明桿往下移動以保證Op2A2=L,并取縱坐標值小的那個點計算----
      a2=a1-Op1Op2*(a1-x2(1))/Op1A2;          %----利用等比定理計算Op2點的坐標----
      b2=b1-Op1Op2*(b1-y2(1))/Op1A2;
  elseif (y2(1)>y2(2))&&(Op1A2>L)
      a2=a1-Op1Op2*(a1-x2(2))/Op1A2;
      b2=b1-Op1Op2*(b1-y2(2))/Op1A2;
  elseif (y2(1)<y2(2))&&(Op1A2<L) %----如果Op1A2小于L，說明桿往上移動以保證Op2A2=L,并取縱坐標值小的那個點計算----
      a2=(a1*L-Op1Op2*x2(1))/(L-Op1Op2);     %----利用等比定理計算Op2點的坐標----
      b2=(b1*L-Op1Op2*y2(1))/(L-Op1Op2);
  else 
      a2=(a1*L-Op1Op2*x2(2))/(L-Op1Op2);
      b2=(b1*L-Op1Op2*y2(2))/(L-Op1Op2);
  end
end

switch alphi_2                     %----三個特殊的角度的(a2,b2)坐標計算----
    case 0                         %----alphi_2 為0度----
    Op1A2=b1+sqrt(r^2-a1^2);
    Op1Op2=abs(Op1A2-L);
    a2=a1;                         %----橫坐標不變----
    y2(1)=-sqrt(r^2-a1^2);         %----求出交點的縱坐標----
    if Op1A2>L                     %----如果Op1A2大于L，說明桿往下移動以保證Op2A2=L----
        b2=b1-Op1Op2*(b1-y2(1))/Op1A2;
    else
        b2=(b1*L-Op1Op2*y2(1))/(L-Op1Op2); %----如果Op1A2小于L，說明桿往上移動以保證Op2A2=L----  
    end
    case 90                        %----alphi_2 為90度----
    Op1A2=a1+sqrt(r^2-b1^2);
    Op1Op2=abs(Op1A2-L);
    b2=b1;                         %----縱坐標不變----
    x2(1)=-sqrt(r^2-b1^2);         %----求出交點的橫坐標----
    if Op1A2>L                     %----原理同上----
       a2=a1-Op1Op2*(a1-x2(1))/Op1A2;
    else
       a2=(a1*L-Op1Op2*x2(1))/(L-Op1Op2);
    end
    case -90                       %----alphi_2 為-90度----
    Op1A2=sqrt(r^2-b1^2)-a1;
    Op1Op2=abs(Op1A2-L);
    b2=b1;                         %----縱坐標不變----
    x2(1)=sqrt(r^2-b1^2);          %----求出交點的橫坐標----
    if Op1A2>L                     %----原理同上----
       a2=a1-Op1Op2*(a1-x2(1))/Op1A2;
    else
       a2=(a1*L-Op1Op2*x2(1))/(L-Op1Op2);
    end
end

    a2=a-a2
    b2=b-b2


%------對第三個點給出alphi_3,求出相應的a3,b3------
alphi_3=input('alphi_3=')
if (alphi_3~=0)&&(alphi_3~=90)&&(alphi_3~=-90)
  if alphi_3>0
      k3=tan(pi/2-alphi_3*pi/180);
  else
      k3=tan(pi/2+abs(alphi_3*pi/180));
  end
  %f3=sym('y-b1=k3*(x-a1)');  %----經過Op3點的直線方程m3----
  %[x,y]=solve(f,f3);         %----算出m3與圓的交點A3有兩個，需要相應的判斷----
  x3=[1/2/(1+k3^2)*(2*k3^2*a1-2*b1*k3+2*(-b1^2-k3^2*a1^2+2*b1*k3*a1+r^2+k3^2*r^2)^(1/2)),1/2/(1+k3^2)*(2*k3^2*a1-2*b1*k3-2*(-b1^2-k3^2*a1^2+2*b1*k3*a1+r^2+k3^2*r^2)^(1/2))];
  y3=[b1+1/2*k3/(1+k3^2)*(2*k3^2*a1-2*b1*k3+2*(-b1^2-k3^2*a1^2+2*b1*k3*a1+r^2+k3^2*r^2)^(1/2))-k3*a1,b1+1/2*k3/(1+k3^2)*(2*k3^2*a1-2*b1*k3-2*(-b1^2-k3^2*a1^2+2*b1*k3*a1+r^2+k3^2*r^2)^(1/2))-k3*a1];
  double(x3);
  double(y3);
  if y3(1)<y3(2)                 %----取縱坐標值小的那個點計算----
      Op1A3=sqrt((x3(1)-a1)^2+(y3(1)-b1)^2);%----Op1點到A3點的距離----
  else
      Op1A3=sqrt((x3(2)-a1)^2+(y3(2)-b1)^2);
  end
  Op1Op3=abs(Op1A3-L);                    %----兩個psd中心的距離----
  if (y3(1)<y3(2))&&(Op1A3>L)     %----如果Op1A3大于L，說明桿往下移動以保證Op3A3=L,并取縱坐標值小的那個點計算----
      a3=a1-Op1Op3*(a1-x3(1))/Op1A3;      %----利用等比定理計算Op3點的坐標----
      b3=b1-Op1Op3*(b1-y3(1))/Op1A3;
  elseif (y3(1)>y3(2))&&(Op1A3>L)
      a3=a1-Op1Op3*(a1-x3(2))/Op1A3;
      b3=b1-Op1Op3*(b1-y3(2))/Op1A3;
  elseif (y3(1)<y3(2))&&(Op1A3<L) %----如果Op1A3小于L，說明桿往上移動以保證Op3A3=L,并取縱坐標值小的那個點計算----
      a3=(a1*L-Op1Op3*x3(1))/(L-Op1Op3);  %----利用等比定理計算Op3點的坐標----
      b3=(b1*L-Op1Op3*y3(1))/(L-Op1Op3);
  else 
      a3=(a1*L-Op1Op3*x3(2))/(L-Op1Op3);
      b3=(b1*L-Op1Op3*y3(2))/(L-Op1Op3);
  end
end
switch alphi_3
    case 0                        %----三個特殊的角度的坐標計算----
    Op1A3=b1+sqrt(r^2-a1^2);
    Op1Op3=abs(Op1A3-L);
    a3=a1;
    y3(1)=-sqrt(r^2-a1^2);
    if Op1A3>L
        b3=b1-Op1Op3*(b1-y3(1))/Op1A3;
    else
        b3=(b1*L-Op1Op3*y3(1))/(L-Op1Op3);
    end
    case 90
    Op1A3=a1+sqrt(r^2-b1^2);
    Op1Op3=abs(Op1A3-L);
    b3=b1;
    x3(1)=-sqrt(r^2-b1^2);
    if Op1A3>L
       a3=a1-Op1Op3*(a1-x3(1))/Op1A3;
    else
       a3=(a1*L-Op1Op3*x3(1))/(L-Op1Op3);
    end
    case -90
    Op1A3=sqrt(r^2-b1^2)-a1;
    Op1Op3=abs(Op1A3-L);
    b3=b1;
    x3(1)=sqrt(r^2-b1^2);
    if Op1A3>L
       a3=a1-Op1Op3*(a1-x3(1))/Op1A3;
    else
       a3=(a1*L-Op1Op3*x3(1))/(L-Op1Op3);
    end
end

    a3=a-a3
    b3=b-b3


%------對第四個點給出alphi_4,求出相應的a4,b4------
alphi_4=input('alphi_4=')
if (alphi_4~=0)&&(alphi_4~=90)&&(alphi_4~=-90)
  if alphi_4>0
      k4=tan(pi/2-alphi_4*pi/180);
  else
      k4=tan(pi/2+abs(alphi_4*pi/180));
  end
  %f4=sym('y-b1=k4*(x-a1)');  %----經過Op4點的直線方程m4----
  %[x,y]=solve(f,f4);         %----算出m4與圓的交點A4有兩個，需要相應的判斷----
  x4=[1/2/(1+k4^2)*(2*k4^2*a1-2*b1*k4+2*(-b1^2-k4^2*a1^2+2*b1*k4*a1+r^2+k4^2*r^2)^(1/2)),1/2/(1+k4^2)*(2*k4^2*a1-2*b1*k4-2*(-b1^2-k4^2*a1^2+2*b1*k4*a1+r^2+k4^2*r^2)^(1/2))];
  y4=[b1+1/2*k4/(1+k4^2)*(2*k4^2*a1-2*b1*k4+2*(-b1^2-k4^2*a1^2+2*b1*k4*a1+r^2+k4^2*r^2)^(1/2))-k4*a1,b1+1/2*k4/(1+k4^2)*(2*k4^2*a1-2*b1*k4-2*(-b1^2-k4^2*a1^2+2*b1*k4*a1+r^2+k4^2*r^2)^(1/2))-k4*a1];
  double(x4);
  double(y4);
  if y4(1)<y4(2)                 %----取縱坐標值小的那個點計算----
      Op1A4=sqrt((x4(1)-a1)^2+(y4(1)-b1)^2);%----Op1點到A4點的距離----
  else
      Op1A4=sqrt((x4(2)-a1)^2+(y4(2)-b1)^2);
  end
  Op1Op4=abs(Op1A4-L);                    %----兩個psd中心的距離----
  if (y4(1)<y4(2))&&(Op1A4>L)     %----如果Op1A4大于L，說明桿往下移動以保證Op4A4=L,并取縱坐標值小的那個點計算----
      a4=a1-Op1Op4*(a1-x4(1))/Op1A4;      %----利用等比定理計算Op4點的坐標----
      b4=b1-Op1Op4*(b1-y4(1))/Op1A4;
  elseif (y4(1)>y4(2))&&(Op1A4>L)
      a4=a1-Op1Op4*(a1-x4(2))/Op1A4;
      b4=b1-Op1Op4*(b1-y4(2))/Op1A4;
  elseif (y4(1)<y4(2))&&(Op1A4<L) %----如果Op1A4小于L，說明桿往上移動以保證Op4A4=L,并取縱坐標值小的那個點計算----
      a4=(a1*L-Op1Op4*x4(1))/(L-Op1Op4);  %----利用等比定理計算Op4點的坐標----
      b4=(b1*L-Op1Op4*y4(1))/(L-Op1Op4);
  else 
      a4=(a1*L-Op1Op4*x4(2))/(L-Op1Op4);
      b4=(b1*L-Op1Op4*y4(2))/(L-Op1Op4);
  end
end

switch alphi_4
    case 0                       %----三個特殊的角度的坐標計算----
    Op1A4=b1+sqrt(r^2-a1^2);
    Op1Op4=abs(Op1A4-L);
    a4=a1;
    y4(1)=-sqrt(r^2-a1^2);
    if Op1A4>L
        b4=b1-Op1Op4*(b1-y4(1))/Op1A4;
    else
        b4=(b1*L-Op1Op4*y4(1))/(L-Op1Op4);
    end
    case 90
    Op1A4=a1+sqrt(r^2-b1^2);
    Op1Op4=abs(Op1A4-L);
    b4=b1;
    x4(1)=-sqrt(r^2-b1^2);
    if Op1A4>L
       a4=a1-Op1Op4*(a1-x4(1))/Op1A4;
    else
       a4=(a1*L-Op1Op4*x4(1))/(L-Op1Op4);
    end
    case -90
    Op1A4=sqrt(r^2-b1^2)-a1;
    Op1Op4=abs(Op1A4-L);
    b4=b1;
    x4(1)=sqrt(r^2-b1^2);
    if Op1A4>L
       a4=a1-Op1Op4*(a1-x4(1))/Op1A4;
    else
       a4=(a1*L-Op1Op4*x4(1))/(L-Op1Op4);
    end
end
 
    a4=a-a4
    b4=b-b4


%----利用原來的方程驗證求出的a2,b2,a3,b3,a4,b4能否滿足方程-----
r1=sqrt((L*sin(alphi_1*pi/180)-a1)^2+(L*cos(alphi_1*pi/180)-b1)^2)
r2=sqrt((L*sin(alphi_2*pi/180)+a2-a)^2+(L*cos(alphi_2*pi/180)+b2-b)^2)
r3=sqrt((L*sin(alphi_3*pi/180)+a3-a)^2+(L*cos(alphi_3*pi/180)+b3-b)^2)
r4=sqrt((L*sin(alphi_4*pi/180)+a4-a)^2+(L*cos(alphi_4*pi/180)+b4-b)^2)