求助：matlab里球面坐标系下的非圆齿轮渐开线命令

爱无止境

我想根据一篇文献里的关系式重现齿轮圆锥曲线，前面公式及命令如下：
# I+ O2 P/ O. J& k+ Q+ ]fi1=[0:0.1:2*pi];
c=0.16;
1 K. }- s+ b. x& U0 |! ER=43.6;
: G& z' m" Y4 x8 [9 O. Ofi2=0.75*fi1+(c/3)*0.75*cos(3*fi1+3*pi/2);      %0.75=z1/z2=3/4,我直接换成0.75代入关系式了。
: ~8 N! O( u7 [$ `6 qdelta1=atan(0.75*(1-c*sin(3*fi1+3*pi/2)));
delta2=acot(0.75*(1-c*sin(3*fi1+3*pi/2)));
) x0 e0 i4 s) @$ f- _$ d
球面方程如下：
0 W+ N- N1 a2 sx1=R*(sin(delta1).*cos(fi1));
y1=R*sin(delta1).*sin(fi1);
7 v3 r6 F0 _+ v3 U# |7 \z1=R*cos(delta1);
- r: E8 U5 y0 J# E' ur1=x1*i+y1*j+z1*k     %这里的i,j,k前的系数是节曲线上的点在三个坐标上的坐标。写到程序里面的不是这样的，下面也是一样。这样写出来大家理解直观一些。
7 S7 ~+ w1 N. N. ?5 D7 Jx2=R*(sin(delta2).*cos(fi2));
: X* {9 k8 d+ L7 m# ^# z7 l# yy2=R*sin(delta2).*sin(fi2);
) n1 ^4 E, w- h+ Dz2=R*cos(delta2);
  x6 h, Z4 r( xr2=x2*i+y2*j+z2*k
$ l" h5 N2 P7 e% @  F
& q* ^( |! ?) j, N' v5 Kr1和r2的值验算出来的都是43.6，没有问题。但是我是不知道用什么样的命令得到下面的球面图形。由于本人对matlab不很了解，到这里用命令老是出错，各位英雄女侠帮帮忙。感激不尽！！！