If condition returns no output
Afficher commentaires plus anciens
function sukh1
eps=0.;
beta_T=2e-4;
lambda=77.6;mu=38.6;C_e=383;k=400; rho=8920;
T_0=318;tau_0=1e-12; tau_1=2*tau_0;
bulk=lambda+2*mu;
m=5;
omg=2*pi*(2^m);
alpha_T=(3*lambda+2*mu)*beta_T;
A=bulk*k;
B=rho*C_e*(bulk)*(tau_0+1i/omg)+rho*k+(alpha_T^2)*T_0*(1-1i*omg(1-eps)*tau_1)*(eps*tau_0+1i/omg);
C=(rho*rho*C_e*(tau_0+i/omg));
%Instead of alpha, alpha^2 is used everywhere. so we will directly find
%alpha^2 and denote it by alf2.
alf2=roots([C -B A]);
%dialational wave with velocity alpha1(alpha2)is called qP(q(T))
nu1=(alpha_T*T_0*omg*omg*(eps*tau_0+i/omg))/(rho*C_e*alf2(1)*(tau_0+i/omg)-k);
nu2=(alpha_T*T_0*omg*omg*(eps*tau_0+i/omg))/(rho*C_e*alf2(2)*(tau_0+i/omg)-k);
eta=nu1/nu2;
%gm stands for gamma;
gm1=bulk/rho*alf2(1);
gm2=bulk/rho*alf2(2);
a=(gm1-eta*gm2); b=1-eta;
%beta=sqrt(mu/rho), bets=beta^2
bets=(mu/rho);
e1=bets/alf2(1);e2=bets/alf2(2);
d=a*b-(e1+eta*eta*e2);
g=a*a+b*b+4*d; es=e1+e2; ep=e1*e2;
beta=sqrt(bets);
C0=1024*eta*(d+eta*es);
C1=256*(d*d-eta*(g+4*d))-1024*eta*eta*(ep+2*es);
C2=128*(2*eta*a*(a+b)-(d-2*eta)*g)+1024*(eta^2)*(2*ep+es);
C3=16*((g^2)+8*a*(a+b)*(d-2*eta)-4*eta*(a^2))-1024*(eta^2)*ep;
C4=-32*a*(a*(d-2*eta)+(a+b)*g);
C5=8*(a^2)*(3*g+4*a*b-8*d);
C6=-8*(a^3)*(a+b);
C7=a^4;
h=roots([C7 C6 C5 C4 C3 C2 C1 C0]);
for j=1:7;
hj=h(j);
%pv is the phase velocity denoted by c in the paper and h=(c^2)/bets
pvs=hj*bets; pv=sqrt(pvs);
q1=sqrt((pvs/alf2(1))-1);
q2=sqrt((pvs/alf2(2))-1);
q3=sqrt(hj-1);
DE=(2-hj)*((2-gm1*hj)-eta*(2-gm2*hj))+4*(q1-eta*q2)*q3;
abs(DE);
V=beta*abs(hj)/real(sqrt(hj));
if abs(DE)==1.1921e-07;
disp(hj)
end
end
end
Réponse acceptée
Plus de réponses (0)
Catégories
En savoir plus sur Gamma Functions dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!