Mie scattering with scattnlay

[paulmueller]

This probably yields faster results than bare Mie simulations and thus might become useful at some point.

Implementation might look like this:

    if distance < radius:
        propdist = radius
    else:
        propdist = distance

    if isinstance(extent, _types_indexable):
        extent = np.array(extent)
    else:
        extent = np.array((extent,extent))
    
    n1 = n_object
    nm = n_medium

    twopi = 2*np.pi*nm

    x = np.ones((1, 1), dtype = np.float64)
    x[0, 0] = twopi*radius
    
    m = np.ones((1, 1), dtype = np.complex128)
    m[0, 0] = n1/nm
    
    nptsx = extent[0]*resolution
    nptsy = extent[1]*resolution

    assert nptsx == np.int(nptsx), "grid size must compute to integer"
    assert nptsy == np.int(nptsy), "grid size must compute to integer"
    nptsx = int(nptsx)
    nptsy = int(nptsy)

    # If you are wondering about the 2PIs, see:
    # https://github.com/ovidiopr/scattnlay/issues/6
    
    scanx = (np.linspace(-extent[0]/2, extent[0]/2, nptsx, endpoint=True) + offset_x)*twopi
    scany = (np.linspace(-extent[1]/2, extent[1]/2, nptsy, endpoint=True) + offset_y)*twopi
    
    coordX, coordY = np.meshgrid(scany, scanx)
    coordX.resize(nptsx*nptsy)
    coordY.resize(nptsx*nptsy)
    coordZ = np.ones(nptsx*nptsy, dtype=np.float64)*distance*twopi
    
    coord = np.vstack((coordY, coordX, coordZ)).transpose()
    
    _terms, E, _H = scattnlay.fieldnlay(x, m, coord)
    
    Ex = E[:,:,0].reshape(nptsx, nptsy)
    background = np.exp(1j*2*np.pi*propdist*n_medium)
    
    field = Ex/background

    # numerically refocus if required
    if distance < radius:
        field = nrefocus.refocus(field,
                                 d=-(radius-distance)*resolution,
                                 nm=n_medium,
                                 res=resolution)

    # coordinates:
    # excerpt from gmm_field documentation:
    #    Here, we are using a cubic grid with 101 x 101 x 101 = 10**6
    #    points that extends from -20 nm to +20 nm in each dimension.
    if isinstance(extent, (list, tuple, np.ndarray)):
        x = np.linspace(-extent[0]/2, extent[0]/2, int(np.round(resolution*extent[0])), endpoint=True)
        y = np.linspace(-extent[1]/2, extent[1]/2, int(np.round(resolution*extent[1])), endpoint=True)
        cx, cy = np.meshgrid(x,y)
    else:    
        x = np.linspace(-extent/2, extent/2, int(np.round(resolution*extent)))
        cy, cx = np.meshgrid(x,x)

    cx += offset_x
    cy += offset_y

    return field, cx, cy