2020-10-24

记录求解spherical quadrilateral two sets of angles

matlab小tips
1文件名仅仅包含数字，字母，下划线。
2 圆括号的使用用来正常的数值计算，中括号常见用来表示矩阵。
3 粗心求平方与*和^符号搞混。
4. matlab 进行三角函数计算时往往都是弧度制。而且针对于tan的反函数求出的值只在1，4象限，需要利用if语句来确定角度定值。
利用z12, x12+y12,x21+y21

m=pi/180;
a12=75*m;a23=110*m;a34=60*m;a41=80*m;
th1=120*m;
%1.caculate th2a, th2b
z12=cos(a34);
x1=sin(a41)*sin(th1); y1=-(sin(a12)*cos(a41)+cos(a12)*sin(a41)*cos(th1));
z1=cos(a12)*cos(a41)-sin(a12)*sin(a41)*cos(th1);
A=y1*sin(a23);
B=sin(a23)*x1;
D=cos(a23)*z1-cos(a34);
%[th2a,th2b]=solve_trig(A,B,D);
X1=(-B+(B^2-D^2+A^2)^0.5)/(D-A);
X2=(-B-(B^2-D^2+A^2)^0.5)/(D-A);
th2a=2*atan(X1)*180/pi;
th2b=2*atan(X2)*180/pi;
th2a_r=2*atan(X1);
th2b_r=2*atan(X2);

%2.caculate th3a,th3b by using X12=s34*s3,Y12=s34*c3
s3a=(x1*cos(th2a_r)-y1*sin(th2a_r))/sin(a34);
c3a=(cos(a23)*(x1*sin(th2a_r)+y1*cos(th2a_r))-sin(a23)*z1)/sin(a34);
th3a_r=atan(s3a/c3a);
th3a=atan(s3a/c3a)*180/pi+180;
s3b=(x1*cos(th2b_r)-y1*sin(th2b_r))/sin(a34);
c3b=(cos(a23)*(x1*sin(th2b_r)+y1*cos(th2b_r))-sin(a23)*z1)/sin(a34);
th3b_r=atan(s3b/c3b);
th3b=atan(s3b/c3b)*180/pi-180;

%caculate th4a,th4b by using X21=s34*s4, Y21=s34*c4
x2_b=sin(a23)*sin(th2a_r);
y2_b=-(sin(a12)*cos(a23)+cos(a12)*sin(a23)*cos(th2a_r));
z2_b=cos(a12)*cos(a23)-sin(a12)*sin(a23)*cos(th2a_r);
s4a=(x2_b*cos(th1)-y2_b*sin(th1))/sin(a34);
c4a=(cos(a41)*(x2_b*sin(th1)+y2_b*cos(th1))-sin(a41)*z2_b)/sin(a34);
th4a=atan(s4a/c4a)*180/pi;

x2_b=sin(a23)*sin(th2b_r);
y2_b=-(sin(a12)*cos(a23)+cos(a12)*sin(a23)*cos(th2b_r));
z2_b=cos(a12)*cos(a23)-sin(a12)*sin(a23)*cos(th2b_r);
s4b=(x2_b*cos(th1)-y2_b*sin(th1))/sin(a34);
c4b=(cos(a41)*(x2_b*sin(th1)+y2_b*cos(th1))-sin(a41)*z2_b)/sin(a34);
th4b=atan(s4b/c4b)*180/pi-180;

display(th2a);display(th2b);display(th3a);display(th3b);display(th4a);display(th4b);