from running code find diff plot vs variable

1 view (last 30 days)
function mm1
clear all
close all
n=linspace(1.43,1.50,10);
D1= linspace(1e-9,3e-7,15);
P=zeros(numel(n),numel(D1));
for j=1:numel(D1)
for k=1: numel(n)
da = D1(j);
na=n(k);
p0 = 1;
p1 = 1.5;
p2 = 2;
TOL = 10^-8;
N0 = 100; format long
h1 = p1 - p0;
h2 = p2 - p1;
DELTA1 = (f(p1,na,da) - f(p0,na,da ))/h1;
DELTA2 = (f(p2,na,da ) - f(p1,na,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=3;
while i <= N0
b = DELTA2 + h2*d;
D = (b^2 - 4*f(p2,na,da )*d)^(1/2);
if abs(b-D) < abs(b+D)
E = b + D;
else
E = b - D;
end
h = -2*f(p2,na,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,na,da ) - f(p0,na,da ))/h1;
DELTA2 = (f(p2,na,da ) - f(p1,na,da ))/h2;
d = (DELTA2 - DELTA1)/(h2 + h1);
i=i+1;
end
end
%pv(j)=p;
P(j)=real(p);
%R(j)=real(r);
end
plot(D1,P)
end
function y=f(x,na,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;
%kx=(2*pi/l)*np*sin(a);
kc=k0*sqrt(x^2-ec);
ka=k0*sqrt(ea-x^2);
kf=k0*sqrt(ef-x^2);
km=k0*sqrt(em-x^2);
ks=k0*sqrt(x^2-es);
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
pl help to plot between dy/dna vs da

Answers (1)

Walter Roberson
Walter Roberson on 28 Jan 2022
n=linspace(1.43,1.50,10);
function y=f(x,na,da)
ea=na.^2;
n will eventually be 1.50 because of the linspace(), so ea = na.^2 will eventually be 2.25
ka=k0*sqrt(ea-x^2);
Suppose that x happens to be 1.5... perhaps because
p1 = 1.5;
then ea-x^2 would be 1.50^2 - 1.5^2 --> 0 so ka would become 0
rac=(ka-kc)/(ka+kc);
ka is 0, so rac is -kc/kc which will be -1 unless kc is exactly 0.
f=(1-rac*exp(2*1i*da*ka))/(1+rac*exp(2*1i*da*ka));
Remember that ka is 0, so 2*1i*da*ka is 0 unless da is infinite or nan. exp(0) is 1. Remember that rac is -1 because it is -kc/kc when ka is 0. So you now have -1*1 which is 1. And then you have 1+(-1) which is 0. So you have a division by 0. Which gives you NaN, which poisons the rest of your calculation.
P(j)=real(p);
You need to double-check that. You created P as a 2D array

Tags

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by