If condition returns no output

Question

Sukhversha Luthra 2018-1-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/375324-if-condition-returns-no-output

评论： Sukhversha Luthra 2018-1-3

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