program is not showing graph
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
function metal3
theta1=linspace(10,70,20);
%T3=1e-9:1e-9:1e-6;
for j=1:numel(theta1)
theta=theta1(j);
p0 = 0.5;
p1 = 1;
p2 = 1.5;
TOL = 10^-8;
N0 = 100;
format long;
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,theta) - f(p0,theta))/h1;
DELTA2 = (f(p2,theta) - f(p1,theta))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,theta)*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,theta)/E;
p = p2 + h;
if abs(h) < TOL
p
disp(p)
break
end
p0 = p1;
p1 = p2;
p2 = p;
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,theta) - f(p0,theta))/h1;
DELTA2 = (f(p2,theta) - f(p1,theta))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i = i + 1;
end
if i > N0
formatSpec = string('The method failed after N0 iterations,N0= %d \n');
end
P(j)=real(p);
end
plot(theta ,r)
end
function y=f(x,theta)
k0=(2*pi/0.6328)*1e6;
n0=1.3707;
n1=1.30;
n2=1.59;
n3=0.065+4*1i;
n4=3.48;
na=1.36;
t1=2.10e-6;
t2=0.198e-6;
t3=1e-6;
t4=0.002e-6;
x=n0*k0*sin(theta);
k1=k0*sqrt(n1.^2-x.^2);
k2=k0*sqrt(n2.^2-x.^2);
k3=k0*sqrt(n3.^2-x.^2);
k4=k0*sqrt(n4.^2-x.^2);
m11= cos(t1*k1)*cos(t2*k2)-(k2/k1)*sin(t1*k1)*sin(t2*k2);
m12=(1/k2)*(cos(t1*k1)*sin(t2*k2)*1i) +(1/k1)*(cos(t2*k2)*sin(t1*k1)*1i);
m21= (k1)*cos(t2*k2)*sin(t1*k1)*1i +(k2)*cos(t1*k1)*sin(t2*k2)*1i;
m22=cos(t1*k1)*cos(t2*k2)-(k1/k2)*sin(t1*k1)*sin(t2*k2);
m34= cos(t3*k3)*cos(t4*k4)-(k4/k3)*sin(t3*k3)*sin(t4*k4);
m32=(1/k4)*(cos(t3*k3)*sin(t4*k4)*1i) +(1/k3)*(cos(t4*k4)*sin(t3*k3)*1i);
m23= (k3)*cos(t4*k4)*sin(t3*k3)*1i +(k4)*cos(t3*k3)*sin(t4*k4)*1i;
m33=cos(t3*k3)*cos(t4*k4)-(k3/k4)*sin(t3*k3)*sin(t4*k4);
M11=m11*m34+m12*m23;
M12=m11*m32+m12*m33;
M21=m21*m34+m22*m23;
M22=m21*m32+m22*m33;
g0=sqrt(x.^2-n0.^2)*k0;
ga= sqrt(x.^2-na.^2)*k0;
y=(na^2*g0*M11-n0^2*ga*M22+g0*ga*M12-na^2*n0^2*M21)/(na^2*g0*M11+n0^2*ga*M22+g0*ga*M12+na^2*n0^2*M21);
r=y.^2;
end
Réponses (1)
Walter Roberson
le 11 Jan 2022
You have
function y=f(x,theta)
You are passing in a variable named x
k0=(2*pi/0.6328)*1e6;
That's a value in the millions
n0=1.3707;
n1=1.30;
That's a value near 1
n2=1.59;
n3=0.065+4*1i;
n4=3.48;
na=1.36;
t1=2.10e-6;
t2=0.198e-6;
t3=1e-6;
t4=0.002e-6;
x=n0*k0*sin(theta);
You have not used the x that you passed in. Instead you have overwritten x with a value that is determined in part by k0, which is in the millions, so x is going to end up in the millions.
k1=k0*sqrt(n1.^2-x.^2);
Remember n1 is near 1, and you are subtracting a million squared, and taking the square root. So you are going to end up with a value that is roughly complex 1 million.
m11= cos(t1*k1)*cos(t2*k2)-(k2/k1)*sin(t1*k1)*sin(t2*k2);
cos of complex 1 million is infinite . Remember that cos of an imaginary number has to do with exponentials .
0 commentaires
Voir également
Catégories
En savoir plus sur Loops and Conditional Statements 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!