肉毒素维持时间:各位大侠，我这有一个matlab程序，需要转换成vb或vc程序，请大家帮忙看看！谢谢！！！

function ans=calOneZR(f,f0,zs,ps)
% 求单个零点的相应函数值
% f0 腔体谐振频率点
% zs,ps 原有的零点和极点
ans=prod(f*f-[ps f0].*[ps f0])-prod(f*f-zs.*zs);

function ans=calP(f0,z,p,f,m)
% function calP(f0,z,p,f,m)
% 得到i－1阶P的函数值，为求其零点使用
num=size(z,2);
N=size(f,2);
if (num==1)
ans=(f0*f0-f(N-1)^2)*(f0*f0-z*z)-m(N-1)^2*f0*f0;
return;
end
ans=(f0*f0-f(N-num)*f(N-num))*prod(f0*f0-z.*z)-m(N-num)^2*f0*f0*prod(f0*f0-p.*p);

function [za,pa]=findNewZP(x,y,z)
% 由i阶输入阻抗的零点和极点得到i+1阶的零点和极点
% 新的零点和极点的数量均减少1
% z,p 新的零点和极点
% x,y,z i阶零点、极点和i腔谐振频率

za=y;
num=size(za,2);
pa=zeros(1,num-1);
for i=1:num-1
pa(i)=fzero(@calOneZR,[za(i) za(i+1)],[],z,x,y);
end

unction ans=findZinZeros(z,p,f,m)
% function ans=findZinZeros(p,f,m)
% 得到i－1阶输入阻抗的零点
% z,p i阶的零点和极点
% f,m 谐振频率和耦合强度

num=size(z,2); % i阶零点数量

delta=0.01;
while (calP(z(1),z,p,f,m)*calP(z(1)*(1-delta),z,p,f,m)>0)
delta=delta*2;
end
ans(1)=fzero('calP',[z(1)*(1-delta) z(1)],[],z,p,f,m);

for i=1:num-1
ans(i+1)=fzero('calP',[z(i) z(i+1)],[],z,p,f,m);
end

delta=0.01;
while (calP(z(num),z,p,f,m)*calP(z(num)*(1+delta),z,p,f,m)>0)
delta=delta*2;
end
ans(num+1)=fzero('calP',[z(num) z(num)*(1+delta)],[],z,p,f,m);

function ans=findZinZeros(z,p,f,m)
% function ans=findZinZeros(p,f,m)
% 得到i－1阶输入阻抗的零点
% z,p i阶的零点和极点
% f,m 谐振频率和耦合强度

num=size(z,2); % i阶零点数量

delta=0.01;
while (calP(z(1),z,p,f,m)*calP(z(1)*(1-delta),z,p,f,m)>0)
delta=delta*2;
end
ans(1)=fzero('calP',[z(1)*(1-delta) z(1)],[],z,p,f,m);

for i=2:num
ans(i)=fzero('calP',[z(i) z(i+1)],[],z,p,f,m);
delta=0.01;
while (calP(z(num),z,p,f,m)*calP(z(num)*(1+delta),z,p,f,m)>0)
delta=delta*2;
end
ans(num+1)=fzero('calP',[z(num) z(num)*(1+delta)],[],z,p,f,m);

function [f,m]=solveFM(z,p)
% function [f,m]=solveFM(z,p)
% function [f,m]=solveFM(z,p)
% 由输入阻抗零点和极点得到谐振频率和耦合强度
% f,m 谐振频率和耦合强度
% z,p 零点和极点位置

num=size(z,2);
if (size(p,2)~=num-1)
'Error in numbers of zeors or poles';
return;
end

zt=z;
pt=p;
f=zeros(1,num);
m=zeros(1,num-1);

for i=1:num-1
f(i)=sqrt(prod(zt.*zt)/prod(pt.*pt));
m(i)=sqrt(sum(zt.*zt)-sum(pt.*pt)-f(i)*f(i));
[zt,pt]=findNewZP(zt,pt,f(i));
end

f(num)=sqrt(prod(zt.*zt)/prod(pt.*pt));

function [z,p]=solveZP(f,m)
% function [z,p]=solveZP(f,m)
% 由谐振频率和耦合强度得到输入阻抗零点和极点
% f,m 谐振频率和耦合强度
% z,p 零点和极点位置

num=size(f,2);
if (size(m,2)~=num-1)
'Error in numbers of frequeies or coupling densities';
return;
end

zt=f(num);
pt=1;

for i=num-1:-1:1
temp=zt;
zt=findZinZeros(zt,pt,f,m);
pt=temp;
end

z=zt;
p=pt;