function sukh1eps=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)endendend
MATLAB: If condition returns no output
problem related to if condition in the program
Best Answer