syms px py pz theta1 theta2 a1 a2 A r t phi;
%t=0:0.0001:pi/40;
sd=0.2-(t-sin(40*t)/40)/4;
px=0.3,py=sd,pz=0.1,a1=0.3,a2=0.25;
A=(px.^2+py.^2+a1.^2-a2.^2)/(2.*a1.*sqrt(px.^2+py.^2));
phi=atan(px./(py+eps));
r=sqrt(px.^2+py.^2);
%theta1=atan(A/sqrt(1-A^2))-phi;
%theta2=(atan(r.*cos(theta1+phi)./(r.*sin(theta1+phi)-a1)))./pi.*180;
theta1=(atan(A/sqrt(1-A^2))-phi)/pi*180;
df=diff(theta1,1);
t=0:0.0001:pi/40;
plot(t,df)

画不了图，求解救