;+ ; 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