SHARE
TWEET

Untitled

a guest Dec 23rd, 2017 110 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ;+
  2. ; NAME:
  3. ;    sphere_coefficients
  4. ;
  5. ; PURPOSE:
  6. ;    Calculates the Mie scattering coefficients
  7. ;    for a multilayered sphere illuminated
  8. ;    by a coherent plane wave linearly polarized in the x direction.
  9. ;
  10. ; CATEGORY:
  11. ;    Holography, light scattering, microscopy
  12. ;
  13. ; CALLING SEQUENCE:
  14. ;    ab = sphere_coefficients(ap, np, nm, lambda)
  15. ;
  16. ; INPUTS:
  17. ;    ap: [nlayers] radii of layered sphere [micrometers]
  18. ;        NOTE: ap and np are reordered automatically so that
  19. ;        ap is in ascending order.
  20. ;
  21. ;    np: [nlayers] (complex) refractive indexes of sphere's layers
  22. ;
  23. ;    nm: (complex) refractive index of medium
  24. ;
  25. ;    lambda: wavelength of light [micrometers]
  26. ;
  27. ; OUTPUTS:
  28. ;    ab: [2,nc] complex a and b scattering coefficients.
  29. ;
  30. ; REFERENCES:
  31. ; 1. Adapted from Chapter 8 in
  32. ;    C. F. Bohren and D. R. Huffman,
  33. ;    Absorption and Scattering of Light by Small Particles,
  34. ;    (New York, Wiley, 1983).
  35. ;
  36. ; 2. W. Yang,
  37. ;    Improved recursive algorithm for
  38. ;    light scattering by a multilayered sphere,
  39. ;    Applied Optics 42, 1710--1720 (2003).
  40. ;
  41. ; 3. O. Pena, U. Pal,
  42. ;    Scattering of electromagnetic radiation by a multilayered sphere,
  43. ;    Computer Physics Communications 180, 2348--2354 (2009).
  44. ;    NB: Equation numbering follows this reference.
  45. ;
  46. ; 4. W. J. Wiscombe,
  47. ;    Improved Mie scattering algorithms,
  48. ;    Applied Optics 19, 1505-1509 (1980).
  49. ;
  50. ; EXAMPLE:
  51. ;  ab = sphere_coefficients([1.0, 1.5], [1.45, 1.56], 1.3326, 0.6328)
  52. ;
  53. ; MODIFICATION HISTORY:
  54. ; Replaces sphere_coefficients.pro, which calculated scattering
  55. ;    coefficients for a homogeneous sphere.  This earlier version
  56. ;    was written by Bo Sun and David G. Grier, and was based on
  57. ;    algorithms by W. J. Wiscombe (1980) and W. J. Lentz (1976).
  58. ;
  59. ; 10/31/2010 Written by F. C. Cheong, New York University
  60. ; 11/02/2010 David G. Grier (NYU) Formatting.  Explicit casts
  61. ;    to double precision.  Use complex functions.
  62. ; 04/10/2011 DGG Cleaned up Nstop code.  Use array math rather than
  63. ;    iteration where convenient.  Eliminate an and bn subroutines.
  64. ;    Simplify function names.  Use IDL 8 array notation.  Fixed
  65. ;    bug in bn coefficients for multilayer spheres.
  66. ; 05/31/2011 DGG Fixed corner condition in Nstop code.
  67. ;
  68. ; Copyright (c) 2010-2011, F. C. Cheong and D. G. Grier
  69. ;
  70.  
  71. ;;;;;
  72. ;
  73. ; Nstop
  74. ;
  75. ; Number of terms to keep in the partial wave expansion
  76. ;
  77. function Nstop, x, m
  78.  
  79. COMPILE_OPT IDL2, HIDDEN
  80.  
  81. ;;; Wiscombe (1980)
  82. xl = x[-1]
  83. if xl lt 8. then $
  84.   ns = floor(xl + 4. * xl^(1./3.) + 1.) $
  85. else if xl lt 4200. then $
  86.   ns = floor(xl + 4.05 * xl^(1./3.) + 2.) $
  87. else if xl gt 4199. then $
  88.   ns = floor(xl + 4. * xl^(1./3.) + 2.)
  89.  
  90. ;;; Yang (2003) Eq. (30)
  91. return, double(floor(max([ns, abs(x*m), abs(shift(x, -1)*m)])) + 15)
  92.  
  93. end
  94.  
  95. ;;;;;
  96. ;
  97. ; sphere_coefficients
  98. ;
  99. function sphere_coefficients, ap, np, nm, lambda
  100.  
  101. COMPILE_OPT IDL2
  102.  
  103. nlayers = n_elements(ap)
  104.  
  105. if n_elements(np) ne nlayers then $
  106.   message, "Error Warning: ap and np must have the same number of elements"
  107.  
  108. ; arrange shells in size order
  109. if nlayers gt 1 then begin
  110.   order = sort(ap)
  111.   ap = ap[order]
  112.   np = np[order]
  113. endif
  114.  
  115. x = abs(2.d * !dpi * nm * ap / lambda) ; size parameter [array]
  116. m = dcomplex(np/nm)                    ; relative refractive index [array]
  117.  
  118. nmax = Nstop(x, m)              ; number of terms in partial-wave expansion
  119.  
  120. ci = dcomplex(0, 1)             ; imaginary unit
  121.  
  122. ; arrays for storing results
  123. ab = dcomplexarr(2, nmax+1, /nozero)
  124.  
  125. D1     = dcomplexarr(nmax+2)
  126. D1_a   = dcomplexarr(nlayers, nmax+2)
  127. D1_am1 = dcomplexarr(nlayers, nmax+2)
  128.  
  129. D3     = dcomplexarr(nmax+1)
  130. D3_a   = dcomplexarr(nlayers, nmax+1)
  131. D3_am1 = dcomplexarr(nlayers, nmax+1)
  132.  
  133. Psi         = dcomplexarr(nmax+1)
  134. Zeta        = dcomplexarr(nmax+1)
  135. PsiZeta     = dcomplexarr(nmax+1)
  136. PsiZeta_a   = dcomplexarr(nlayers, nmax+1)
  137. PsiZeta_am1 = dcomplexarr(nlayers, nmax+1)
  138.  
  139. Q  = dcomplexarr(nlayers, nmax+1)
  140. Ha = dcomplexarr(nlayers, nmax+1)
  141. Hb = dcomplexarr(nlayers, nmax+1)
  142.  
  143. ; Calculate D1, D3 and PsiZeta for Z1 in the first layer
  144. z1 = x[0] * m[0]
  145. ; D1_a[0, nmax + 1] = dcomplex(0) ; Eq. (16a)
  146. for n = nmax + 1.d, 1.d, -1.d do $ ; downward recurrence Eq. (16b)
  147.   D1_a[0, n-1] = n/z1 - 1.d/(D1_a[0, n] + n/z1)
  148.  
  149. PsiZeta_a[0, 0] = 0.5d * (1.d - exp(2.d * ci * z1)) ; Eq. (18a)
  150. D3_a[0, 0] = ci                                     ; Eq. (18a)
  151. for n = 1.d, nmax do begin      ;upward recurrence Eq. (18b)
  152.   PsiZeta_a[0, n] = PsiZeta_a[0, n-1] * $
  153.                     (n/z1 - D1_a[0, n-1]) * (n/z1 - D3_a[0, n-1])
  154.   D3_a[0, n] = D1_a[0, n] + ci/PsiZeta_a[0, n]
  155. endfor
  156.  
  157. ; Ha and Hb in the core
  158. Ha[0, *] = D1_a[0, 0:-2]     ; Eq. (7a)
  159. Hb[0, *] = D1_a[0, 0:-2]     ; Eq. (8a)
  160.  
  161. ; Iterate from layer 2 to layer L
  162. for ii = 1, nlayers - 1 do begin
  163.   z1 = x[ii] * m[ii]
  164.   z2 = x[ii-1] * m[ii]
  165.   ; Downward recurrence for D1, Eqs. (16a) and (16b)
  166. ;   D1_a[ii, nmax+1]   = dcomplex(0)      ; Eq. (16a)
  167. ;   D1_am1[ii, nmax+1] = dcomplex(0)
  168.   for n = nmax + 1.d, 1.d, -1.d do begin ; Eq. (16 b)
  169.      D1_a[ii, n-1]   = n/z1 - 1.d/(D1_a[ii, n]   + n/z1)
  170.      D1_am1[ii, n-1] = n/z2 - 1.d/(D1_am1[ii, n] + n/z2)
  171.   endfor
  172.  
  173.   ; Upward recurrence for PsiZeta and D3, Eqs. (18a) and (18b)
  174.   PsiZeta_a[ii, 0]   = 0.5d * (1.d - exp(2.d * ci * z1)) ; Eq. (18a)
  175.   PsiZeta_am1[ii, 0] = 0.5d * (1.d - exp(2.d * ci * z2))
  176.   D3_a[ii, 0]   = ci
  177.   D3_am1[ii, 0] = ci
  178.   for n = 1.d, nmax do begin   ; Eq. (18b)
  179.      PsiZeta_a[ii, n]   = PsiZeta_a[ii, n-1] * $
  180.                           (n/z1 -  D1_a[ii, n-1]) * $
  181.                           (n/z1 -  D3_a[ii, n-1])
  182.      PsiZeta_am1[ii, n] = PsiZeta_am1[ii, n-1] * $
  183.                           (n/z2 - D1_am1[ii, n-1]) * $
  184.                           (n/z2 - D3_am1[ii, n-1])
  185.      D3_a[ii, n]   = D1_a[ii, n]   + ci/PsiZeta_a[ii, n]
  186.      D3_am1[ii, n] = D1_am1[ii, n] + ci/PsiZeta_am1[ii, n]
  187.   endfor
  188.  
  189.   ; Upward recurrence for Q
  190.   Q[ii, 0] = (exp(-2.d * ci * z2) - 1.d) / (exp(-2.d * ci * z1) - 1.d)
  191.   for n = 1.d, nmax do begin
  192.      Num = (z1 * D1_a[ii, n]   + n) * (n - z1 * D3_a[ii, n-1])
  193.      Den = (z2 * D1_am1[ii, n] + n) * (n - z2 * D3_am1[ii, n-1])
  194.      Q[ii, n] = (x[ii-1]/x[ii])^2 * Q[ii, n-1] * Num/Den
  195.   endfor
  196.  
  197.   ; Upward recurrence for Ha and Hb, Eqs. (7b), (8b) and (12) - (15)
  198.   for n = 1.d, nmax do begin
  199.      G1 = m[ii] * Ha[ii-1, n] - m[ii-1] * D1_am1[ii, n]
  200.      G2 = m[ii] * Ha[ii-1, n] - m[ii-1] * D3_am1[ii, n]
  201.      Temp = Q[ii, n] * G1
  202.      Num = G2 * D1_a[ii, n] - Temp * D3_a[ii, n]
  203.      Den = G2 - Temp
  204.      Ha[ii, n] = Num/Den
  205.  
  206.      G1 = m[ii-1] * Hb[ii-1, n] - m[ii] * D1_am1[ii, n]
  207.      G2 = m[ii-1] * Hb[ii-1, n] - m[ii] * D3_am1[ii, n]
  208.      Temp = Q[ii, n] * G1
  209.      Num = G2 * D1_a[ii, n] - Temp * D3_a[ii, n]
  210.      Den = G2 - Temp
  211.      Hb[ii, n] = Num/Den
  212.   endfor
  213. endfor                          ;ii (layers)
  214.  
  215. z1 = dcomplex(x[-1])
  216. ; Downward recurrence for D1, Eqs. (16a) and (16b)
  217. ; D1[nmax+1] = dcomplex(0)   ; Eq. (16a)
  218. for n = nmax, 1.d, -1.d do $ ; Eq. (16b)
  219.   D1[n-1] = n/z1 - (1.d/(D1[n] + n/z1))
  220.  
  221. ; Upward recurrence for Psi, Zeta, PsiZeta and D3, Eqs. (18a) and (18b)
  222. Psi[0]     = sin(z1)       ; Eq. (18a)
  223. Zeta[0]    = -ci * exp(ci * z1)
  224. PsiZeta[0] = 0.5d * (1.d - exp(2.d * ci * z1))
  225. D3[0] = ci
  226. for n = 1.d, nmax do begin ; Eq. (18b)
  227.   Psi[n]  = Psi[n-1]  * (n/z1 - D1[n-1])
  228.   Zeta[n] = Zeta[n-1] * (n/z1 - D3[n-1])
  229.   PsiZeta[n] = PsiZeta[n-1] * (n/z1 -D1[n-1]) * (n/z1 - D3[n-1])
  230.   D3[n] = D1[n] + ci/PsiZeta[n]
  231. endfor
  232.  
  233. ; Scattering coefficients, Eqs. (5) and (6)
  234. n = dindgen(nmax + 1)
  235. ab[0, *]  = (Ha[-1, *]/m[-1] + n/x[-1]) * Psi  - shift(Psi,  1) ; Eq. (5)
  236. ab[0, *] /= (Ha[-1, *]/m[-1] + n/x[-1]) * Zeta - shift(Zeta, 1)
  237. ab[1, *]  = (Hb[-1, *]*m[-1] + n/x[-1]) * Psi  - shift(Psi,  1) ; Eq. (6)
  238. ab[1, *] /= (Hb[-1, *]*m[-1] + n/x[-1]) * Zeta - shift(Zeta, 1)
  239. ab[*, 0]  = dcomplex(0)
  240.  
  241. return, ab
  242. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top