Refractive index profile calculation

Hello all,

I'm tryng to calculate a Refractive index profile from an article (Ref 1):

But my code does not reproduce smothly as it was suppose to do. Goes to 1.48 (last value as the curve above) but as a square:

The thing is that as you can see in my code below I have N layers (this means N values of refrective index) and when I change the values below 1000 I get some weird shapes. The right curve (first curve attached in this text) has 1000 of values, I mean they (arthurs of that curve) used N = 1000 because converge to a continuos curve. I think I'm very near to the response but I don't know what could be wrong.

If I change N = 10, for example, look the shape:

I mean is not continuos, of course, because it is no the right value that the artours used to had that curve but with N = 1000 it is like a square. And I'm using all the equations to calculate right. Check this:

1. From Ref 1 they said:

2. And the figure that I'm studying to reproduce is (From Ref. 1):

As you can see in my code I calculate theta correct and also Lf up to Lf maximum (this case is 10um (as you can see in the article curve and in my code)).

I really appreciate any help !!

NOTE: THE ARTICLE HAS OPEN ACESS FROM OSA WEBSITE. IT IS NOT RESTRICTED SO ANYONE CAN SEE.

References

1. Qian Wang, Yingyan Huang, Ter-Hoe Loh, Doris Keh Ting Ng, and Seng-Tiong Ho, "Thin-film stack based integrated GRIN coupler with aberration-free focusing and super-high NA for efficient fiber-to-nanophotonic-chip coupling," Opt. Express 18, 4574-4589 (2010).

[madhan ravi]

It’s better post your code as text than as an image.

[darova]

I don't any expression for n

only and

[darova]

According to these formulas

n should be this

for i = 2:n
    nn = n(1:i-1)/n(i);
    Lf = h/sum(sqrt(nn.^2-1));
    n(i) = Lf/sqrt(Lf^2+h^2)*n(i-1);
end

[Stephen23]

Note to Tay: All content you submit to this forum is covered by the Creative Commons Attribution Share Alike 3.0 license:

https://www.mathworks.com/matlabcentral/content/terms-of-use.html#answers

so you have already given everyone on the planet permission to copy and distribute your question.

Original question retrieved from Google Cache:

https://webcache.googleusercontent.com/search?q=cache:gb3GyY1cR6oJ:https://www.mathworks.com/matlabcentral/answers/525655-refractive-index-profile-calculation

https://web.archive.org/web/20211124061959/https://webcache.googleusercontent.com/search?q=cache%3Agb3GyY1cR6oJ%3Ahttps%3A%2F%2Fwww.mathworks.com%2Fmatlabcentral%2Fanswers%2F525655-refractive-index-profile-calculation

https://web.archive.org/web/20211124062650/https://www.mathworks.com/matlabcentral/answers/525655-refractive-index-profile-calculation

Refractive index profile calculation

Hello all,

I'm tryng to calculate a Refractive index profile from an article (Ref 1):

NOTE: THE ARTICLE HAS OPEN ACESS FROM OSA WEBSITE. IT IS NOT RESTRICTED SO ANYONE CAN SEE.

References

1. Qian Wang, Yingyan Huang, Ter-Hoe Loh, Doris Keh Ting Ng, and Seng-Tiong Ho, "Thin-film stack based integrated GRIN coupler with aberration-free focusing and super-high NA for efficient fiber-to-nanophotonic-chip coupling," Opt. Express 18, 4574-4589 (2010).

[Tay]

Thanks Stephen :))

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[David Goodmanson]

Hi Tay,

The following code addresses (1) in the text. For Lf = 1, D = 1 and 1000 layers the resulting plot reproduces the one in the text pretty well.

When D > Lf you can get some funky results with the refractive index being less than 1.

Lf = 1;
D = 1;
N = 1000;
h = D/N;
Lfh = Lf/h;
n1 = 3.5;       % refractive index of layer 1
nvec = n1;      % vector of refractive indices

for k = 1:N
  ni = fzero(@(ni) fun(ni,nvec,Lfh),[.1, nvec(end)-1e-10]);   % have to be careful with upper limit
  nvec = [nvec ni];
end
plot(nvec)
grid on

function y = fun(ni,nvec,Lfh)
  y = sum(1./sqrt(nvec.^2/ni^2-1)) - Lfh;
end

[David Goodmanson]

Hi Tay,

the idea is that you want to find the zero of a function involving the sum of terms of the form

(1./sqrt(nk.^2/ni^2-1))

That means that ni has to be smaller than all of the previously calculated nk, which are stored in nvec. If ni were larger than any of them, then in the expression above you would end up taking the square root of a negative number.

The numbers in nvec are decreasing, so ni will be less than all of them if it's less than nvec(end). The e-10 makes it smaller by just a little bit, so .1 and nvec(end)-1e-10 are suitable lower and upper limits and the value of ni does end up being in between those limits. The fzero function finds a zero crossing between the two limits and can find the zero.

[David Goodmanson]

Hello Tay,

could you tell me why you want the answer to be deleted?

[Walter Roberson]

The question and answer were available for 18 months; what remaining damage could there be?

[David Goodmanson]

Hello Tay,

Deleting a question that has been answered is not considered accceptable. One reason is that the site is a source of information and solutions, and deleting a question takes away potentially useful information. A more personal reason is that having put in time on the answer, I would not like to see it go away. In this case the question concerns verifying a process involving Snell's law from a freely available 11 year old paper, so I do not yet understand what harm there might be. One idea is, if you wanted to rephrase the question in some manner, I would certainly be amenable to modifying the answer accordingly, and then deleting the comments that have to do with this situation. Any thoughts?

[Tay]

Thank you David. this was supposed to be simple, but it started to be an forum

no problem, it's ok hahaha