program is not working
1 view (last 30 days)
Show older comments
function shiv5
K0= (2*pi/1530)*1e6:(2*pi/10)*1e6:(2*pi/1630)*1e6;
for j=1:numel(K0)
k0 = K0(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,k0) - f(p0,k0))/h1;
DELTA2 = (f(p2,k0) - f(p1,k0))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,k0)*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,k0)/E;
p = p2 + h;
if abs(h) < TOL
disp(p)
break
end
p0 = p1;
p1 = p2;
p2 = p;
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,k0) - f(p0,k0))/h1;
DELTA2 = (f(p2,k0) - f(p1,k0))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=i+1;
end
P(j)=real(p);
R(j)=real(r)
end
plot(K0,R)
end
function y=f(x,k0)
n0=1.3707;
n1=1.3;
n2=1.59;
n3=1.45;
n4=3.46;
na=1.36;
t1=2.10e-6;
t2=0.198e-6;
t3=0.0012e-6;
t4=0.002e-6;
%x=1.36-1i*0.0123;
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= 1i*(g0*M11/n0^2+ga*M22/na^2)-M21+(g0*ga/n0^2*na^2)*M12 ;
r=(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);
end
5 Comments
Accepted Answer
More Answers (1)
KSSV
on 6 Jan 2022
This line will not creat any array, it is empty. So loop will not run and the required variables in the plot command are not defined.
K0= (2*pi/1530)*1e6:(2*pi/10)*1e6:(2*pi/1630)*1e6;
Consider replacing the line with:
K0= linspace((2*pi/1530)*1e6,(2*pi/1630)*1e6,10);
But again the variable r, in the line:
R(j)=real(r)
is not defined anywhere. Please rethink on your code.
2 Comments
See Also
Categories
Find more on Adding custom doc in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!