function y=sp1ld(f)
%考慮2部分,連接處位移,力矩,剪切力，位移對r偏導(dǎo)連續(xù)，支柱自由
%kp    平面機(jī)電耦和系數(shù)
kp=0.6;
%g31
g31=-11.4*10^(-3);
%h      PZT厚度
h=3*10^(-4);
%s11d
s11d=14.4*10^(-12);
%s12d
s12d=-7.71*10^(-12);
%s33e
s33e=18.8*10^(-12);
%s33d
s33d=9.46*(10)^(-12);
%d33
d33=374*10^(-12);
%h1    黃銅皮厚度
h1=3.5*10^(-4);
%p1    黃銅皮密度
p1=8.9*10^3;
p1=8*10^3;
%p     PZT密度
p=7600;
p=8.29*10^3;
%a     PZT 簡支半徑
a=5.4*10^(-3);
%t     PZT外徑
t=a;
%e     銅片外徑
e=6.3*10^(-3);
%E 黃銅片楊式模量
E=12.3*10^(10);
%bm   黃銅片泊松比
bm=0.35;  
%V   電壓
V=10;
%以下為求得的間接值
%zm  等效電路電容與電感的阻抗和
%M   等效電路電感
%K   等效電路電順
%Uaz 純彈性能總和
%Tz  純動(dòng)能總和
%we 均值
%b
b=-s12d/s11d;
%B33T
B33T=(s33e-s33d)/d33/d33;
%B33
B33=B33T/(1-kp^2);
%V   電壓
%A
A=-h^3/(3*s11d*(1-b^2));
%B
B=g31*h^2/(2*s11d*(1-b));
%C
C=h*g31/(2*s11d*(1-b)*B33);
%At
At=-E*(h1)^3/3/(1-bm^2);
%D
D=-A-At-B*C;
%F
F=-A*b-At*bm-B*C;
%G
G=B*V/B33/h;
%L
L=At*(1-bm)+A*(1-b);
%k
k=((p1*h1+p*h)/D*(2*pi.*f)^2)^(1/4);
ka=k*a;
kt=k*t;
%D
D3=-At;
%F
F3=-At*bm;
%L
L3=At*(1-bm);
%k
k3=((p1*h1)/D3*(2*pi.*f)^2)^(1/4);
k3t=k3*t;
k3e=k3*e;
x=rand(6,6);
x(1,1)=0;
x(1,2)=0;
x(1,3)=D3*k3^2*(besselj(0,k3e)-1/k3e*besselj(1,k3e))+F3*k3/e*besselj(1,k3e);
x(1,4)=D3*k3^2*(bessely(0,k3e)-1/k3e*bessely(1,k3e))+F3*k3/e*bessely(1,k3e);
x(1,5)=-D3*k3^2*(besseli(0,k3e)-1/k3e*besseli(1,k3e))-F3*k3/e*besseli(1,k3e);
x(1,6)=D3*k3^2*(-besselk(0,k3e)-1/k3e*besselk(1,k3e))+F3*k3/e*besselk(1,k3e);
x(2,1)=0;
x(2,2)=0;
x(2,3)=-D3*k3^3*besselj(1,k3e);
x(2,4)=-D3*k3^3*bessely(1,k3e);
x(2,5)=-D3*k3^3*besseli(1,k3e);
x(2,6)=D3*k3^3*besselk(1,k3e);
x(3,1)=D*k^2*(besselj(0,kt)-1/kt*besselj(1,kt))+F*k/t*besselj(1,kt);
x(3,2)=-D*k^2*(besseli(0,kt)-1/kt*besseli(1,kt))-F*k/t*besseli(1,kt);
x(3,3)=-(D3*k3^2*(besselj(0,k3t)-1/k3t*besselj(1,k3t))+F3*k3/t*besselj(1,k3t));
x(3,4)=-(D3*k3^2*(bessely(0,k3t)-1/k3t*bessely(1,k3t))+F3*k3/t*bessely(1,k3t));
x(3,5)=-(-D3*k3^2*(besseli(0,k3t)-1/k3t*besseli(1,k3t))-F3*k3/t*besseli(1,k3t));
x(3,6)=-(D3*k3^2*(-besselk(0,k3t)-1/k3t*besselk(1,k3t))+F3*k3/t*besselk(1,k3t));
x(4,1)=-D*k^3*besselj(1,kt);
x(4,2)=-D*k^3*besseli(1,kt);
x(4,3)=-(-D3*k3^3*besselj(1,k3t));
x(4,4)=-(-D3*k3^3*bessely(1,k3t));
x(4,5)=-(-D3*k3^3*besseli(1,k3t));
x(4,6)=-(D3*k3^3*besselk(1,k3t));
x(5,1)=besselj(0,kt);
x(5,2)=besseli(0,kt);
x(5,3)=-besselj(0,k3t);
x(5,4)=-bessely(0,k3t);
x(5,5)=-besseli(0,k3t);
x(5,6)=-besselk(0,k3t);
x(6,1)=-k*besselj(1,kt);
x(6,2)=k*besseli(1,kt);
x(6,3)=-(-k3*besselj(1,k3t));
x(6,4)=-(-k3*bessely(1,k3t));
x(6,5)=-(k3*besseli(1,k3t));
x(6,6)=-(-k3*besselk(1,k3t));
x1=[0,0,-G,0,0,0]';
x2=x\x1;
A1=x2(1);
A3=x2(2);
B1=0;
B2=0;
B3=0;
B4=0;
C1=x2(3);
C2=x2(4);
C3=x2(5);
C4=x2(6);
Tz1=-pi*(p1*h1+p*h)*(2*pi*f)^2/(k^2)*...
    (A1^2*ka^2/2*(besselj(0,ka)^2+besselj(1,ka)^2)+...
     2*A1*A3*ka/2*(besseli(0,ka)*besselj(1,ka)+besselj(0,ka)*besseli(1,ka))+...
     A3^2*ka^2/2*(besseli(0,ka)^2-besseli(1,ka)^2)+...
     B1^2*(kt^2/2*(besselj(0,kt)^2+besselj(1,kt)^2)-ka^2/2*(besselj(0,ka)^2+besselj(1,ka)^2))+...
     B2^2*(kt^2/2*(bessely(0,kt)^2+bessely(1,kt)^2)-ka^2/2*(bessely(0,ka)^2+bessely(1,ka)^2))+...
     B3^2*(kt^2/2*(besseli(0,kt)^2-besseli(1,kt)^2)-ka^2/2*(besseli(0,ka)^2-besseli(1,ka)^2))+...
     B4^2*(kt^2/2*(besselk(0,kt)^2-besselk(1,kt)^2)-ka^2/2*(besselk(0,ka)^2-besselk(1,ka)^2))+...
     2*B1*B2*(kt^2/2*(besselj(0,kt)*bessely(0,kt)+besselj(1,kt)*bessely(1,kt))-ka^2/2*(besselj(0,ka)*bessely(0,ka)+besselj(1,ka)*bessely(1,ka)))+...
     2*B1*B3*(kt/2*(besseli(0,kt)*besselj(1,kt)+besselj(0,kt)*besseli(1,kt))-ka/2*(besseli(0,ka)*besselj(1,ka)+besselj(0,ka)*besseli(1,ka)))+...
     2*B1*B4*(kt/2*(besselk(0,kt)*besselj(1,kt)-besselj(0,kt)*besselk(1,kt))-ka/2*(besselk(0,ka)*besselj(1,ka)-besselj(0,ka)*besselk(1,ka)))+...
     2*B2*B3*(kt/2*(besseli(0,kt)*bessely(1,kt)+bessely(0,kt)*besseli(1,kt))-ka/2*(besseli(0,ka)*bessely(1,ka)+bessely(0,ka)*besseli(1,ka)))+...
     2*B2*B4*(kt/2*(besselk(0,kt)*bessely(1,kt)-bessely(0,kt)*besselk(1,kt))-ka/2*(besselk(0,ka)*bessely(1,ka)-bessely(0,ka)*besselk(1,ka)))+...
     2*B3*B4*(kt^2/2*(besseli(0,kt)*besselk(0,kt)+besseli(1,kt)*besselk(1,kt))-ka^2/2*(besseli(0,ka)*besselk(0,ka)+besseli(1,ka)*besselk(1,ka)))...
    );
Tz3=-pi*(p1*h1)*(2*pi*f)^2/(k3^2)*...
    (C1^2*(k3e^2/2*(besselj(0,k3e)^2+besselj(1,k3e)^2)-k3t^2/2*(besselj(0,k3t)^2+besselj(1,k3t)^2))+...
     C2^2*(k3e^2/2*(bessely(0,k3e)^2+bessely(1,k3e)^2)-k3t^2/2*(bessely(0,k3t)^2+bessely(1,k3t)^2))+...
     C3^2*(k3e^2/2*(besseli(0,k3e)^2-besseli(1,k3e)^2)-k3t^2/2*(besseli(0,k3t)^2-besseli(1,k3t)^2))+...
     C4^2*(k3e^2/2*(besselk(0,k3e)^2-besselk(1,k3e)^2)-k3t^2/2*(besselk(0,k3t)^2-besselk(1,k3t)^2))+...
     2*C1*C2*(k3e^2/2*(besselj(0,k3e)*bessely(0,k3e)+besselj(1,k3e)*bessely(1,k3e))-k3t^2/2*(besselj(0,k3t)*bessely(0,k3t)+besselj(1,k3t)*bessely(1,k3t)))+...
     2*C1*C3*(k3e/2*(besseli(0,k3e)*besselj(1,k3e)+besselj(0,k3e)*besseli(1,k3e))-k3t/2*(besseli(0,k3t)*besselj(1,k3t)+besselj(0,k3t)*besseli(1,k3t)))+...
     2*C1*C4*(k3e/2*(besselk(0,k3e)*besselj(1,k3e)-besselj(0,k3e)*besselk(1,k3e))-k3t/2*(besselk(0,k3t)*besselj(1,k3t)-besselj(0,k3t)*besselk(1,k3t)))+...
     2*C2*C3*(k3e/2*(besseli(0,k3e)*bessely(1,k3e)+bessely(0,k3e)*besseli(1,k3e))-k3t/2*(besseli(0,k3t)*bessely(1,k3t)+bessely(0,k3t)*besseli(1,k3t)))+...
     2*C2*C4*(k3e/2*(besselk(0,k3e)*bessely(1,k3e)-bessely(0,k3e)*besselk(1,k3e))-k3t/2*(besselk(0,k3t)*bessely(1,k3t)-bessely(0,k3t)*besselk(1,k3t)))+...
     2*C3*C4*(k3e^2/2*(besseli(0,k3e)*besselk(0,k3e)+besseli(1,k3e)*besselk(1,k3e))-k3t^2/2*(besseli(0,k3t)*besselk(0,k3t)+besseli(1,k3t)*besselk(1,k3t)))...
    );
Tz=Tz1+Tz3;
Uaz1=pi*D*k^2*(...
    A1^2*ka^2/2*(besselj(0,ka)^2+besselj(1,ka)^2)+...
    A3^2*ka^2/2*(besseli(0,ka)^2-besseli(1,ka)^2)-...
    2*A1*A3*ka/2*(besseli(0,ka)*besselj(1,ka)+besselj(0,ka)*besseli(1,ka)));
Uaz2=pi*L*k^2*(-A1*besselj(1,ka)+A3*besseli(1,ka))^2;
Uaz3=pi*D*k^2*(...
     B1^2*(kt^2/2*(besselj(0,kt)^2+besselj(1,kt)^2)-ka^2/2*(besselj(0,ka)^2+besselj(1,ka)^2))+...
     B2^2*(kt^2/2*(bessely(0,kt)^2+bessely(1,kt)^2)-ka^2/2*(bessely(0,ka)^2+bessely(1,ka)^2))+...
     B3^2*(kt^2/2*(besseli(0,kt)^2-besseli(1,kt)^2)-ka^2/2*(besseli(0,ka)^2-besseli(1,ka)^2))+...
     B4^2*(kt^2/2*(besselk(0,kt)^2-besselk(1,kt)^2)-ka^2/2*(besselk(0,ka)^2-besselk(1,ka)^2))+... 
     2*B1*B2*(kt^2/2*(besselj(0,kt)*bessely(0,kt)+besselj(1,kt)*bessely(1,kt))-ka^2/2*(besselj(0,ka)*bessely(0,ka)+besselj(1,ka)*bessely(1,ka)))-...
     2*B1*B3*(kt/2*(besseli(0,kt)*besselj(1,kt)+besselj(0,kt)*besseli(1,kt))-ka/2*(besseli(0,ka)*besselj(1,ka)+besselj(0,ka)*besseli(1,ka)))-...
     2*B1*B4*(kt/2*(besselk(0,kt)*besselj(1,kt)-besselj(0,kt)*besselk(1,kt))-ka/2*(besselk(0,ka)*besselj(1,ka)-besselj(0,ka)*besselk(1,ka)))-...
     2*B2*B3*(kt/2*(besseli(0,kt)*bessely(1,kt)+bessely(0,kt)*besseli(1,kt))-ka/2*(besseli(0,ka)*bessely(1,ka)+bessely(0,ka)*besseli(1,ka)))-...
     2*B2*B4*(kt/2*(besselk(0,kt)*bessely(1,kt)-bessely(0,kt)*besselk(1,kt))-ka/2*(besselk(0,ka)*bessely(1,ka)-bessely(0,ka)*besselk(1,ka)))+...
     2*B3*B4*(kt^2/2*(besseli(0,kt)*besselk(0,kt)+besseli(1,kt)*besselk(1,kt))-ka^2/2*(besseli(0,ka)*besselk(0,ka)+besseli(1,ka)*besselk(1,ka)))...
    );
Uaz4=pi*L*(-B1*k*besselj(1,kt)-B2*k*bessely(1,kt)+B3*k*besseli(1,kt)-B4*k*besselk(1,kt))^2-...
     pi*L*(-B1*k*besselj(1,ka)-B2*k*bessely(1,ka)+B3*k*besseli(1,ka)-B4*k*besselk(1,ka))^2;
Uaz5=pi*D3*k3^2*(...
     C1^2*(k3e^2/2*(besselj(0,k3e)^2+besselj(1,k3e)^2)-k3t^2/2*(besselj(0,k3t)^2+besselj(1,k3t)^2))+...
     C2^2*(k3e^2/2*(bessely(0,k3e)^2+bessely(1,k3e)^2)-k3t^2/2*(bessely(0,k3t)^2+bessely(1,k3t)^2))+...
     C3^2*(k3e^2/2*(besseli(0,k3e)^2-besseli(1,k3e)^2)-k3t^2/2*(besseli(0,k3t)^2-besseli(1,k3t)^2))+...
     C4^2*(k3e^2/2*(besselk(0,k3e)^2-besselk(1,k3e)^2)-k3t^2/2*(besselk(0,k3t)^2-besselk(1,k3t)^2))+... 
     2*C1*C2*(k3e^2/2*(besselj(0,k3e)*bessely(0,k3e)+besselj(1,k3e)*bessely(1,k3e))-k3t^2/2*(besselj(0,k3t)*bessely(0,k3t)+besselj(1,k3t)*bessely(1,k3t)))-...
     2*C1*C3*(k3e/2*(besseli(0,k3e)*besselj(1,k3e)+besselj(0,k3e)*besseli(1,k3e))-k3t/2*(besseli(0,k3t)*besselj(1,k3t)+besselj(0,k3t)*besseli(1,k3t)))-...
     2*C1*C4*(k3e/2*(besselk(0,k3e)*besselj(1,k3e)-besselj(0,k3e)*besselk(1,k3e))-k3t/2*(besselk(0,k3t)*besselj(1,k3t)-besselj(0,k3t)*besselk(1,k3t)))-...
     2*C2*C3*(k3e/2*(besseli(0,k3e)*bessely(1,k3e)+bessely(0,k3e)*besseli(1,k3e))-k3t/2*(besseli(0,k3t)*bessely(1,k3t)+bessely(0,k3t)*besseli(1,k3t)))-...
     2*C2*C4*(k3e/2*(besselk(0,k3e)*bessely(1,k3e)-bessely(0,k3e)*besselk(1,k3e))-k3t/2*(besselk(0,k3t)*bessely(1,k3t)-bessely(0,k3t)*besselk(1,k3t)))+...
     2*C3*C4*(k3e^2/2*(besseli(0,k3e)*besselk(0,k3e)+besseli(1,k3e)*besselk(1,k3e))-k3t^2/2*(besseli(0,k3t)*besselk(0,k3t)+besseli(1,k3t)*besselk(1,k3t)))...
    );
Uaz6=pi*L3*(-C1*k3*besselj(1,k3e)-C2*k3*bessely(1,k3e)+C3*k3*besseli(1,k3e)-C4*k3*besselk(1,k3e))^2-...
     pi*L3*(-C1*k3*besselj(1,k3t)-C2*k3*bessely(1,k3t)+C3*k3*besseli(1,k3t)-C4*k3*besselk(1,k3t))^2;
Uaz=Uaz1+Uaz2+Uaz3+Uaz4+Uaz5+Uaz6;
Ueaz=-pi*C*V*(-A1*ka*besselj(1,ka)+A3*ka*besseli(1,ka))...
     -pi*C*V*(-B1*kt*besselj(1,kt)-B2*kt*bessely(1,kt)+B3*kt*besseli(1,kt)-B4*kt*besselk(1,kt))...
     +pi*C*V*(-B1*ka*besselj(1,ka)-B2*ka*bessely(1,ka)+B3*ka*besseli(1,ka)-B4*ka*besselk(1,ka));
we=2/(e^2)*(...
1/(k3^2)*(C1*k3e*besselj(1,k3e)+C2*k3e*bessely(1,k3e)+C3*k3e*besseli(1,k3e)-C4*k3e*besselk(1,k3e)-C1*k3t*besselj(1,k3t)-C2*k3t*bessely(1,k3t)-C3*k3t*besseli(1,k3t)+C4*k3t*besselk(1,k3t))+...
1/(k^2)*(A1*ka*besselj(1,ka)+A3*ka*besseli(1,ka)+B1*kt*besselj(1,kt)+B2*kt*bessely(1,kt)+B3*kt*besseli(1,kt)-B4*kt*besselk(1,kt)-B1*ka*besselj(1,ka)-B2*ka*bessely(1,ka)-B3*ka*besseli(1,ka)+B4*ka*besselk(1,ka)));
%al 鋁錐輻射半徑
a0=0.02*10^(-3);
al=4.48*10^(-3);
Rl=4.95*10^(-3);
l=(Rl^2-(al-a0)^2)^(1/2);
rr=2*pi*f/340;
Rml=400*pi*al^2*(1-2*besselj(1,2*rr*al)/(2*rr*al));
Xml=400*pi*al^2*2*besselk(1,2*rr*al)/((2*rr*al)^2);
Zml=Rml+i*Xml;
Zal=Zml/(al^2);
Sl=pi*al^2;
S0=pi*a0^2;
BbA=exp(i*(-2)*rr*l)*(Zal-400/Sl)/(Zal+400/Sl);
Za0=400/S0*(Zal*cos(rr*l)+i*400/Sl*sin(rr*l))/(400/Sl*cos(rr*l)+i*Zal*sin(rr*l));
Zm0=Za0*S0^2;
n=-Ueaz*2/(we)/V;
K=2*Uaz/we^2;
M=-2*Tz/((2*pi*f)^2*we^2)+0.00002;
zm=abs(i*2*pi*f*M+K/(i*2*pi*f)+Zm0+Zml);
vj=abs(n*V/zm);
vj0=vj/(we)*(A1+A3);
A=vj0*400/(1-BbA);
pl=A*(exp(-i*rr*l)+BbA*exp(i*rr*l));
ul=pl/Zal;
vl=ul/(Sl);
%y=20*log10(2*400*2*pi*f*(we)*pi*t^2/(340/f)/2/0.3)+94;
y=20*log10(2*400*vl*pi*al^2/(340/f)/2/0.3)+94;