请教精通“亮点模型”的高手看一下程序
 

小弟是初学者，遇到一个问题如下：
现有脉冲宽度是5ms，频率是30kHz的正弦信号入射到椭球体上，照射角phi（椭球轴线与入射方向的夹角）取0，椭球亮点幅度反射因子A=(a^2*sin(theta)^2+b^2*cos(theta)^2)/(2*a),其中a,b为椭球体长，短半轴的长度，tan(theta)=b/a*tan(phi),绘制回波信号pb的幅值随时间变化图形？
我编的程序是：a=70;
b=8;
r=1000;
phi=0;
c=1500;
f=30000;
fs=30000*10;
k=2*pi*f/c;
t=0:1/fs:0.005;
d=b/a*tan(phi);
theta=atan(d);
p0=sin(2*pi*f*t);
A=(a^2*sin(theta)^2+b^2*cos(theta)^2)/(2*a) ;
pb=(exp(i*k*r)*conj(exp(i*k*r)))^.5/r*p0*A
plot(t,pb)
但结果图形好像不对，请高手指正错误，等待中！