If condition returns no output

1 vue (au cours des 30 derniers jours)
Sukhversha Luthra
Sukhversha Luthra le 3 Jan 2018
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

Birdman
Birdman le 3 Jan 2018
Replace this
if abs(DE)==1.1921e-07
disp(hj)
end
with this
if abs(DE)-1.1921e-07<1e-4
disp(hj)
end
  1 commentaire
Sukhversha Luthra
Sukhversha Luthra le 3 Jan 2018
Thanks Birdman. It worked.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Gamma Functions dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by