program is not working

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

[shiv gaur]

pl help to plot the graph between k0 vs R I have some problem that ie relate with programming

[Image Analyst]

I know saying simply "it's not working" is one of the "don'ts" in the guidelines. Please read them (again):

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[KSSV]

Why do you think it is not working? Any error? What is your problem?

[shiv gaur]

sir value of r is not showing plot is also pproblem

[shiv gaur]

if any one help

Answer 1

[shiv gaur]

this is the answer

function shiv5

%K0= (2*pi/1530)*1e6:(2*pi/10)*1e6:(2*pi/1630)*1e6;

K0= linspace((2*pi/1530)*1e6,(2*pi/1630)*1e6,10);

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(2*pi./K0,P)

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);

r1=conj(r);

y=r*r1;

end

Answer 2

[KSSV]

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.

[shiv gaur]

the value of r from above equation ie

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);

if we take value the value of x any one then plot b/w k0 vs r

pl help

[Walter Roberson]

You define r within f but you do not return it from f.