基于[(R-h1)sina1-x1]^2+[(R-h1)cosa1+z1]^2=r^2 [(R-h2)sina2-x1]^2+[(R-h2)cosa2+z1]^2=r^2 用matlab绘制轴心轨迹

首先，我们需要将两个方程组合起来，消去变量 $z_1$，并且将其转化为 MATLAB 中的可计算形式。做如下变换：
$$
\begin{cases}
(R-h_1)^2\sin^2 a_1 + (R-h_1)^2\cos^2 a_1 + 2(R-h_1)\cos a_1 z_1 + z_1^2 = r^2 \\
(R-h_2)^2\sin^2 a_2 + (R-h_2)^2\cos^2 a_2 + 2(R-h_2)\cos a_2 z_1 + z_1^2 = r^2 \\
\end{cases}
$$

将两个方程相减，可得：
$$
\begin{aligned}
&(R-h_1)^2\sin^2 a_1 - (R-h_2)^2\sin^2 a_2 + (R-h_1)^2\cos^2 a_1 - (R-h_2)^2\cos^2 a_2 + \\
&2(R-h_1)\cos a_1 z_1 - 2(R-h_2)\cos a_2 z_1 = 0 \\
\end{aligned}
$$

整理一下，可得：
$$
z_1 = \frac{(R-h_2)^2\sin^2 a_2 - (R-h_1)^2\sin^2 a_1 - (R-h_2)^2\cos^2 a_2 + (R-h_1)^2\cos^2 a_1}{2(R-h_2)\cos a_2 - 2(R-h_1)\cos a_1}
$$

现在，我们可以使用 MATLAB 画出轴心轨迹。假设 $R=1, h_1=0.5, h_2=1.5, a_1=\pi/4, a_2=\pi/6, r=0.5$，则代码如下：

R = 1;
h1 = 0.5;
h2 = 1.5;
a1 = pi/4;
a2 = pi/6;
r = 0.5;

x1 = 0:0.01:2;
z1 = ((R-h2)^2*sin(a2)^2 - (R-h1)^2*sin(a1)^2 - (R-h2)^2*cos(a2)^2 + (R-h1)^2*cos(a1)^2) ./ (2*(R-h2)*cos(a2) - 2*(R-h1)*cos(a1));

y1 = sqrt(r^2 - (x1-(R-h1)*sin(a1)).^2) + (R-h1)*cos(a1) + z1;
y2 = -sqrt(r^2 - (x1-(R-h1)*sin(a1)).^2) + (R-h1)*cos(a1) + z1;
y3 = sqrt(r^2 - (x1-(R-h2)*sin(a2)).^2) + (R-h2)*cos(a2) + z1;
y4 = -sqrt(r^2 - (x1-(R-h2)*sin(a2)).^2) + (R-h2)*cos(a2) + z1;

plot3(x1, y1, z1);
hold on;
plot3(x1, y2, z1);
plot3(x1, y3, z1);
plot3(x1, y4, z1);
grid on;
xlabel('x');
ylabel('y');
zlabel('z');

运行代码，可得以下轴心轨迹的图像：