;+
; NAME:
; sphericalfield
;
; PURPOSE:
; Calculates the complex electric field defined by an array of scattering
; coefficients.
;
; CATEGORY:
; Holography, light scattering, microscopy
;
; CALLING SEQUENCE:
; field = sphericalfield(x, y, z, a, lambda, $
; mpp = mpp)
;
; INPUTS:
; x: [npts] array of pixel coordinates [pixels]
; y: [npts] array of pixel coordinates [pixels]
; z: If field is required in a single plane, then
; z is the plane's distance from the sphere's center
; [pixels].
; Otherwise, z is an [npts] array of coordinates.
;
; NOTE: Ideally, x, y and z should be double precision.
; This is left to the calling program for efficiency.
;
; a: [2,nc] array of a and b scattering coefficients, where
; nc is the number of terms required for convergence.
;
; lambda: wavelength of light in medium [pixels]
;
; KEYWORD FLAGS:
; cartesian: If set, return field components in Cartesian
; coordinates. Default: Spherical polar coordinates
;
; OUTPUTS:
; field: [3,npts] complex values of field at the positions r.
; [0,*]: r component
; [1,*]: theta component
; [2,*]: phi component
;
; If CARTESIAN is set:
; [0,*]: x component (incident polarization)
; [1,*]: y component (transverse component)
; [2,*]: z component (axial component, relative to
; incident beam).
;
; REFERENCE:
; 1. Adapted from Chapter 4 in
; C. F. Bohren and D. R. Huffman,
; Absorption and Scattering of Light by Small Particles,
; (New York, Wiley, 1983).
; 2. W. J. Wiscombe,
; Improved Mie scattering algorithms,
; Applied Optics 19, 1505-1509 (1980).
;
; MODIFICATION HISTORY:
; Written by David G. Grier, New York University, 5/2007
; 6/9/2007: DGG finally read Section 4.8 in Bohren and Huffman about
; numerical stability of the recursions used to compute the scattering
; coefficients. Feh. Result is a total rewrite.
; 6/20/2007: DGG Calculate \tau_n(\cos\theta) and \pi_n(\cos\theta)
; according to recurrence relations in
; W. J. Wiscombe, Appl. Opt. 19, 1505-1509 (1980).
; This is supposed to improve numerical accuracy.
; 2/8/2008: DGG. Replaced single [3,npts] array of input coordinates
; with two [npts] arrays for x and y, and a separate input for z.
; Eliminated double() call for coordinates. Z may have 1 element or
; npts elements. Small documentation fixes.
; 4/3/2008: Bo Sun (Sephiroth), NYU: Calculate Lorenz-Mie a and b
; coefficients using continued fractions rather than recursion.
; Osman Akcakir from Arryx pointed out that the results are
; more accurate in extreme cases. Method described in
; William J. Lentz, "Generating Bessel functions in Mie scattering
; calculations using continued fractions," Appl. Opt. 15, 668-671
; (1976).
; 4/4/2008: DGG small code clean-ups and documentation. Added
; RECURSIVE keyword for backward compatibility in computing a and b
; coefficients.
; 4/11/2008: Sephiroth: Corrected small error in jump code for
; repeated fractions in Mie coefficients.
; 6/25/2008: DGG Don't clobber x coordinate input values.
; 10/9/2008: DGG adapted from SPHEREFIELD by separating out
; calculation of scattering coefficients, a_n and b_n. This
; is therefore more general, and can be replaced more
; readily with a GPU-accelerated version.
; 10/13/2008: DGG eliminated RECURSIVE keyword.
; 06/18/2010: DGG Added COMPILE_OPT.
; 09/04/2011: DGG Compute -\tau_n rather than \tau_n
;
; Copyright (c) 2007-2011, Bo Sun and David G. Grier
;
; UPDATES:
; The most recent version of this program may be obtained from
; http://physics.nyu.edu/grierlab/software.html
;
; LICENSE:
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License as
; published by the Free Software Foundation; either version 2 of the
; License, or (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
; General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
; 02111-1307 USA
;
; If the Internet and WWW are still functional when you are using
; this, you should be able to access the GPL here:
; http://www.gnu.org/copyleft/gpl.html
;-
function sphericalfield, x_, y_, z_, ab, lambda, $
cartesian = cartesian ; project to cartesian coordinates
COMPILE_OPT IDL2
npts = n_elements(x_)
nc = n_elements(ab[0,*])-1 ; number of terms required for convergence
k = 2.d * !dpi / lambda ; wavenumber in medium [pixel^-1]
ci = dcomplex(0,1)
; convert to spherical coordinates centered on the sphere.
; (r, theta, phi) is the spherical coordinate of the pixel
; at (x,y) in the imaging plane at distance z from the
; center of the sphere.
rho = sqrt(x_^2 + y_^2)
r = sqrt(rho^2 + z_^2)
theta = atan(rho, z_)
phi = atan(y_, x_)
costheta = cos(theta)
sintheta = sin(theta)
cosphi = cos(phi)
sinphi = sin(phi)
kr = k*r ; reduced radial coordinate
; starting points for recursive function evaluation ...
; ... Riccati-Bessel radial functions, page 478
sinkr = sin(kr)
coskr = cos(kr)
xi_nm2 = dcomplex(coskr, sinkr) ; \xi_{-1}(kr)
xi_nm1 = dcomplex(sinkr,-coskr) ; \xi_0(kr)
; ... angular functions (4.47), page 95
pi_nm1 = 0.d ; \pi_0(\cos\theta)
pi_n = 1.d ; \pi_1(\cos\theta)
; storage for vector spherical harmonics: [r,theta,phi]
Mo1n = dcomplexarr(3, npts)
Ne1n = dcomplexarr(3, npts)
; storage for scattered field
Es = dcomplexarr(3, npts)
; Compute field by summing multipole contributions
for n = 1.d, nc do begin
; upward recurrences ...
; ... Legendre factor (4.47)
; Method described by Wiscombe (1980)
swisc = pi_n * costheta
twisc = swisc - pi_nm1
tau_n = pi_nm1 - n * twisc ; -\tau_n(\cos\theta)
; ... Riccati-Bessel function, page 478
xi_n = (2.d*n - 1.d) * xi_nm1 / kr - xi_nm2 ; \xi_n(kr)
; vector spherical harmonics (4.50)
; Mo1n[0,*] = 0.d ; no radial component
Mo1n[1,*] = pi_n * xi_n ; ... divided by cosphi/kr
Mo1n[2,*] = tau_n * xi_n ; ... divided by sinphi/kr
dn = (n * xi_n)/kr - xi_nm1
Ne1n[0,*] = n*(n + 1.d) * pi_n * xi_n ; ... divided by cosphi sintheta/kr^2
Ne1n[1,*] = tau_n * dn ; ... divided by cosphi/kr
Ne1n[2,*] = pi_n * dn ; ... divided by sinphi/kr
; prefactor, page 93
En = ci^n * (2.d*n + 1.d) / n / (n + 1.d)
; the scattered field in spherical coordinates (4.45)
Es += En * (ci * ab[0,n] * Ne1n - ab[1,n] * Mo1n)
; upward recurrences ...
; ... angular functions (4.47)
; Method described by Wiscombe (1980)
pi_nm1 = pi_n
pi_n = swisc + (n + 1.d) * twisc / n
; ... Riccati-Bessel function
xi_nm2 = xi_nm1
xi_nm1 = xi_n
endfor
; geometric factors were divided out of the vector
; spherical harmonics for accuracy and efficiency ...
; ... put them back at the end.
Es[0,*] *= cosphi * sintheta / kr^2
Es[1,*] *= cosphi / kr
Es[2,*] *= sinphi / kr
; By default, the scattered wave is returned in spherical
; coordinates. Project components onto Cartesian coordinates.
; Assumes that the incident wave propagates along z and
; is linearly polarized along x
if keyword_set(cartesian) then begin
Ec = Es
Ec[0,*] = Es[0,*] * sintheta * cosphi
Ec[0,*] += Es[1,*] * costheta * cosphi
Ec[0,*] -= Es[2,*] * sinphi
Ec[1,*] = Es[0,*] * sintheta * sinphi
Ec[1,*] += Es[1,*] * costheta * sinphi
Ec[1,*] += Es[2,*] * cosphi
Ec[2,*] = Es[0,*] * costheta - Es[1,*] * sintheta
return, Ec
endif
return, Es
end