ricktu288 / ray-optics

A web app for creating and simulating 2D geometric optical scenes, with a gallery of (interactive) demos.
https://phydemo.app/ray-optics/
Apache License 2.0
1.31k stars 220 forks source link

Gradient-index optics #69

Closed ricktu288 closed 1 year ago

ricktu288 commented 2 years ago

There have been several gradient-index optics examples made by @StasFainer:

They are made by and simulated with brute force stacking of a lot of single-index refractors, so the performance is very poor and a lot of artifacts are generated.

I think there should be a much better and efficient way to simulate them, such as having a tool that takes an arbitrary n(x,y) and simulate by the following: If a beginning point p1 of a ray is inside the relevant region, then rayIntersection will always return a point p such that the distance between p1 and p is some small distance eps. And then shot will give the best estimation of the new direction of the ray.

However, I am not very familiar with optics higher than the level of freshman physics (my research area is condensed matter theory). So I don't know the proper numerical algorithm to estimate the direction of the ray based on n(x,y) so that the numerical error is minimized. Also, I don't know how to deal with the reflected rays if the function is not smooth in general.

It would be great if someone familiar in this area can create such a tool, or at least provide a proper algorithm.

StasFainer commented 2 years ago

This should be much more efficient, and a really interesting tool.

Regarding how to estimate the direction of the ray based on n(x,y) - the path of the ray is governed by the Eikonal equation, which can be solved numerically, from what I read. Maybe I'll try to make such a tool in the future, and in the meantime I hope that this suggestion helps somewhat.

StasFainer commented 1 year ago

@ricktu288 I'm considering trying to implement this tool, and wanted to ask if you worked on it thus far? (so that we are not repeating the same work)

ricktu288 commented 1 year ago

Great! I haven't worked on it. Please implement that.

StasFainer commented 1 year ago

@ricktu288 I attached a zip below, with the code and four examples – the difference between the four "circlelens_grin.js" files, is only in the properties of the object returned by the "create" function.

For now, I only implemented one numerical solver (the simple Euler method).

In addition, for now, the "draw" function in circlelens_grin.js works as if the grin lens has a constant refractive index, equal to 2.

Before I implement additional features, such as surface merging, wavelength dependent refractive index, other glass forms (half-plane, free-shape, …) etc., I wanted to ask what do you think about the current solution/what changes should be made to it?

By the way, about the implementation of the toolbar for this tool, I thought about adding an additional user text input to the circle glass toolbar, which will contain the refractive index function, and something like a checkbox, that will switch between the two modes – grin refractive index, and constant refractive index. What are your thoughts about the toolbar implementation for this tool?

Finally, there was an issue with ensuring (in an “if” statement) the condition that a point is inside/outside a circle, or on its boundary, due to numerical inaccuracies, which forced me to add the “eps” property, and it seems to solve it. I think it might be a good idea to add a user input for the value of the "eps" and "step_size" properties, in the toolbar for this tool, because the range of useful values for these properties, may depend on the scenario.

Simply put, the condition for ensuring that a point ray.p1 is inside the circle “obj” (the case for outside the circle is very similar), is: (graphs.length_squared(obj.p1, ray.p1) - R_squared - obj.eps < 0) && (graphs.length_squared(obj.p1, ray.p1) - R_squared + obj.eps < 0)

And for ensuring that a point ray.p1 is on the boundary of the circle “obj”, the condition is:

(graphs.length_squared(obj.p1, ray.p1) - R_squared - obj.eps < 0) && (graphs.length_squared(obj.p1, ray.p1) - R_squared + obj.eps > 0)

It’s worth noting, that in practice, it seems that the tool works properly even when omitting/changing some of these conditions, however, I obviously can’t check every possible scenario, and I feel more confident with the current conditions in the code.

And speaking of the "step_size" property, the length of the first ray segment which is refracted by the "refract" function (from the "refractor" object), is of length 1 (line 861 in "refractor.js"), so I guess that in general it will be better that a grin object (like "circlelens_grin") will have its own "refract" function, where the length of this refracted segment is "step_size" (specifically in the attached code, step_size = 1), do you agree?

code+examples.zip

Luneburg lens - luneburg

Maxwell fisheye lens - fisheye

Logarithmic spiral - logarithmic_spiral

Constant refractive index (n=2) comparison - n=2

ricktu288 commented 1 year ago

I think the current solution inside the glass is OK, but the reflected rays on the surface are missing currently. One complication is that if the user enters a discontinuous function, in theory there should also be reflected rays from the discontinuity (and need surfaceMerging. But note that currently surfaceMerging is also not done properly for the custom equation glass anyway). I don't know how easy is to implement this properly.

Using an additional option for the circlelens tool may be a good idea (I think that should be an "advanced option" which is hidden by default, as in the beam tool.) However, if you also add the rectangular GRIN slab tool, there is no direct constant-index counterpart in the toolbar currently. Combining the circular one with an existing tool but isolating the rectangular one seems inconsistent. One way to solve this is to add a rectangular glass which is constant-index by default. Also, I think the equation field should be consistent with the "custom equation" tool, that is, use MathQuill.

The "length" of the "ray" object has no meaning, so it should not matter what size you set. Setting it to the step size is a good idea, but the rayIntersection must work even if the length is not the step size (for example, the ray may interact with another tool inside the

ricktu288 commented 1 year ago

If the refractive index is not restricted to having a particular symmetry, then maybe add the options to the custom shape/equation tool is also fine, but the coordinates may be more complicated.

Also, I think the derivatives should be calculated automatically. Maybe math.js is a good solution for symbolic derivative. Note however that it appear to not integrated with MathQuill natively. One way to bridge these two libraries is to use tex-math-parser.

StasFainer commented 1 year ago

In the Luneburg and Maxwell fisheye lenses examples, the reflected rays are missing, as they should be, because in these two specific examples, there is no discontinuity in the refractive index – the radius of these lenses according to the latex equation in the “p” property, is 100, but if you decrease the actual radius of the circle in the simulator, below 100, the reflected rays will appear due to a discontinuity in the refractive index, as in the picture below –

image

You can also see reflected rays in the picture of the fourth example (where n=2).

About what I wrote regarding the “step_size” property, I previously thought that out of all the ray segments which are drawn inside the circle, the first is of length 1 and the rest are of length step_size, which is not consistent, because if we set the numerical step accuracy to be step_size, the length of all the ray segments inside the circle, should be equal to step_size. However, the “rayIntersection” function in the “circlelens_grin” object takes care of it (all the ray segments inside the circle, are of length step_size, from what I checked), so no problem there.

I’ll keep working on it, and maybe the next thing that I’ll implement is another GRIN lens shape.

StasFainer commented 1 year ago

The new refractor_grin tool, is a simple polygon GRIN lens.

Note however, that the current implementation also has the circular arcs functionality (as the original refractor tool), which shouldn’t be used as part of this new tool – I guess I can take out the circular arcs functionality from the refractor_grin tool, but still one should know to use this new tool with simple polygons only (not self-intersecting and without holes), as the current GRIN implementation supports only that.

Below is the GRIN slab example, using this new tool.

I added the “origin” property to both of the current GRIN lenses (it represents the origin for the refractive index function n(x,y)), which the user should input into the toolbar (which I’ll implement in the future).

I also made changes to the circlelens_grin.js structure – so that the main differences between different GRIN lenses (so far – circle and simple polygon GRIN lenses), are in the implementation of the functions – isInsideGlass, isOutsideGlass and isOnBoundary (the shot and rayIntersection functions in the new tool, are “inherited” from the circlelens_grin tool).

The code is simpler now, but obviously the price is that it’s not as efficient. However, I’m not sure if the alternative, more efficient solution, will make the execution noticeably faster, for many examples. In any case, I’ll leave it as is for now, and if needed, more efficient implementations could be done in the future (for example, by implementing custom shot and rayIntersection functions for specific GRIN lenses, instead of using the generic inherited ones).

I also added auxiliary functions to circlelens_grin.js and refractor_grin.js, to make the code more readable.

Now, I’ll try to implement surface merging for these two GRIN lenses.

GRIN_slab_JSON.zip

GRIN_slab

StasFainer commented 1 year ago

Unfortunately, I partially began working on the grin surface merging, just this week. I'll try to finish it by the next weekend.

StasFainer commented 1 year ago

I added the “body merging/overlapping” capability (overlapping grin objects function as expected now). “surface merging” works as before (no additional code was needed for that).

Below is an example of seismic P wave(rays) shadow zones(for more information, you can see this). This is only to demonstrate the body merging function, and is only qualitative (I approximated the mantle, outer core and inner core [white disc in the example below] refractive indices, by the seismic wave velocity from here, using a piecewise linear function – it will probably be easier to do, once I’ll implement the toolbar for these grin lenses). Note that I added a circular blocker near the boundary of the mantle (from the inside), to eliminate many reflected waves, otherwise the shadow zones are not as pronounced (the P wave shadow zones, by definition [from what I understood], refer to direct P waves only, without reflections/refractions – see the previous clip). Still, some rays enter the shadow zones in this example, due to reflections from the core.

image

The next two examples also demonstrate the body merging function: The refractive index of the grin circle lens, is the following (arbitrary) positive and bounded function: $1+exp\left(-x^2\right)$, and the grin rectangle has a refractive index equal to the reciprocal of the previous function. As can be seen from the picture with them overlapping, versus the picture with the rectangle alone, the overlapping region acts like air (n=1), as expected.

image

image

Now, I’ll implement the toolbar for these grin lenses (including symbolic derivatives, etc.).

PS – I thought about maybe adding an absorption feature for this tool, in the future, which I noticed someone requested today:

seismic_rays + bodyMerging examples.zip

JeremyLeaf commented 1 year ago

PS – I thought about maybe adding an absorption feature for this tool, in the future, which I noticed someone requested today.

That would be amazing if you could! Thanks for considering it 🙂. I was thinking that I was incorrect in my original post about the gradient intensity of the beam. It would actually be an exponential decay or something within the glass. Might be a bit harder to plot.

StasFainer commented 1 year ago

Update - I wasn't able to work on the toolbar until now. I'll try to finish it until the next weekend (or the weekend afterwards).

StasFainer commented 1 year ago

I added the GRIN circle and free-shape objects to the toolbar.

The integration of the tex to math parser was annoying - it works in a Node environment rather than a web environment, so certain changes had to be made to this parser. The thing is that prior to this, I didn't know anything about these subjects, and after trying different solutions (like using module bundlers, such as Webpack), the solution I went with was to simply make manual changes to the source code - which turned out to be not many modifications.

Since now the partial derivatives are done automatically, I added below another example of a GRIN lens, from this article. Note that this "cloaking lens" has a refractive index with a singularity at its center - the function itself appears in (3.1) and (3.2) in this article - so like in the logarithmic spiral example, I added a circular blocker around the center. Also, the default step size had to be changed from 1 to 0.1 to accurately simulate this example.

image cloaking_lens.zip

Tell me if you want me to change/add something before deploying these tools.

In addition, I just now noticed a problem with the body merging functionality - the current implementation creates a bug when there're multiple GRIN objects with refractive index functions that are not defined everywhere (this problem doesn't exist without this functionality). In the following example, the circle is a Luneburg lens (its refractive index is defined only in a circle with radius 100) and the rectangle has a constant refractive index of 1. I know where the problem is, and I'll think about how to fix this (modify the current implementation of the body merging function, or replace the whole implementation).

image

ricktu288 commented 1 year ago

I think one thing to add is to draw the glass according to the refractive index function, rather than a uniform gray color (note that for the original non-GRIN tools, even if you stack multiple glasses, the resulting gray color always corresponds to the effective refractive index). I don't know if that will affect the performance a lot. If it becomes too slow, a workaround may be to just draw the border or use a new kind of representation (so at least not conflicting with existing visual representation). The current representation is confusing because it makes people think the inside is uniform and the border is always discontinuous.

As for body merging, i think it may not be as important as surface merging, so if it is too difficult to be fully implemented, you may just add a warning, etc., as you currently do for the requirement of simple polygon

ricktu288 commented 1 year ago

圖片 Also, a newly created GRIN glass appears not working properly.

StasFainer commented 1 year ago

圖片 Also, a newly created GRIN glass appears not working properly.

It does work, you just need to set the origin of n(x,y) accordingly (copy the coordinates of the center of your lens, to the origin of n(x,y) input box).

In the beginning, I thought about maybe setting the origin of n(x,y) for the grin circle lens, to be equal to the center of the lens once it's created, in the 'create' function itself, but then it's not that clear what to do in the free-shape case, since there isn't a "clear" choice for the center in that case. Also, should the origin of the grin circle lens be constantly updated with the center of the grin circle lens? I think it shouldn't..

In any case, I can set the origin of n(x,y) for the grin circle lens, to be equal to the center of the lens once it's created, such that consequent modifications of the center will not effect the origin of n(x,y), but then the user should remember to manually change the origin of n(x,y) accordingly, after each modification of the center of the existing grin circle lens.

StasFainer commented 1 year ago

I think one thing to add is to draw the glass according to the refractive index function, rather than a uniform gray color (note that for the original non-GRIN tools, even if you stack multiple glasses, the resulting gray color always corresponds to the effective refractive index). I don't know if that will affect the performance a lot. If it becomes too slow, a workaround may be to just draw the border or use a new kind of representation (so at least not conflicting with existing visual representation). The current representation is confusing because it makes people think the inside is uniform and the border is always discontinuous.

As for body merging, i think it may not be as important as surface merging, so if it is too difficult to be fully implemented, you may just add a warning, etc., as you currently do for the requirement of simple polygon

Sure, I'll try it, check the performance, and continue accordingly.

About the body merging, I think it's important because it gives you flexibility to create different shapes, that you are unable to create otherwise (without adding new tools every time), like the concentric circles example of seismic rays (which would require an annulus tool without the body merging). I think it's still worth trying to fix the bug with body merging.

ricktu288 commented 1 year ago

It does work, you just need to set the origin of n(x,y) accordingly (copy the coordinates of the center of your lens, to the origin of n(x,y) input box).

Maybe a better way is to use a default n(x,y) which is defined on the entire plane so that a new user will not get confused.

StasFainer commented 1 year ago

It does work, you just need to set the origin of n(x,y) accordingly (copy the coordinates of the center of your lens, to the origin of n(x,y) input box).

Maybe a better way is to use a default n(x,y) which is defined on the entire plane so that a new user will not get confused.

Yes, that might be better. I can set the default to simply be n=1.5 like in the circle lens tool.

ricktu288 commented 1 year ago

Or maybe something like a periodic n(x,y) (is there a physically relevant example?) so that a new user can play with it non-trivially without need to enter an equation.

ricktu288 commented 1 year ago

Sure, I'll try it, check the performance, and continue accordingly.

One way to improve the performance (if the direct pixel-wise drawing is too slow) is to draw it on a virtual canvas and cache the image, and that every time draw() is called just draw that image on the canvas (also note the current "hack" of automatic red color for absolute n<1).

Also note that the use of getImageData/putImageData may interfere with the hardware accelerated drawing on some hardware and may decrease performance and image quality (e.g. rasterization bias) in the entire simulator.

I guess using getImageData/putImageData on a virtual canvas may be OK (without decreasing performance/quality), but I never tried that before.

ricktu288 commented 1 year ago

I estimated the performance of the point-wise calculation (for drawing the gray color) on an nxn grid. However, the performance starts to become poor for n>100. I think if we stick with drawing an n=100 grid then some cases (especially with singularity) the appearance on screen may become low-quality.

Although caching the image will help, this does not prevent re-calculation when the user change the shape or possibly adjusting the parameters if #105 is going to be implemented.

Possible solutions are using WebGL, workers (also very complicated), or the kind of setInterval drawing as what is currently done for ray tracing, but those are all quite completed and need to rewrite a lot of core code.

Therefore, I think the best solution may be to just draw with another appearance without showing the index. I suggest filling with a constant dark sky-blue color.

So I think what's left is to fix the body-merging bug.

ricktu288 commented 1 year ago

Adaptive sampling may help, but I don't know how complicated it is to implement that. Since both of us does not have much free time, I think the best way is to make a minimal implementation and merge the PR. Making the appearance intuitive may be left for future improvement.

ricktu288 commented 1 year ago

For the examples in the Gallery, I suggest putting two versions together, one is a rough approximation using a series of constant refractor where the distance between the surfaces is human-readable, and another version is using the GRIN refractor.

If #105 is to be implemented, then in the future the rough approximation one may be replaced by an interactive, adjustable numbers or slices.

StasFainer commented 1 year ago

Sure, I'll work on fixing the body-merging bug first and foremost. I'm for the minimal implementation using constant filling, since like you said, both of us don't have much free time at the moment. After that, I'll add the examples to the gallery (including some that I posted in the current thread).

StasFainer commented 1 year ago

To fix the body merging bug, I currently added another implementation of the body merging functionality, based on symbolic computation, as opposed to the previous implementation which was mainly numerical.

The grin objects have a toggle to switch between these two implementations – like I wrote in the note popover, the symbolic implementation is slower but doesn't have this bug(note that this bug arises in very specific scenarios, as written in the note popover).

I considered other numerical implementations, as to fix this bug without the slower performance. However, every time I found a specific scenario in which this implementation fails (with regards to this bug), whereas the symbolic one remains robust (the slowness of the symbolic implementation stems, from what I checked, from the use of the symbolic derivative function – 'math.derivative' ).

I changed the constant color of the GRIN objects – feel free to change it or to suggest a different color.

The grin circlelens object is now initialized with a specific gaussian like function, and the grin refractor with a specific periodic function. Both of these functions are defined in the entire plane, and the numerical error associated with ray tracing them, can be handled within the prescribed limits of the step size (0.1 to 1), from what I checked.

image image

examples.zip

The periodic sine and cosine based refractive index functions, seem to be quite susceptible to numerical errors using the current Euler method. It would be interesting to add additional numerical algorithms in the future, such that the user could choose one from a defined list, that best suits his/her needs.

ricktu288 commented 1 year ago

Great! I will merge the PR in about 1~2 weeks. Thanks for your contribution!

StasFainer commented 1 year ago

No problem :) I'll add the examples to the gallery after you merge the PR (unless you prefer otherwise).

ricktu288 commented 1 year ago

Now I've merged the PR.

For the examples, I think my previous suggestion of putting together a human-readable rough approximation alone with the GRIN glass one may not be good in some cases, especially those susceptible to numerical inaccuracy. So just put whatever you think is suitable.