possible numerical problems
Closed this issue · 6 comments
dvincentwest commented
I am trying to use scattnlay, to reproduce figure 4.6 from Bohren & Huffman looking at the scattering properties of water droplets.
I pulled optical properties of water from Wikipedia:
https://en.wikipedia.org/wiki/Optical_properties_of_water_and_ice
There appear to be some numerical errors in calculating the Qext. I am wondering if this is common / expected or if it is a bug or if I am doing something wrong:
The following code was run in a Jupyter Notebook:
%pylab inline
import sys
sys.path.append('/home/vince/Coding/scattnlay/')
import scattnlay
wl_um = \
array([0.2 , 0.225, 0.25 , 0.275, 0.3 , 0.325, 0.35 , 0.375, 0.4 ,
0.425, 0.45 , 0.475, 0.5 , 0.525, 0.55 , 0.575, 0.6 , 0.625,
0.65 , 0.675, 0.7 , 0.725, 0.75 , 0.775, 0.8 , 0.825, 0.85 ,
0.875, 0.9 , 0.925, 0.95 , 0.975, 1. , 1.2 , 1.4 , 1.6 ,
1.8 , 2. , 2.2 , 2.4 , 2.6 , 2.65 , 2.7 , 2.75 , 2.8 ,
2.85 , 2.9 , 2.95 , 3. , 3.05 , 3.1 , 3.15 , 3.2 , 3.25 ,
3.3 , 3.35 , 3.4 , 3.45 , 3.5 , 3.6 , 3.7 , 3.8 , 3.9 ,
4. , 4.1 , 4.2 , 4.3 , 4.4 , 4.5 , 4.6 , 4.7 , 4.8 ,
4.9 , 5. , 5.1 , 5.2 , 5.3 , 5.4 , 5.5 , 5.6 , 5.7 ,
5.8 , 5.9 , 6. , 6.1 , 6.2 , 6.3 , 6.4 , 6.5 , 6.6 ,
6.7 , 6.8 , 6.9 , 7. ])
n_water = \
array([1.396, 1.373, 1.362, 1.354, 1.349, 1.346, 1.343, 1.341, 1.339,
1.338, 1.337, 1.336, 1.335, 1.334, 1.333, 1.333, 1.332, 1.332,
1.331, 1.331, 1.331, 1.33 , 1.33 , 1.33 , 1.329, 1.329, 1.329,
1.328, 1.328, 1.328, 1.327, 1.327, 1.327, 1.324, 1.321, 1.317,
1.312, 1.306, 1.296, 1.279, 1.242, 1.219, 1.188, 1.157, 1.142,
1.149, 1.201, 1.292, 1.371, 1.426, 1.467, 1.483, 1.478, 1.467,
1.45 , 1.432, 1.42 , 1.41 , 1.4 , 1.385, 1.374, 1.364, 1.357,
1.351, 1.346, 1.342, 1.338, 1.334, 1.332, 1.33 , 1.33 , 1.33 ,
1.328, 1.325, 1.322, 1.317, 1.312, 1.305, 1.298, 1.289, 1.277,
1.262, 1.248, 1.265, 1.319, 1.363, 1.357, 1.347, 1.339, 1.334,
1.329, 1.324, 1.321, 1.317])
k_water = \
array([1.10e-07, 4.90e-08, 3.35e-08, 2.35e-08, 1.60e-08, 1.08e-08,
6.50e-09, 3.50e-09, 1.86e-09, 1.30e-09, 1.02e-09, 9.35e-10,
1.00e-09, 1.32e-09, 1.96e-09, 3.60e-09, 1.09e-08, 1.39e-08,
1.64e-08, 2.23e-08, 3.35e-08, 9.15e-08, 1.56e-07, 1.48e-07,
1.25e-07, 1.82e-07, 2.93e-07, 3.91e-07, 4.86e-07, 1.06e-06,
2.93e-06, 3.48e-06, 2.89e-06, 9.89e-06, 1.38e-04, 8.55e-05,
1.15e-04, 1.10e-03, 2.89e-04, 9.56e-04, 3.17e-03, 6.70e-05,
1.90e-02, 5.90e-02, 1.15e-01, 1.85e-01, 2.68e-01, 2.98e-01,
2.72e-01, 2.40e-01, 1.92e-01, 1.35e-01, 9.24e-02, 6.10e-02,
3.68e-02, 2.61e-02, 1.95e-02, 1.32e-02, 9.40e-03, 5.15e-03,
3.60e-03, 3.40e-03, 3.80e-03, 4.60e-03, 5.62e-03, 6.88e-03,
8.45e-03, 1.03e-02, 1.34e-02, 1.47e-02, 1.57e-02, 1.50e-02,
1.37e-02, 1.24e-02, 1.11e-02, 1.01e-02, 9.80e-03, 1.03e-02,
1.16e-02, 1.42e-02, 2.03e-02, 3.30e-02, 6.22e-02, 1.07e-01,
1.31e-01, 8.80e-02, 5.70e-02, 4.49e-02, 3.92e-02, 3.56e-02,
3.37e-02, 3.27e-02, 3.22e-02, 3.20e-02])
xvals = 1 / wl_um
fig, (ax_n, ax_k) = subplots(2, 1, figsize=(6,6))
ax_n.plot(xvals, n_water)
ax_k.plot(xvals, k_water, 'r')
fig.tight_layout()
ax_k.set_xlabel('inverse microns')
ax_k.set_ylabel('k')
ax_n.set_ylabel('n')
ax_n.set_title('refractive index of water')
radius = 1.0 # um
N_medium = 1 # air
x_water = (2 * pi * N_medium * radius / wl_um).reshape(-1, 1)
m = ((n_water + 1j * k_water) / N_medium).reshape(-1, 1)
results = scattnlay.scattnlay(x_water, m)
Qext = results[1]
# create a finer grid for interpolation of refractive index properties
wl_fine = linspace(wl_um.min(), wl_um.max(), 1000)
n_water_fine = interp(wl_fine, wl_um, n_water)
k_water_fine = interp(wl_fine, wl_um, k_water)
x_water_fine = (2 * pi * N_medium * radius / wl_fine).reshape(-1, 1)
m_fine = ((n_water_fine + 1j * k_water_fine) / N_medium).reshape(-1, 1)
results_fine = scattnlay.scattnlay(x_water_fine, m_fine)
Qext_fine = results_fine[1]
fig, (ax1, ax2) = subplots(1, 2, figsize=(9, 4), dpi=90, facecolor='white')
ax1.plot(1/wl_fine, Qext_fine, label='interpolated')
ax1.plot(1/wl_um, Qext, label='table')
ax2.semilogy(1/wl_fine, Qext_fine, label='interpolated')
ax2.semilogy(1/wl_um, Qext, 'o', label='table')
for ax in (ax1, ax2):
ax.set_xlabel('inverse microns')
ax.set_ylabel('Excinction efficiency')
ax.legend()
fig.tight_layout()
idx = argmin(Qext)
ax1.text(1/wl_um[idx], Qext[idx],
f"neg Qext @ {wl_um[idx]} microns",
horizontalalignment='left')
ovidiopr commented
Dear Vince,
Thank you for the report. The calculations in the book by Bohren and
Huffman were done with a dielectric function that is reported in the same
book (Fig. 10.3). At first sight they look different than the ones that you
have used. Before going into the code, could you please check it?
Best regards,
Ovidio
…On Wed, Jun 26, 2019 at 5:29 PM Vince W. ***@***.***> wrote:
I am trying to use scattnlay, to reproduce figure 4.6 from Bohren &
Huffman
<https://www.wiley.com/en-us/Absorption+and+Scattering+of+Light+by+Small+Particles-p-9780471293408>
looking at the scattering properties of water droplets.
I pulled optical properties of water from Wikipedia:
https://en.wikipedia.org/wiki/Optical_properties_of_water_and_ice
There appear to be some numerical errors in calculating the Qext. I am
wondering if this is common / expected or if it is a bug or if I am doing
something wrong:
The following code was run in a Jupyter Notebook:
%pylab inline
import sys
sys.path.append('/home/vince/Coding/scattnlay/')import scattnlay
wl_um = \
array([0.2 , 0.225, 0.25 , 0.275, 0.3 , 0.325, 0.35 , 0.375, 0.4 ,
0.425, 0.45 , 0.475, 0.5 , 0.525, 0.55 , 0.575, 0.6 , 0.625,
0.65 , 0.675, 0.7 , 0.725, 0.75 , 0.775, 0.8 , 0.825, 0.85 ,
0.875, 0.9 , 0.925, 0.95 , 0.975, 1. , 1.2 , 1.4 , 1.6 ,
1.8 , 2. , 2.2 , 2.4 , 2.6 , 2.65 , 2.7 , 2.75 , 2.8 ,
2.85 , 2.9 , 2.95 , 3. , 3.05 , 3.1 , 3.15 , 3.2 , 3.25 ,
3.3 , 3.35 , 3.4 , 3.45 , 3.5 , 3.6 , 3.7 , 3.8 , 3.9 ,
4. , 4.1 , 4.2 , 4.3 , 4.4 , 4.5 , 4.6 , 4.7 , 4.8 ,
4.9 , 5. , 5.1 , 5.2 , 5.3 , 5.4 , 5.5 , 5.6 , 5.7 ,
5.8 , 5.9 , 6. , 6.1 , 6.2 , 6.3 , 6.4 , 6.5 , 6.6 ,
6.7 , 6.8 , 6.9 , 7. ])
n_water = \
array([1.396, 1.373, 1.362, 1.354, 1.349, 1.346, 1.343, 1.341, 1.339,
1.338, 1.337, 1.336, 1.335, 1.334, 1.333, 1.333, 1.332, 1.332,
1.331, 1.331, 1.331, 1.33 , 1.33 , 1.33 , 1.329, 1.329, 1.329,
1.328, 1.328, 1.328, 1.327, 1.327, 1.327, 1.324, 1.321, 1.317,
1.312, 1.306, 1.296, 1.279, 1.242, 1.219, 1.188, 1.157, 1.142,
1.149, 1.201, 1.292, 1.371, 1.426, 1.467, 1.483, 1.478, 1.467,
1.45 , 1.432, 1.42 , 1.41 , 1.4 , 1.385, 1.374, 1.364, 1.357,
1.351, 1.346, 1.342, 1.338, 1.334, 1.332, 1.33 , 1.33 , 1.33 ,
1.328, 1.325, 1.322, 1.317, 1.312, 1.305, 1.298, 1.289, 1.277,
1.262, 1.248, 1.265, 1.319, 1.363, 1.357, 1.347, 1.339, 1.334,
1.329, 1.324, 1.321, 1.317])
k_water = \
array([1.10e-07, 4.90e-08, 3.35e-08, 2.35e-08, 1.60e-08, 1.08e-08,
6.50e-09, 3.50e-09, 1.86e-09, 1.30e-09, 1.02e-09, 9.35e-10,
1.00e-09, 1.32e-09, 1.96e-09, 3.60e-09, 1.09e-08, 1.39e-08,
1.64e-08, 2.23e-08, 3.35e-08, 9.15e-08, 1.56e-07, 1.48e-07,
1.25e-07, 1.82e-07, 2.93e-07, 3.91e-07, 4.86e-07, 1.06e-06,
2.93e-06, 3.48e-06, 2.89e-06, 9.89e-06, 1.38e-04, 8.55e-05,
1.15e-04, 1.10e-03, 2.89e-04, 9.56e-04, 3.17e-03, 6.70e-05,
1.90e-02, 5.90e-02, 1.15e-01, 1.85e-01, 2.68e-01, 2.98e-01,
2.72e-01, 2.40e-01, 1.92e-01, 1.35e-01, 9.24e-02, 6.10e-02,
3.68e-02, 2.61e-02, 1.95e-02, 1.32e-02, 9.40e-03, 5.15e-03,
3.60e-03, 3.40e-03, 3.80e-03, 4.60e-03, 5.62e-03, 6.88e-03,
8.45e-03, 1.03e-02, 1.34e-02, 1.47e-02, 1.57e-02, 1.50e-02,
1.37e-02, 1.24e-02, 1.11e-02, 1.01e-02, 9.80e-03, 1.03e-02,
1.16e-02, 1.42e-02, 2.03e-02, 3.30e-02, 6.22e-02, 1.07e-01,
1.31e-01, 8.80e-02, 5.70e-02, 4.49e-02, 3.92e-02, 3.56e-02,
3.37e-02, 3.27e-02, 3.22e-02, 3.20e-02])
xvals = 1 / wl_um
fig, (ax_n, ax_k) = subplots(2, 1, figsize=(6,6))
ax_n.plot(xvals, n_water)
ax_k.plot(xvals, k_water, 'r')
fig.tight_layout()
ax_k.set_xlabel('inverse microns')
ax_k.set_ylabel('k')
ax_n.set_ylabel('n')
ax_n.set_title('refractive index of water')
[image: image]
<https://user-images.githubusercontent.com/9872343/60193264-f6391680-97fc-11e9-8e14-963ae6dba3a9.png>
radius = 1.0 # um
N_medium = 1 # air
result_keys = "terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2".split(", ")
idx = result_keys.index("Qext")
# create a finer grid for interpolation of refractive index properties
wl_fine = linspace(wl_um.min(), wl_um.max(), 1000)
n_water_fine = interp(wl_fine, wl_um, n_water)
k_water_fine = interp(wl_fine, wl_um, k_water)
results = list()for wl, n, k in zip(wl_um, n_water, k_water):
N_particle = n + 1j * k
m = N_particle / N_medium
x = 2 * pi * N_medium * radius / wl
x = np.array(x).reshape(1, 1)
m = np.array(m).reshape(1)
results.append(scattnlay.scattnlay(x, m))
Qext = [vals[idx] for vals in results]
results_fine = list()for wl, n, k in zip(wl_fine, n_water_fine, k_water_fine):
N_particle = n + 1j * k
m = N_particle / N_medium
x = 2 * pi * N_medium * radius / wl
x = np.array(x).reshape(1, 1)
m = np.array(m).reshape(1)
results_fine.append(scattnlay.scattnlay(x, m))
Qext_fine = [vals[idx] for vals in results_fine]
fig, (ax1, ax2) = subplots(1, 2, figsize=(9, 4), dpi=90, facecolor='white')
ax1.plot(1/wl_fine, Qext_fine, label='interpolated')
ax1.plot(1/wl_um, Qext, label='table')
ax2.semilogy(1/wl_fine, Qext_fine, label='interpolated')
ax2.semilogy(1/wl_um, Qext, 'o', label='table')for ax in (ax1, ax2):
ax.set_xlabel('inverse microns')
ax.set_ylabel('Excinction efficiency')
ax.legend()
fig.tight_layout()
idx = argmin(Qext)
ax1.text(1/wl_um[idx], Qext[idx],
f"neg Qext @ {wl_um[idx]} microns",
horizontalalignment='left')
[image: image]
<https://user-images.githubusercontent.com/9872343/60193336-17016c00-97fd-11e9-900d-aa72eeb197b5.png>
—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
<#17?email_source=notifications&email_token=ACNOR2O3QDOLWBFVS6SVQG3P4ODO3A5CNFSM4H3TP5X2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4G325K7Q>,
or mute the thread
<https://github.com/notifications/unsubscribe-auth/ACNOR2OZJIBMXMISFI2B6LTP4ODO3ANCNFSM4H3TP5XQ>
.
dvincentwest commented
Dear Ovidio, I will post the further code below, but it may be that there are instabilities near certain multiples of pi. See this quick test code snippet:
x_test = array([[5 * pi]]) # large problems for radius = 1.0um and wavelength = 0.4um where x = 5 * pi
m_test = array([[1.4]])
print(scattnlay.scattnlay(x_test, m_test)) # Qext is the -42.85114173 entry belowprints:
(array([43]),
array([-42.85114173]), # <-- Qext entry
array([378.14033685]),
array([-420.99147859]),
array([130.11363197]),
array([-361.44011008]),
array([0.84251517]),
array([-8.82451019]),
array([], shape=(1, 0), dtype=complex128),
array([], shape=(1, 0), dtype=complex128))
dvincentwest commented
Dear Ovidio, please see the updated code below. OCR was a bit fun for this, as the n/k values were extracted from the original paper referenced in B&H:
Optical Constants of Water in the 200-nm to 200-μm Wavelength Region
George M. Hale and Marvin R. Querry
Applied Optics Vol. 12, Issue 3, pp. 555-563 (1973)
https://doi.org/10.1364/AO.12.000555
wl_um = \
array([ 0.2 , 0.225, 0.25 , 0.275, 0.3 , 0.325, 0.35 ,
0.375, 0.4 , 0.425, 0.45 , 0.475, 0.5 , 0.525,
0.55 , 0.575, 0.6 , 0.625, 0.65 , 0.675, 0.7 ,
0.725, 0.75 , 0.775, 0.8 , 0.825, 0.85 , 0.875,
0.9 , 0.925, 0.95 , 0.975, 1. , 1.2 , 1.4 ,
1.6 , 1.8 , 2. , 2.2 , 2.4 , 2.6 , 2.65 ,
2.7 , 2.75 , 2.8 , 2.85 , 2.9 , 2.95 , 3. ,
3.05 , 3.1 , 3.15 , 3.2 , 3.25 , 3.3 , 3.35 ,
3.4 , 3.45 , 3.5 , 3.6 , 3.7 , 3.8 , 3.9 ,
4. , 4.1 , 4.2 , 4.3 , 4.4 , 4.5 , 4.6 ,
4.7 , 4.8 , 4.9 , 5. , 5.1 , 5.2 , 5.3 ,
5.4 , 5.5 , 5.6 , 5.7 , 5.8 , 5.9 , 6. ,
6.1 , 6.2 , 6.3 , 6.4 , 6.5 , 6.6 , 6.7 ,
6.8 , 6.9 , 7. , 7.1 , 7.2 , 7.3 , 7.4 ,
7.5 , 7.6 , 7.7 , 7.8 , 7.9 , 8. , 8.2 ,
8.4 , 8.6 , 8.8 , 9. , 9.2 , 9.4 , 9.6 ,
9.8 , 10. , 10.5 , 11. , 11.5 , 12. , 12.5 ,
13. , 13.5 , 14. , 14.5 , 15. , 15.5 , 16. ,
16.5 , 17. , 17.5 , 18. , 18.5 , 19. , 19.5 ,
20. , 21. , 22. , 23. , 24. , 25. , 26. ,
27. , 28. , 29. , 30. , 32. , 34. , 36. ,
38. , 40. , 42. , 44. , 46. , 48. , 50. ,
60. , 70. , 80. , 90. , 100. , 110. , 120. ,
130. , 140. , 150. , 160. , 170. , 180. , 190. ,
200. ])
n_water = \
array([1.396, 1.373, 1.362, 1.354, 1.349, 1.346, 1.343, 1.341, 1.339,
1.338, 1.337, 1.336, 1.335, 1.334, 1.333, 1.333, 1.332, 1.332,
1.331, 1.331, 1.331, 1.33 , 1.33 , 1.33 , 1.329, 1.329, 1.329,
1.328, 1.328, 1.328, 1.327, 1.327, 1.327, 1.324, 1.321, 1.317,
1.312, 1.306, 1.296, 1.279, 1.242, 1.219, 1.188, 1.157, 1.142,
1.149, 1.201, 1.292, 1.371, 1.426, 1.467, 1.483, 1.478, 1.467,
1.45 , 1.432, 1.42 , 1.41 , 1.4 , 1.385, 1.374, 1.364, 1.357,
1.351, 1.346, 1.342, 1.338, 1.334, 1.332, 1.33 , 1.33 , 1.33 ,
1.328, 1.325, 1.322, 1.317, 1.312, 1.305, 1.298, 1.289, 1.277,
1.262, 1.248, 1.265, 1.319, 1.363, 1.357, 1.347, 1.339, 1.334,
1.329, 1.324, 1.321, 1.317, 1.314, 1.312, 1.309, 1.307, 1.304,
1.302, 1.299, 1.297, 1.294, 1.291, 1.286, 1.281, 1.275, 1.269,
1.262, 1.255, 1.247, 1.239, 1.229, 1.218, 1.185, 1.153, 1.126,
1.111, 1.123, 1.146, 1.177, 1.21 , 1.241, 1.27 , 1.297, 1.325,
1.351, 1.376, 1.401, 1.423, 1.443, 1.461, 1.476, 1.48 , 1.487,
1.5 , 1.511, 1.521, 1.531, 1.539, 1.545, 1.549, 1.551, 1.551,
1.546, 1.536, 1.527, 1.522, 1.519, 1.522, 1.53 , 1.541, 1.555,
1.587, 1.703, 1.821, 1.886, 1.924, 1.957, 1.966, 2.004, 2.036,
2.056, 2.069, 2.081, 2.094, 2.107, 2.119, 2.13 ])
k_water = \
array([1.10e-07, 4.90e-08, 3.35e-08, 2.35e-08, 1.60e-08, 1.08e-08,
6.50e-09, 3.50e-09, 1.86e-09, 1.30e-09, 1.02e-09, 9.35e-10,
1.00e-09, 1.32e-09, 1.96e-09, 3.60e-09, 1.09e-08, 1.39e-08,
1.64e-08, 2.23e-08, 3.35e-08, 9.15e-08, 1.56e-07, 1.48e-07,
1.25e-07, 1.82e-07, 2.93e-07, 3.91e-07, 4.86e-07, 1.06e-06,
2.93e-06, 3.48e-06, 2.89e-06, 9.89e-06, 1.38e-04, 8.55e-05,
1.15e-04, 1.10e-03, 2.89e-04, 9.56e-04, 3.17e-03, 6.70e-03,
1.90e-02, 5.90e-02, 1.15e-01, 1.85e-01, 2.68e-01, 2.98e-01,
2.72e-01, 2.40e-01, 1.92e-01, 1.35e-01, 9.24e-02, 6.10e-02,
3.68e-02, 2.61e-02, 1.95e-02, 1.32e-02, 9.40e-03, 5.15e-03,
3.60e-03, 3.40e-03, 3.80e-03, 4.60e-03, 5.62e-03, 6.88e-03,
8.45e-03, 1.03e-02, 1.34e-02, 1.47e-02, 1.57e-02, 1.50e-02,
1.37e-02, 1.24e-02, 1.11e-02, 1.01e-02, 9.80e-03, 1.03e-02,
1.16e-02, 1.42e-02, 2.03e-02, 3.30e-02, 6.22e-02, 1.07e-01,
1.31e-01, 8.80e-02, 5.70e-02, 4.49e-02, 3.92e-02, 3.56e-02,
3.37e-02, 3.27e-02, 3.22e-02, 3.20e-02, 3.20e-02, 3.21e-02,
3.22e-02, 3.24e-02, 3.26e-02, 3.28e-02, 3.31e-02, 3.35e-02,
3.39e-02, 3.43e-02, 3.51e-02, 3.61e-02, 3.72e-02, 3.85e-02,
3.99e-02, 4.15e-02, 4.33e-02, 4.54e-02, 4.79e-02, 5.08e-02,
6.62e-02, 9.68e-02, 1.42e-01, 1.99e-01, 2.59e-01, 3.05e-01,
3.43e-01, 3.70e-01, 3.88e-01, 4.02e-01, 4.14e-01, 4.22e-01,
4.28e-01, 4.29e-01, 4.29e-01, 4.26e-01, 4.21e-01, 4.14e-01,
4.04e-01, 3.93e-01, 3.82e-01, 3.73e-01, 3.67e-01, 3.61e-01,
3.56e-01, 3.50e-01, 3.44e-01, 3.38e-01, 3.33e-01, 3.28e-01,
3.24e-01, 3.29e-01, 3.43e-01, 3.61e-01, 3.85e-01, 4.09e-01,
4.36e-01, 4.62e-01, 4.88e-01, 5.14e-01, 5.87e-01, 5.76e-01,
5.47e-01, 5.36e-01, 5.32e-01, 5.31e-01, 5.26e-01, 5.14e-01,
5.00e-01, 4.95e-01, 4.96e-01, 4.97e-01, 4.99e-01, 5.01e-01,
5.04e-01])
xvals = wl_um
fig, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, figsize=(10, 5), facecolor='white')
ax1.semilogx(wl_um, n_water)
ax3.semilogx(wl_um, k_water, 'r')
ax2.plot(1/wl_um, n_water)
ax4.plot(1/wl_um, k_water, 'r')
ax1.set_ylabel('n')
ax3.set_ylabel('k')
ax3.set_xlabel('wavelength (microns)')
ax4.set_xlabel('inverse microns')
fig.tight_layout()
radius = 1.0 # um
N_medium = 1 # air
x_water = (2 * pi * N_medium * radius / wl_um).reshape(-1, 1)
m = ((n_water + 1j * k_water) / N_medium).reshape(-1, 1)
results = scattnlay.scattnlay(x_water, m)
Qext = results[1]
# create a finer grid for interpolation of refractive index properties
wl_fine = linspace(wl_um.min(), wl_um.max(), 1000)
n_water_fine = interp(wl_fine, wl_um, n_water)
k_water_fine = interp(wl_fine, wl_um, k_water)
x_water_fine = (2 * pi * N_medium * radius / wl_fine).reshape(-1, 1)
m_fine = ((n_water_fine + 1j * k_water_fine) / N_medium).reshape(-1, 1)
results_fine = scattnlay.scattnlay(x_water_fine, m_fine)
Qext_fine = results_fine[1]
fig, (ax1, ax2) = subplots(1, 2, figsize=(9, 4), dpi=90, facecolor='white')
ax1.plot(1/wl_fine, Qext_fine, label='interpolated')
ax1.plot(1/wl_um, Qext, label='table')
ax2.semilogy(1/wl_fine, Qext_fine, label='interpolated')
ax2.semilogy(1/wl_um, Qext, 'o', label='table')
for ax in (ax1, ax2):
ax.set_xlabel('inverse microns')
ax.set_ylabel('Excinction efficiency')
ax.legend()
fig.tight_layout()
idx = argmin(Qext)
ax1.text(1/wl_um[idx], Qext[idx],
f"neg Qext @ {wl_um[idx]} microns",
horizontalalignment='left')
kostyfisik commented
The particle is large, so you need to use mp feature. Note that data starts from 0.2 mkm, so you cannot go above 5 mkm^-1. Same for *.yml data file from refractive index info. My result looks to be OK, full example https://drive.google.com/open?id=1GZ_RcFVPw2Oy0TGb17BUt6lVwuFJgH82.
dvincentwest commented
@kostyfisik thanks very much. A lot of work to realize my mistake was just needing to pass a single boolean value, but I am very appreciative of your time.
kostyfisik commented
@dvincentwest You are welcome :)