Incident polarisation implementation may be unphysical

[golovart]

Results for simulations with rotations of incident polarisation seem to be physically questionable. As a small demonstration:

Two silver spheres (50nm radius, complex refraction index) at a large distance (bigger than the wavelength).

s1 = Sphere(r = 0.05, center = (0, -0.3, 0))
s2 = Sphere(r = 0.05, center = (0, 0.3, 0))
detector = hp.detector_grid(100, .01)
bisph = Scatterer(lambda point: np.logical_or(s1.contains(point), s2.contains(point)), n = 0.05+1j*2.1, center = (0.5, 0.5, 0))
len_dda = Lens(1.38, DDA())

# changing bounds by hand, since there is nothing in the center
# and default algorithm will fail
bond = np.array(bisph.bounds)
c = bisph.center
for i,b in enumerate(bond):
    bond[i,0] = min(s1.bounds[i][0], s2.bounds[i][0])
    bond[i,1] = max(s1.bounds[i][1], s2.bounds[i][1])
bisph.indicators.bound = bond

fig, axes = plt.subplots(1, 3, figsize=(23, 6))
axes = np.ravel(axes)
for i,th in enumerate([0,45,90]):
    th *= np.pi/180
    holo = calc_field(detector, bisph, medium_index=1., illum_wavelen=0.4, illum_polarization=(np.cos(th), np.sin(th)), theory=len_dda)
    holo = (np.abs(holo)**2).sum(dim='vector')
    holo.plot(ax=axes[i])
    axes[i].set_title(str(np.around(th*180/np.pi))+' degree')
plt.show()

The spheres are always on the same position, but rotation of polarisation rotates the image!

A possible reason for that can be coming from mieangfuncs.calc_scat_field. It has a step of rotating incident (Ex, Ey) before multiplying with amplitude scattering matrix.

The physical reason for that is unclear, since default perpendicular and parallel polarisations are along X and Y.

If we parametrize polarisation vector with alpha - it looks like it is similar to rotating phi on the corresponding angle. The transformation from mieangfuncs.incfield:

Therefore, this might be the reason image rotates, when we change the incident polarisation vector.

Why do we need the mieangfuncs.incfield step in the mieangfuncs.calc_scat_field?

[jeromefung]

In the formulation of Mie scattering we're using, the amplitude scattering matrix acts on components of the incident and scattered fields relative to the scattering plane (which is defined as containing the scatterer and the detector location). So we absolutely need to take the spherical coordinate phi of the detector location into account to the polarization vector (E_x, E_y) in the lab frame to components (E_par, E_perp) relative to the scattering plane, at least in order to implement Lorenz-Mie scattering correctly. This consideration also applies to the multisphere code, which is formulated in a similar way.

Whether doing this is correct in the case of connecting to ADDA is a different question, however, and I haven't worked with ADDA in HoloPy in detail. I would agree that something looks fishy with what you showed.

Best, Jerome

On Fri, Sep 25, 2020 at 11:56 AM golovart notifications@github.com wrote:

Results for simulations with rotations of incident polarisation seem to be physically questionable. As a small demonstration:

Two silver spheres (50nm radius, complex refraction index) at a large distance (bigger than the wavelength).

s1 = Sphere(r = 0.05, center = (0, -0.3, 0)) s2 = Sphere(r = 0.05, center = (0, 0.3, 0)) detector = hp.detector_grid(100, .01) bisph = Scatterer(lambda point: np.logical_or(s1.contains(point), s2.contains(point)), n = 0.05+1j*2.1, center = (0.5, 0.5, 0)) len_dda = Lens(1.38, DDA())

changing bounds by hand, since there is nothing in the center

and default algorithm will fail

bond = np.array(bisph.bounds) c = bisph.center for i,b in enumerate(bond): bond[i,0] = min(s1.bounds[i][0], s2.bounds[i][0]) bond[i,1] = max(s1.bounds[i][1], s2.bounds[i][1]) bisph.indicators.bound = bond

fig, axes = plt.subplots(1, 3, figsize=(23, 6)) axes = np.ravel(axes) for i,th in enumerate([0,45,90]): th *= np.pi/180 holo = calc_field(detector, bisph, medium_index=1., illum_wavelen=0.4, illum_polarization=(np.cos(th), np.sin(th)), theory=len_dda) holo = (np.abs(holo)*2).sum(dim='vector') holo.plot(ax=axes[i]) axes[i].set_title(str(np.around(th180/np.pi))+' degree') plt.show()

[image: image] https://user-images.githubusercontent.com/32298047/94285294-c6267900-ff53-11ea-868b-300829a1c07e.png

The spheres are always on the same position, but rotation of polarisation rotates the image!

A possible reason for that can be coming from mieangfuncs.calc_scat_field. It has a step of rotating incident (Ex, Ey) before multiplying with amplitude scattering matrix.

The physical reason for that is unclear, since default perpendicular and parallel polarisations are along X and Y.

If we parametrize polarisation vector with alpha - it looks like it is similar to rotating phi on the corresponding angle. The transformation from mieangfuncs.incfield: [image: image] https://user-images.githubusercontent.com/32298047/94288152-63cf7780-ff57-11ea-9e5b-5ba1dc31e776.png

Therefore, this might be the reason image rotates, when we change the incident polarisation vector.

Why do we need the mieangfuncs.incfield step in the mieangfuncs.calc_scat_field?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/manoharan-lab/holopy/issues/371, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB3MN4CNLSJJWRMCZ5IGQHDSHS4UFANCNFSM4RZWQ4DA .

[jeromefung]

(I should add that the incfield subroutine directly implements the 1st two equations on p. 62 of Bohren & Huffman, which has a figure illustrating the definition of the scattering plane.)

Jerome

On Fri, Sep 25, 2020 at 2:19 PM Jerome Fung jerome.fung@gmail.com wrote:

In the formulation of Mie scattering we're using, the amplitude scattering matrix acts on components of the incident and scattered fields relative to the scattering plane (which is defined as containing the scatterer and the detector location). So we absolutely need to take the spherical coordinate phi of the detector location into account to the polarization vector (E_x, E_y) in the lab frame to components (E_par, E_perp) relative to the scattering plane, at least in order to implement Lorenz-Mie scattering correctly. This consideration also applies to the multisphere code, which is formulated in a similar way.

Whether doing this is correct in the case of connecting to ADDA is a different question, however, and I haven't worked with ADDA in HoloPy in detail. I would agree that something looks fishy with what you showed.

Best, Jerome

On Fri, Sep 25, 2020 at 11:56 AM golovart notifications@github.com wrote:

Results for simulations with rotations of incident polarisation seem to be physically questionable. As a small demonstration:

Two silver spheres (50nm radius, complex refraction index) at a large distance (bigger than the wavelength).

s1 = Sphere(r = 0.05, center = (0, -0.3, 0)) s2 = Sphere(r = 0.05, center = (0, 0.3, 0)) detector = hp.detector_grid(100, .01) bisph = Scatterer(lambda point: np.logical_or(s1.contains(point), s2.contains(point)), n = 0.05+1j*2.1, center = (0.5, 0.5, 0)) len_dda = Lens(1.38, DDA())

changing bounds by hand, since there is nothing in the center

and default algorithm will fail

bond = np.array(bisph.bounds) c = bisph.center for i,b in enumerate(bond): bond[i,0] = min(s1.bounds[i][0], s2.bounds[i][0]) bond[i,1] = max(s1.bounds[i][1], s2.bounds[i][1]) bisph.indicators.bound = bond

fig, axes = plt.subplots(1, 3, figsize=(23, 6)) axes = np.ravel(axes) for i,th in enumerate([0,45,90]): th *= np.pi/180 holo = calc_field(detector, bisph, medium_index=1., illum_wavelen=0.4, illum_polarization=(np.cos(th), np.sin(th)), theory=len_dda) holo = (np.abs(holo)*2).sum(dim='vector') holo.plot(ax=axes[i]) axes[i].set_title(str(np.around(th180/np.pi))+' degree') plt.show()

[image: image] https://user-images.githubusercontent.com/32298047/94285294-c6267900-ff53-11ea-868b-300829a1c07e.png

The spheres are always on the same position, but rotation of polarisation rotates the image!

A possible reason for that can be coming from mieangfuncs.calc_scat_field . It has a step of rotating incident (Ex, Ey) before multiplying with amplitude scattering matrix.

The physical reason for that is unclear, since default perpendicular and parallel polarisations are along X and Y.

If we parametrize polarisation vector with alpha - it looks like it is similar to rotating phi on the corresponding angle. The transformation from mieangfuncs.incfield: [image: image] https://user-images.githubusercontent.com/32298047/94288152-63cf7780-ff57-11ea-9e5b-5ba1dc31e776.png

Therefore, this might be the reason image rotates, when we change the incident polarisation vector.

Why do we need the mieangfuncs.incfield step in the mieangfuncs.calc_scat_field?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/manoharan-lab/holopy/issues/371, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB3MN4CNLSJJWRMCZ5IGQHDSHS4UFANCNFSM4RZWQ4DA .

[barkls]

It turns out the behaviour isn't because of DDA, but the new Lens class. Here's proof: Top row uses Lens(Multisphere), bottom row uses Multisphere only. Mie seems better behaved but I don't understand why. Any ideas @ralex0 or @briandleahy ?

code to generate image:

import numpy as np
from matplotlib import pyplot as plt

import holopy as hp
from holopy.scattering import Sphere, Spheres, Multisphere, Mie, calc_intensity
from holopy.scattering.theory import Lens # this should probably be moved up to scattering namespace

detector = hp.detector_grid(100, .01)

s1 = Sphere(r = 0.05, center = (0.5, 0.8, 3), n=0.05+1j*2.1)
s2 = Sphere(r = 0.05, center = (0.5, 0.2, 3), n=0.05+1j*2.1)
scatterer = Spheres([s1, s2])

no_lens_theory = Multisphere()
lens_theory = Lens(1.4, no_lens_theory)

fig, axes = plt.subplots(2, 3, figsize=(20, 10))
axes = np.ravel(axes)
for i, th in enumerate([0, 45, 90]):
    for j, theory in enumerate([lens_theory, no_lens_theory]):
        theta = th * np.pi/180
        polarization=(np.cos(theta), np.sin(theta))
        holo = calc_intensity(detector, scatterer, medium_index=1., illum_wavelen=0.4,
                                      illum_polarization=polarization, theory=theory)
        holo.plot(ax=axes[i+j*3])
        axes[i].set_title(str(np.around(th))+' degree')
plt.savefig('test.png')
plt.show()

[briandleahy]

Sorry for the delayed response. This looks like a bug to me.

The reason that the code has to take two separate pathways is that a lens rotates the polarization of light. So the Lens class needs to perform some modifications to the coordinates to calculate a hologram. That being said, the polarization rotation is a small effect. It does not cause the image to rotate. So this looks like a bug.

I'll take a look at it quickly and see (i) if the same problem appears for Mielens and (ii) if I can pinpoint where the problem is.

Thanks for pointing this out @golovart

[briandleahy]

Some of this may be due to the pixel scaling being used in the demonstration. When I check this while (i) zooming out and (ii) using the Lens theory & Mie superpostiion, I don't see any rotation of the image, although if you look closely you will see some changes near the center, which are presumably due to polarization effects.

import numpy as np
from matplotlib import pyplot as plt

import holopy as hp
from holopy.scattering import (
    Sphere, Spheres, Multisphere, Mie, calc_holo, calc_intensity)
from holopy.scattering.theory import Lens, MieLens

# units in um
detector = hp.detector_grid(100, 0.1)

s1 = Sphere(r=0.05, center=(4.5, 4.8, 3), n=0.05+1j*2.1)
s2 = Sphere(r=0.05, center=(5.5, 4.2, 3), n=0.05+1j*2.1)
scatterer = Spheres([s1, s2])

no_lens_theory = Mie()  # Multisphere()
lens_theory = Lens(1.4, no_lens_theory)
mielens = MieLens(1.4)

fig, axes = plt.subplots(3, 3, figsize=(20, 10))
axes = np.ravel(axes)
for i, th in enumerate([0, 45, 90]):
    for j, theory in enumerate([lens_theory, no_lens_theory, mielens]):
        theta = th * np.pi/180
        polarization=(np.cos(theta), np.sin(theta))
        holo = calc_intensity(
            detector, scatterer, medium_index=1.0, illum_wavelen=0.4,
            illum_polarization=polarization, theory=theory)
        ax=axes[i+j*3]
        ax.imshow(holo.values.squeeze(), vmin=0)

        title = '{}: {} degree'.format(theory.__class__.__name__, int(th))
        ax.set_title(title)
        ax.set_xticks([])
        ax.set_yticks([])

plt.show()
plt.tight_layout()
plt.savefig('mie-mielens-lensmie.png')

However, the superposition pathway for the calculation is different from the Multisphere or DDA pathway. @barkls, can I ask you to check this zoomed out, with Multisphere or DDA? Right now Mutisphere crashes on my machine and it will take me a little while to fix.

[barkls]

I uncommented the Multisphere in your code and here is the result:

The difference in behaviour between Mie and Multisphere is surely a clue into what's going wrong here...

[briandleahy]

@golovart, this will take me a little while to look into. In the meanwhile, I just want to provide a few suggestions so this doesn't slow you down too much:

• As @barkls has shown, this is only a problem with the new Lens class. You can ignore the effect of the lens by just using DDA() instead of Lens(DDA()).
• The Multisphere calculations are much faster than the DDA calculations for spheres. Multisphere should be more accurate as well.

[golovart]

Well, unfortunately, I am working on simulation of the images of non-default objects (so DDA()) in the focal plane of optical microscope (so Lens()).

But can we be sure that, for instance, the output for polarisation (1,0) is correct? Because then, as a temporary workaround, one can rotate scatterer before passing to ADDA, obtain the result and then rotate the output image back. Does that make any sense? At least while there is an issue with polarisation in Lens()

Another (probably quite stupid/basic) physics question: can we simulate unpolarised incident light with HoloPy? Or would we need to calculate fields for a bunch of different polarisation angles and average the results? (And therefore it would become computationally extremely expensive and long)

Thank you @barkls and @briandleahy for looking into this problem with Lens!

[jeromefung]

We haven't implemented scattering of unpolarized incident light in HoloPy. But it would not be too difficult to do: all of the scattering theories (Lorenz-Mie, multisphere, and DDA) ultimately calculate the amplitude scattering matrix. The scattering response to unpolarized light is given by the S_11 element of the 4 x 4 scattering matrix, which (up to a factor of 1/2) the sum of the squared magnitudes of the 4 elements of the 2 x 2 amplitude scattering matrix. (See Bohren & Huffman, Ch. 3 for details.) Consequently, there is no need to calculate fields for many different polarizations and numerically average the results.

On Tue, Oct 6, 2020 at 7:06 PM golovart notifications@github.com wrote:

Well, unfortunately, I am working on simulation of the images of non-default objects (so DDA()) in the focal plane of optical microscope (so Lens()).

But can we be sure that, for instance, the output for polarisation (1,0) is correct? Because then, as a temporary workaround, one can rotate scatterer before passing to ADDA, obtain the result and then rotate the output image back. Does that make any sense? At least while there is an issue with polarisation in Lens()

Another (probably quite stupid/basic) physics question: can we simulate unpolarised incident light with HoloPy? Or would we need to calculate fields for a bunch of different polarisation angles and average the results? (And therefore it would become computationally extremely expensive and long)

Thank you @barkls https://github.com/barkls and @briandleahy https://github.com/briandleahy for looking into this problem with Lens!

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/manoharan-lab/holopy/issues/371#issuecomment-704598758, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB3MN4BK3FJC2MHYFF65UYTSJOPFPANCNFSM4RZWQ4DA .

[briandleahy]

@golovart -- I am pretty confident that the polarization output for (1,0) is correct. @ralex0 and I checked it numerically against the MieLens theory, which has been a little more thoroughly tested.

Regarding unpolarized incident light -- it's more or less as @jeromefung says, with the additional caveat that the lens changes the polarization, and that the Lens API doesn't directly expose the scattering matrices. That being said, I think you'll be able to evaluate the image under unpolarized light with 1-2 calls to the Lens model. For instance, the image of a sphere under unpolarized light can be calculated as the mean of the image under x-polarized and y-polarized light. I'm not sure of the top of my head if this generalizes to complex scatterers though. I would take a look at Eq. (1) in the Lens paper, which comes from Eq. (6) in the MieLens paper and see if you can show that it holds for general scatterers.

Hope that helps!

[briandleahy]

The gross problem with polarization rotation is now fixed on my fork, under the "fixlenspolarizationrotation" branch. (And the calculation for x-polarized light was correct.)

However, in the process of fixing this, I found that the polarization rotates slightly differently under Lens(Mie) and Mie(Lens). This is captured under the test test_raw_fields_similar_mielens_ypolarizationon the "fixlenspolarizationrotation" branch.

Both theories agree when the polarization is in the x-direction, but when the incident polarization is in the y-direction, the x-component scattered fields differ by a factor of -1. Chasing this down it seems that, in the Lens theory, the perpendicular component of the scattered fields rotates differently than in the MieLens theory. It will take me a little while to figure out why this is the case, which version is correct, and then to solve the issue.

@ralex0 , any thoughts? I can provide a MWE if you'd like.

[golovart]

I discovered another counter-intuitive behavior of simulation with Lens. If I just take a small sphere and look at the pixel with peak intensity, I would expect it to be in the same point (in XY plane) as the center of the sphere. And ADDA without Lens gives this result. However, with Lens it appears significantly shifted and I am not sure about the reason for such symmetry break.

s = Sphere(n=0.05+1j*2.1, r=0.025, center=(1,1,0.5))
detector = hp.detector_grid(200, .01)
len_dda = Lens(1.4, DDA(use_gpu=True, max_dpl_size=0.002))
holo = calc_field(detector, s, medium_index=1., illum_wavelen=0.4, illum_polarization=(1,0), theory=len_dda)
holo = (np.abs(holo)**2).sum(dim='vector')
print('center with lens:', np.unravel_index(np.argmax(holo), (200,200)))

no_len_dda = DDA(use_gpu=True, max_dpl_size=0.002)
holo = calc_field(detector, s, medium_index=1., illum_wavelen=0.4, illum_polarization=(1,0), theory=no_len_dda)
holo = (np.abs(holo)**2).sum(dim='vector')
print('center with no lens:', np.unravel_index(np.argmax(holo), (200,200)))

The output is (in pixels, while detector pixel is 0.01): center with lens: (100, 119) center with no lens: (100, 100)

And this is the default X-polarisation.

[ralex0]

@briandleahy yeah can you post or share a mwe so I can see what you mean?

[briandleahy]

@golovart -- yes, that seems correct. If you look at the holograms, you will see that the brightest pixels are symmetricly located across the x-direction, away from the center of the scatterer, on a set of rings. The rings are there because of diffraction effects through the lens. The central spot is dimmer because the sphere is slightly defocused, by about 1 wavelength. The defocus gives an interference pattern which is darker in the center than elsewhere. The position of the brightest pixel shifts in x as the particle shifts in z; you can see that when the particle is in focus, the brightest pixel is at the hologram's center. There is a bit of a discussion in this in our Mie + Lens paper, albeit for the hologram (including the interference with the incident beam, which you aren't including here) instead of for the intensity.

So I think that the lens calculation is correct here. (Although I have not yet merged in the fix for the polarization rotation.)

[golovart]

@briandleahy yeah, that was a bad example. In fact, my real question was about a sphere in the focus (and I tried rotating it by passing addacmd=['-orient',str(theta),'0','0']). And the problem there was that the peak was fluctuating around center by no more that 1 pixel (Yet it was strange that symmetry was slightly violated). But in focus it is problematic to compare with no lens. So I decided to try putting it out of focus to be able to spot whether it's caused by lens.

Anyway, I tried running with your fix and now peak is always in the center for a sphere in the focus.

[briandleahy]

@ralex0 I've been looking into this, and I think the problem is a little more subtle than I was saying before. You can find a MWE on my branch, under holopy.scattering.tests.test_lens.TestLens.test_transforms_correctly_with_polarization_rotation. I've pasted the MWE code here; it should be runnable:

from holopy.scattering.tests.test_lens import *

# We test that rotating the lab frame correctly rotates
# the polarization.
# If we rotate (x0, y0) -> (y1, -x1), then the polarization
# in the new coordinates should be
# E1x = E0y, E1y = -E1x
scatterer = test_common.sphere
medium_wavevec = 2 * np.pi / test_common.wavelen
medium_index = test_common.index
theory = Lens(
            lens_angle=LENS_ANGLE,
            theory=Mie(False, False),
            quad_npts_theta=200,
            quad_npts_phi=200,
            )

krho = np.linspace(0, 100, 11)
phi_0 = 0 * krho
phi_1 = np.full_like(krho, -np.pi / 2)
kz = np.full_like(krho, 20.0)

pol_0 = xr.DataArray([1.0, 0, 0])
pos_0 = np.array([krho, phi_0, kz])

pol_1 = xr.DataArray([0, -1.0, 0])
pos_1 = np.array([krho, phi_1, kz])

args = (scatterer, medium_wavevec, medium_index)

fields_0 = theory._raw_fields(pos_0, *args, pol_0)
fields_1 = theory._raw_fields(pos_1, *args, pol_1)

tols = {'atol': 1e-5, 'rtol': 1e-5}
assert_allclose(fields_1[0],  fields_0[1], **tols)
assert_allclose(fields_1[1], -fields_0[0], **tols)

You can see the problem more clearly by plotting:

import matplotlib.pyplot as plt
f = plt.figure()
ax1 = f.add_subplot(1, 2, 1)
ax1.plot(fields_1[1].real); ax1.plot(-fields_0[0].real); ax1.set_title("Transformed Raw fields")
ax2 = f.add_subplot(1, 2, 2)
ax2.plot(fields_1[1].real + fields_0[0].real); ax1.set_title("Difference; should be 0")

You'll notice that the difference is on the order of ~1-5/100, when it should be 0. At first I thought this was numerical. However, it does not go away as the number of quadrature points increases. I thought that perhaps it was an issue with the phi quadrature, so I improved the phi quadrature a bit (using trapz-like points, which are optimal for integrals over a preiod, instead of Gauss-Legendre quadrature) but that didn't fix anything.

One clue is that the difference between Lens(Mie) and MieLens is zero for x-polarized light but a few percent for y-polarized light. So it's something with the polarization rotation, that is small.....

Thoughts?

[barkls]

This should be fixed by #373. @briandleahy can we close this?

[briandleahy]

Yes. The issues in this thread are fixed by #373