program is not showing graph

3 次查看(过去 30 天)
shiv gaur
shiv gaur 2022-1-11
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

回答(1 个)

Walter Roberson
Walter Roberson 2022-1-11
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 .

类别

Help CenterFile Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by