Info

This question is closed. Reopen it to edit or answer.

plot is showing blank

1 view (last 30 days)
shiv gaur
shiv gaur on 7 Jan 2022
Closed: Cris LaPierre on 7 Jan 2022
function shiv7
%K0= (2*pi/1530)*1e6:(2*pi/10)*1e6:(2*pi/1630)*1e6;
%K0= linspace((2*pi/1530)*1e6,(2*pi/1630)*1e6,5);
D1=linspace(1,300,10);
for j=1:numel(D1)
da = D1(j);
p0 = 1;
p1 = 1.5;
p2 = 2;
TOL = 10^-8;
N0 = 100; format long
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,da) - f(p0,da ))/h1;
DELTA2 = (f(p2,da ) - f(p1,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,da )*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,da)/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,da ) - f(p0,da ))/h1;
DELTA2 = (f(p2,da ) - f(p1,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=i+1;
end
P(j)=real(p);
%R(j)=real(r);
end
plot(D1*1e-9,P)
end
function y=f(kx,da )
l=6330*(10)^(-10);
k0=2*pi/l;
nc=1.33;
ec=nc^2;
na=1.43;
ea=na^2;
nf=1.59;
ef=nf^2;
nm=0.064+1i*4;
em=nm^2;
ns=1.495;
es=ns^2;
df=0.5e-6;
dm=0.03e-6;
%da=0.06e-6;
%kx=(2*pi/l)*np*sin(a);
kc=sqrt(kx^2-ec*k0^2);
ka=sqrt(ea*k0^2-kx^2);
kf=sqrt(ef*k0^2-kx^2);
km=sqrt(em*k0^2-kx^2);
ks=sqrt(kx^2-es*k0^2);
rfm=(kf-km)/(kf+km);
rms=(km-ks)/(ks+km);
rfa=(kf-ka)/(kf+ka);
rac=(ka-kc)/(ka+kc);
a=(1-rfm)/(1+rfm);
b=(1-rms*exp(2*1i*dm*km))/(1+rms*exp(2*1i*dm*km));
c=(1-rfa)/(1+rfa);
f=(1-rac*exp(2*1i*da*ka))/(1+rac*exp(2*1i*da*ka));
pfms=2*atan(1i*a*b);
pfac=2*atan(1i*c*f);
y=2*kf*df-pfms-pfac;
end
  2 Comments
shiv gaur
shiv gaur on 7 Jan 2022
pl help for plot and the change in na/change in kx vs da
Star Strider
Star Strider on 7 Jan 2022
Please always look at the variables the code is creating and that are to be plotted —
%K0= (2*pi/1530)*1e6:(2*pi/10)*1e6:(2*pi/1630)*1e6;
%K0= linspace((2*pi/1530)*1e6,(2*pi/1630)*1e6,5);
D1=linspace(1,300,10);
for j=1:numel(D1)
da = D1(j);
p0 = 1;
p1 = 1.5;
p2 = 2;
TOL = 10^-8;
N0 = 100; format long
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,da) - f(p0,da ))/h1;
DELTA2 = (f(p2,da ) - f(p1,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,da )*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,da)/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,da ) - f(p0,da ))/h1;
DELTA2 = (f(p2,da ) - f(p1,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=i+1;
end
pv(j) = p;
P(j)=real(p);
%R(j)=real(r);
end
pv % Display 'p' Vector
pv =
NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi NaN + NaNi
p_stats = [nnz(isnan(pv)), numel(pv)] % [Number Of 'NaN' Values, Length Of Vector]
p_stats = 1×2
10 10
plot(D1*1e-9,P)
function y=f(kx,da )
l=6330*(10)^(-10);
k0=2*pi/l;
nc=1.33;
ec=nc^2;
na=1.43;
ea=na^2;
nf=1.59;
ef=nf^2;
nm=0.064+1i*4;
em=nm^2;
ns=1.495;
es=ns^2;
df=0.5e-6;
dm=0.03e-6;
%da=0.06e-6;
%kx=(2*pi/l)*np*sin(a);
kc=sqrt(kx^2-ec*k0^2);
ka=sqrt(ea*k0^2-kx^2);
kf=sqrt(ef*k0^2-kx^2);
km=sqrt(em*k0^2-kx^2);
ks=sqrt(kx^2-es*k0^2);
rfm=(kf-km)/(kf+km);
rms=(km-ks)/(ks+km);
rfa=(kf-ka)/(kf+ka);
rac=(ka-kc)/(ka+kc);
a=(1-rfm)/(1+rfm);
b=(1-rms*exp(2*1i*dm*km))/(1+rms*exp(2*1i*dm*km));
c=(1-rfa)/(1+rfa);
f=(1-rac*exp(2*1i*da*ka))/(1+rac*exp(2*1i*da*ka));
pfms=2*atan(1i*a*b);
pfac=2*atan(1i*c*f);
y=2*kf*df-pfms-pfac;
end
The ‘p’ variable consists entirely of NaN values, and NaN values cannot be plotted. You need to discover the reason for this, and correct the code.
.

Answers (2)

shiv gaur
shiv gaur on 7 Jan 2022
pl help for ploting
  2 Comments
shiv gaur
shiv gaur on 7 Jan 2022
graph should be like this fig
shiv gaur
shiv gaur on 7 Jan 2022
1.graph is showing blank
2.pl plot dkx/dna=change in kx/change in na vs da value
any one pl help

shiv gaur
shiv gaur on 7 Jan 2022
now it is working
function shiv7
%K0= (2*pi/1530)*1e6:(2*pi/10)*1e6:(2*pi/1630)*1e6;
%K0= linspace((2*pi/1530)*1e6,(2*pi/1630)*1e6,5);
%D1=1e-9:1e-9:3e-7;
D1= linspace(1e-9,3e-7,15);
for j=1:numel(D1)
da = D1(j);
p0 = 1;
p1 = 1.5;
p2 = 2;
TOL = 10^-8;
N0 = 100; format long
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,da) - f(p0,da ))/h1;
DELTA2 = (f(p2,da ) - f(p1,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,da )*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,da)/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,da ) - f(p0,da ))/h1;
DELTA2 = (f(p2,da ) - f(p1,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=i+1;
end
%pv(j)=p;
P(j)=real(p)
%R(j)=real(r);
end
k0=2*pi/0.6328*10^-10;
plot(D1,P/k0)
end
function y=f(kx,da )
l=6330*(10)^(-10);
k0=2*pi/l;
nc=1.33;
ec=nc^2;
na=1.43;
ea=na^2;
nf=1.59;
ef=nf^2;
nm=0.064+1i*4;
em=nm^2;
ns=1.495;
es=ns^2;
df=0.5e-6;
dm=0.03e-6;
%da=0.06e-6;
%kx=(2*pi/l)*np*sin(a);
kc=sqrt(kx^2-ec*k0^2);
ka=sqrt(ea*k0^2-kx^2);
kf=sqrt(ef*k0^2-kx^2);
km=sqrt(em*k0^2-kx^2);
ks=sqrt(kx^2-es*k0^2);
rfm=(kf-km)/(kf+km);
rms=(km-ks)/(ks+km);
rfa=(kf-ka)/(kf+ka);
rac=(ka-kc)/(ka+kc);
a=(1-rfm)/(1+rfm);
b=(1-rms*exp(2*1i*dm*km))/(1+rms*exp(2*1i*dm*km));
c=(1-rfa)/(1+rfa);
f=(1-rac*exp(2*1i*da*ka))/(1+rac*exp(2*1i*da*ka));
pfms=2*atan(1i*a*b);
pfac=2*atan(1i*c*f);
y=2*kf*df-pfms-pfac;
end

This question is closed.

Tags

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!