Use anisotropic refractive index with FD modesolver

[DanHickstein]

I would like to use an anisotropic refractive index with the FD modesolver. This seems to basically be implemented, but I can't figure out a practical way to make my "epsfunc" return a 2D array of length-5 tuples. To me, it seems more practical to have the function return a NxMx5 numpy array, but this is not quite compatible with the current code (seems easy to modify).

Before I start hacking the VFDModeSolver function, I thought that I would ask: do you already have an example of how to work with an anisotropic refractive index?

[lbolla]

There is an example of VFDModeSolver in https://github.com/lbolla/EMpy/blob/master/examples/ex_modesolver.py and https://github.com/lbolla/EMpy/blob/master/examples/ex_modesolver_2.py for the case of scalar eps, though. The idea of using 5 scalars rather than a 3x3 matrix for eps is that it would give the impression that eps_xz and other mixed z components are used, when they are not. It should be fairly easy to have an adapter function that takes a full matrix and converts it to the required tuple-based format. I am happy to accept PRs.

[DanHickstein]

I'm fine with using 5 scalars rather than a 3x3 matrix. I was just struggling with how, practically, to make a function (my eps_func) that returns a, say, 100x100 matrix of 5-element tuples. I could generate a 100x100x5 numpy array, but I do not have much experience mixing numpy arrays with tuples.

Perhaps this is easy to do...

[lbolla]

How about something like this?

import numpy

def epsfunc(x_, y_):
    eps = numpy.zeros((len(x_), len(y_), 5))
    for ix, xx in enumerate(x_):
        for iy, yy in enumerate(y_):
            if abs(xx - 1) <= 0.5 and abs(yy - 2) <= 1:
                eps[ix, iy, :] = [1, 2, 3, 4, 5]
            else:
                eps[ix, iy, :] = [5, 4, 3, 2, 1]
    return eps

x = numpy.linspace(0, 3, 10)
y = numpy.linspace(0, 3, 10)
eps = epsfunc(x, y)

print eps

It's not pretty, but it should work. I know understand your question and, indeed, an NxNx5 matrix is equivalent (in the sense that it work work just as well) as a NxN matrix of 5-tuples. In fact, that's what my code above does. The whole point is that each element eps[i, j, :] must be destructible into 5 elements.

[DanHickstein]

Yes, your example makes sense. Using a NxMx5 (anisotropic) or NxM (isotropic) matrix seems like the most intuitive way for the FD modesolver to operate. VFDModeSolver should simply check to see if it's being passed a 2D or 3D matrix.

However, this is currently not how it operates, and I don't think that the epsfunc is actually a valid input to VFDModeSolver (I tried swapping the epsfunc in ex_modesolver.py for the this one that you included above, and got an error). Using a 3D matrix causes an error during the numpy.c_ call on line 337, where it is basically assumed that the matrix is 2D:

 tmp = numpy.c_[tmp[:, 0:1], tmp, tmp[:, -1:]]

Also, I believe that an error would also occur on line 339, where the code explicitly checks if the elements of the tmp matrix are tuples:

if isinstance(tmp, tuple):
    epsxx, epsxy, epsyx, epsyy, epszz = tmp

[lbolla]

@DanHickstein Pull and try now. I've fixed the solver to accept an epsfunc returning an NxMx5 matrix. I've also updated the examples, adding an anisotropic eps and a semi-vectorial solver, as well.

[DanHickstein]

It seems to be working really well! Thanks!