-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathdbr_spherefield.pro
179 lines (166 loc) · 6.91 KB
/
dbr_spherefield.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
;+
; NAME:
; dbr_spherefield
;
; PURPOSE:
; Calculates the complex electric field scattered by a sphere illuminated
; by a plane wave linearly polarized in the x direction.
;
; CATEGORY:
; Holography, light scattering, microscopy
;
; CALLING SEQUENCE:
; field = dbr_spherefield(x, y, z, a, np, $
; nm = nm, lambda = lambda, $
; mpp = mpp, k = k)
;
; 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: radius of sphere [micrometers]
;
; np : (complex) refractive index of sphere
;
; KEYWORD PARAMETERS:
; nm: (complex) refractive index of medium.
; Default: 1.33 (water)
;
; lambda: vacuum wavelength of light [micrometers]
; Default: 0.632 (HeNe)
;
; mpp: Microns per pixel.
; Default: 0.135 (Nikon rig)
;
; precision: precision with which field is calculated.
; Example: precision=1.e-6 allows errors smaller than
; one part in a million.
; Default: Machine precision (effectively)
;
; KEYWORD FLAGS:
; cartesian: If set, return field components in Cartesian
; coordinates. Default: spherical polar coordinates
;
; gpu: If set, use GPU_SPHERICALFIELD to
; compute series of vector spherical harmonics.
; This can substantially improve performance on systems
; with GPUlib installed.
;
; 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).
;
; PROCEDURE:
; Calls SPHERE_COEFFCIENTS to obtain scattering
; coefficients, and then calls SPHERICALFIELD to sum up
; the series of vector spherical harmonics.
;
; 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).
;
; MODIFICATION HISTORY:
; Written by David G. Grier, New York University, 5/2007
; 06/09/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.
; 06/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.
; 02/08/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.
; 04/03/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).
; 04/04/2008: DGG small code clean-ups and documentation. Added
; RECURSIVE keyword for backward compatibility in computing a and b
; coefficients.
; 04/11/2008: Sephiroth: Corrected small error in jump code for
; repeated fractions in Mie coefficients.
; 06/25/2008: DGG Don't clobber x coordinate input values.
; 10/09/2008: DGG converted to a shell for calling SPHERE_COEFFICIENTS
; and SPHERICALFIELD. This greatly simplifies generalizing the
; code to other systems. Added PRECISION keyword.
; 10/12/2008: DGG Added GPU keyword.
; 01/15/2009: DGG check for PRECISION with KEYWORD_SET to avoid
; unnecessary processing when PRECISION=0. Documentation cleanups.
; 06/18/2010: DGG Added COMPILE_OPT.
; 03/19/2014: Adapted from spherefield.pro by David Ruffner,
; New York University
;
; Copyright (c) 2006-2010 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 shold be able to access the GPL here:
; http://www.gnu.org/copyleft/gpl.html
;-
function dbr_spherefield, x, y, z, a, np, $
nm = nm, $ ; refractive index of medium
lambda = lambda, $ ; vacuum wavelength of light
mpp = mpp, $ ; micrometers per pixel
precision = precision, $ ; accuracy of scattering coefficients
k = k, $ ; scaled wavenumber in medium
cartesian = cartesian, $ ; project to cartesian coordinates
gpu = gpu ; use GPU acceleration
COMPILE_OPT IDL2
if n_elements(nm) eq 0 then $
nm = dcomplex(1.3326, 1.5e-8) ; refractive index of medium (water)
if n_elements(lambda) eq 0 then $
lambda = 0.6328d ; wavelength in vacuum [microns]
if n_elements(mpp) eq 0 then $
mpp = 0.101d ; length calibration [microns/pixel]
ab = sphere_coefficients(a, np, nm, lambda)
if n_elements(precision) eq 1 then begin
fac = total(abs(ab[1,*]),1)
w = where(fac gt precision*max(fac))
ab = ab[*,[0,w]] ; retain first coefficient for bookkeeping
endif
lambda_m = lambda / real_part(nm) / mpp ; medium wavelength [pixel]
k = 2.d * !pi / lambda_m ; wavenumber in medium [pixel^{-1}]
field = dbr_sphericalfield(x, y, z, ab, lambda_m, cartesian=cartesian)
return, field
end