Multiparticle Mie C: only fixed polarization implemented

Could you double check this fix? mie.forward_holo returns weird holograms for incident fields that are polarized in the y-direction. Fields that are x-polarized seem to be ok. The fit in test_mie_noisysingle (which assumes y-polarization) now fails.

I'll check it, but for fitting based on the input decks and analyze/fit.py framework, the mie_fortran rather than C-based mie is currently enabled.

yes, but I changed the code to use mie rather than mie_fortran prior to running the fits. Without the polarization correction introduced by the fix, the fits converge to values similar to those given by mie_fortran. With the polarization correction, they do not converge at all. This was just an example of why the proposed fix does not seem to work. For a simpler example, see the attachments below. These were produced with this code:

import numpy as np
import os
import string
import pylab

from nose.tools import raises, assert_raises
from numpy.testing import assert_, assert_equal, assert_array_almost_equal
from nose.tools import with_setup

import holopy
from holopy.model.theory import mie
#from holopy.model.theory import mie_fortran as mie
from holopy.model.theory import Mie
from holopy.model.scatterer import Sphere, SphereCluster
from holopy.third_party import nmpfit
from holopy.process import normalize

# define optical train
wavelen = 658e-9
polarization = [0., 1.0] # y-polarized
#polarization = [1.0, 0.] # x-polarized
divergence = 0
pixel_scale = [.1151e-6, .1151e-6]
index = 1.33

optics = holopy.optics.Optics(wavelen=wavelen, index=index,
                              pixel_scale=pixel_scale,
                              polarization=polarization,
                              divergence=divergence)

# initial guess
scaling_alpha = .6
radius = .85e-6
n_particle_real = 1.59
n_particle_imag = 1e-4
x = .576e-05
y = .576e-05
z = 15e-6

imshape = 128

holo = mie.forward_holo(imshape, optics, n_particle_real,
                            n_particle_imag, radius, x, y, z,
                            scaling_alpha)

It turns out this was a particularly nasty bug, and one that is probably related to why the superposition code was failing sometimes in the fit. Jerome's fix for the polarization is correct, but there was a problem in the code. You have to be very careful about using the /= construction. If you pass an array to a function and then do a /= (or *= or +=) operation, you modify the original array!

In [24]: arr = numpy.ones(4)
In [25]: def f(x):
   ....: x /= 2
   ....:
In [26]: print arr
[ 1. 1. 1. 1.]
In [27]: f(arr)
In [28]: print arr
[ 0.5 0.5 0.5 0.5]

Note that this happens even if you assign a new name to the variable, because this just creates another reference back to the original array:

In [29]: arr = numpy.ones(4)
In [30]: def f(x):
   ....: y = x
   ....: y /= 2
   ....:
In [31]: print arr
[ 1. 1. 1. 1.]
In [32]: f(arr)
In [33]: print arr
[ 0.5 0.5 0.5 0.5]

So I would recommend generally doing y = x.copy() in functions to avoid accidentally modifying the argument to the function. Use /= with extreme caution when it's operating on an argument passed to the function.

The mie.py code for forward_holo(..., x, y, ...) contained lines like

    x /= opt.pixel[0]
    y /= opt.pixel[1]

which ended up modifying the original arrays passed to the function. I think this was responsible for the fitter going awry. It may also be why I was getting some strange holograms for different polarizations and different phase factors.

Fixes coming soon.

There is still a bug in the way the mie cython code computes holograms for y-polarized light. I still get images that look like the attachments (see the sidebar of the bug page for the pictures). The file holopy/tests/test_mie shows the code used to calculate these holograms.

It looks like the problem is that MieFieldExtension.c assumes the incident light is polarized in the x-direction. So Jerome's fix won't work, as the problem is internal to the c code. We could rotate the fields produced by MieFieldExtension, but unless the incident field is purely x- or y-polarized, we'd have to calculate the fields over a larger grid than the image size to be able to rotate it properly.

I would propose that instead of fixing MieFieldExtension, we simply deprecate it and work with mie-fortran. We will need to add a function to mieangfuncs.f90 to return the fields, but once that's in place superposition should be straightforward. Any comments?

Marking this bug as WONTFIX.

Looks like that is probably the best option at this point.

We could rewrite the C code so that it doesn't assume a fixed
polarization (which would probably be a better option than rotating
after the fact), but that would be a fair chunk of work.

The only minor concern with all this fortran code is I don't know who
will maintain it after Jerome graduates, but that shouldn't stop us from
implementing what looks like the right solution at the moment.

On 06/16/2011 05:15 PM, Vinothan N. Manoharan wrote:
> It looks like the problem is that MieFieldExtension.c assumes the
> incident light is polarized in the x-direction. So Jerome's fix won't
> work, as the problem is internal to the c code. We could rotate the
> fields produced by MieFieldExtension, but unless the incident field is
> purely x- or y-polarized, we'd have to calculate the fields over a
> larger grid than the image size to be able to rotate it properly.
>
> I would propose that instead of fixing MieFieldExtension, we simply
> deprecate it and work with mie-fortran. We will need to add a function
> to mieangfuncs.f90 to return the fields, but once that's in place
> superposition should be straightforward. Any comments?
>
> Marking this bug as WONTFIX.
>
> ** Changed in: holopy/stable
> Status: New => Won't Fix
>
> ** Changed in: holopy/dev
> Status: New => Won't Fix
>

OK -- I will take care of tweaking mieangfuncs.f90.

I've added code to mieangfuncs.f90 to directly output the fields. There's now a new function, mie.calc_mie_fields_f(), with identical calling syntax to mie.calc_mie_fields() but with the calculations based on my Fortran codes.

In so doing, I refactored some of mieangfuncs.f90 for clarity and to reduce repetition. I also improved the numerical robustness mieangfuncs.getsphercoords at small theta.

There remain differences between calc_mie_fields() and calc_mie_fields_f(). While trying to track down the differences, I corrected one probable error in void *flds() in MieFieldExtension.c (line 338), where the spherical coordinate phi was being calculated as arctan(x/y). I haven't successfully thoroughly understood void *flds() to know whether or not there are additional issues in that code; for that reason i have not substituted calc_mie_fields_f() for calc_mie_fields() in mie.forward_holo().

There is one difference (at the 0.05% level) between mie.forward_holo() and mie_fortran.forward_holo() stemming from how the codes handle the E^2 hologram term. mie calculates |(Ex, Ey, Ez)|^2 whereas mie_fortran calculates |(Ex, Ey)|^2. My argument for not including E_z^2 was that the camera's pixels should be sensitive to the z component of the Poynting vector, E x B, and the z component of E x B cannot depend on Ez.

mie.py in dev now supports arbitrary polarization through the Fortran code. In stable, mie.py has been renamed to mie_cython.py, which does not support arbitrary polarization. The mie_fortran module does support it, and it now enabled by default for Mie superposition calculations.