lt5m.m

自己编的matlab厂程序。这是将局部坐标系的坐标转化为球面坐标系的坐标

   Am4=[-0.1494;0.0135;0.9710];%在此输入1号点的坐标，其顺序是东坐标，北坐标，高程
   Am6=[0.1232;-0.1406;0.9932];%在此输入2号点的坐标，其顺序同上
   Am5=[0.1227;0.1730;1.0063];%在此输入3号点的坐标，其顺序同上
   T=[0,1,0;-1,0,0;0,0,1];
   rb1=T*Am6.*100;
   rb2=T*Am5.*100;
   rb3=T*Am4.*100;
   Am1=rb1;
   Am2=rb2;
   Am3=rb3;
   a=(rb2(2)-rb1(2))*(rb3(3)-rb1(3))-(rb2(3)-rb1(3))*(rb3(2)-rb1(2));
   b=(rb2(3)-rb1(3))*(rb3(1)-rb1(1))-(rb2(1)-rb1(1))*(rb3(3)-rb1(3));
   c=(rb2(1)-rb1(1))*(rb3(2)-rb1(2))-(rb2(2)-rb1(2))*(rb3(1)-rb1(1));
   d=rb1(3)*(rb2(2)-rb1(2))*(rb3(1)-rb1(1))+rb1(2)*(rb2(1)-rb1(1))*(rb3(3)-rb1(3))+rb1(1)*(rb2(3)-rb1(3))*(rb3(2)-rb1(2))-rb1(1)*(rb2(2)-rb1(2))*(rb3(3)-rb1(3))-rb1(2)*(rb2(3)-rb1(3))*(rb3(1)-rb1(1))-rb1(3)*(rb2(1)-rb1(1))*(rb3(2)-rb1(2));
   d2=d-12.64*(sqrt(a*a+b*b+c*c));
  % [x0,y0,z0]=solve('a*x0+b*y0+c*z0+d2','a/(rb1(1)-x0)-b/(rb1(2)-y0)','c/(rb1(3)-z0)-a/(rb1(1)-x0)','x0','y0','z0')
 %[x0,y0,z0]=solve('-304.7862*x0+(34.7600)*y0+808.0698*z0+(-1.0778e+005)','-304.7862/(38.9243-x0)-(34.7600)/(-9.5513-y0)','808.0698/(135.0723-z0)-(-304.7862)/(38.9243-x0)')%,'x0','y0','z0');
   x01 =-(b*a*rb1(2)-b^2*rb1(1)-c^2*rb1(1)+c*a*rb1(3)+d2*a)/(a^2+b^2+c^2);
   y01 =(a^2*rb1(2)+rb1(2)*c^2-b*rb1(1)*a-b*c*rb1(3)-b*d2)/(a^2+b^2+c^2);
   z01 =(-c*rb1(1)*a-c*b*rb1(2)-c*d2+a^2*rb1(3)+rb1(3)*b^2)/(a^2+b^2+c^2);
   point1=[x01;y01;z01];
    x02 =-(b*a*rb2(2)-b^2*rb2(1)-c^2*rb2(1)+c*a*rb2(3)+d2*a)/(a^2+b^2+c^2);
   y02 =(a^2*rb2(2)+rb2(2)*c^2-b*rb2(1)*a-b*c*rb2(3)-b*d2)/(a^2+b^2+c^2);
   z02 =(-c*rb2(1)*a-c*b*rb2(2)-c*d2+a^2*rb2(3)+rb2(3)*b^2)/(a^2+b^2+c^2);
   point2=[x02;y02;z02];
    x03 =-(b*a*rb3(2)-b^2*rb3(1)-c^2*rb3(1)+c*a*rb3(3)+d2*a)/(a^2+b^2+c^2);
   y03 =(a^2*rb3(2)+rb3(2)*c^2-b*rb3(1)*a-b*c*rb3(3)-b*d2)/(a^2+b^2+c^2);
   z03 =(-c*rb3(1)*a-c*b*rb3(2)-c*d2+a^2*rb3(3)+rb3(3)*b^2)/(a^2+b^2+c^2);
   point3=[x03;y03;z03];
   t=1/3*(point1+point2+point3);
    ox1=(Am1(1)+Am2(1)+Am3(1))/3;
    oy1=(Am1(2)+Am2(2)+Am3(2))/3;
    oz1=(Am1(3)+Am2(3)+Am3(3))/3;
    ro1=[ox1;oy1;oz1];
    %this absulte coordinate vector
    Vx=[1;0;0];Vy=[0;1;0];Vz=[0;0;1];
    V=[Vx';Vy';Vz'];
    %method1
    %Vx1=[-Am1(1)+ox1;-Am1(2)+oy1;-Am1(3)+oz1];
    Vx1=Am1-ro1;
    ax=(Am2(1)-ro1(1));ay=(Am2(2)-ro1(2));az=(Am2(3)-ro1(3));
    bx=(Am3(1)-ro1(1));by=(Am3(2)-ro1(2));bz=(Am3(3)-ro1(3));
    Vz1=[ay*bz-az*by;az*bx-ax*bz;ax*by-ay*bx];
    ax=Vz1(1);ay=Vz1(2);az=Vz1(3);
    bx=Vx1(1);by=Vx1(2);bz=Vx1(3);
    Vy1=[ay*bz-az*by;az*bx-ax*bz;ax*by-ay*bx];
    V11=[Vx1';Vy1';Vz1'];
    A11=V*V11';
    A1=A11';
    Axx=Vx'*Vx1/(sqrt(Vx1'*Vx1));Ayx=Vy'*Vx1/(sqrt(Vx1'*Vx1));Azx=Vz'*Vx1/(sqrt(Vx1'*Vx1));
    Axy=Vx'*Vy1/(sqrt(Vy1'*Vy1));Ayy=Vy'*Vy1/(sqrt(Vy1'*Vy1));Azy=Vz'*Vy1/(sqrt(Vy1'*Vy1));
    Axz=Vx'*Vz1/(sqrt(Vz1'*Vz1));Ayz=Vy'*Vz1/(sqrt(Vz1'*Vz1));Azz=Vz'*Vz1/(sqrt(Vz1'*Vz1));
    A=[Axx,Axy,Axz;Ayx,Ayy,Ayz;Azx,Azy,Azz];
    B=[sqrt(3)/2,0.5,0;-0.5,sqrt(3)/2,0;0,0,1]*A;
    %t=point4-B*[15.588;9;0];
    %o=[0;0;0];
    SSO=t;              %得到中心点的坐标
    %  球角局坐标系标 
    ss1=[20.7846;0;3];ss2=[-10.3923;18;3];ss3=[-10.3923;-18;3];ss4=[1.1238;2.7816;20.7846];ss5=[-2.9708;-0.4175;20.7846];ss6=[1.8470;-2.3640;20.7846];
    %  球角体坐标系标 
    SS1=B*ss1+t;
    SS2=B*ss2+t;
    SS3=B*ss3+t;
    SS4=B*ss4+t;
    SS5=B*ss5+t;
    SS6=B*ss6+t;
    output=[SSO,SS1,SS2,SS3,SS4,SS5,SS6]