Merge branch 'master' of github.com:ovidiopr/scattnlay

examples/H2O/H2O-Hale.yml

@@ -0,0 +1,182 @@
+# this file is part of refractiveindex.info database
+# refractiveindex.info database is in the public domain
+# copyright and related rights waived via CC0 1.0
+
+REFERENCES: "G. M. Hale and M. R. Querry. Optical constants of water in the 200-nm to 200-µm wavelength region, <a href=\"https://doi.org/10.1364/AO.12.000555\"><i>Appl. Opt.</i> <b>12</b>, 555-563 (1973)</a><br>See also <a href=\"https://refractiveindex.info/?shelf=main&book=H2O&page=Segelstein\"><i>Segelstein 1981</i></a> for more recent data from the same group."
+COMMENTS: "Liquid water (H<sub>2</sub>O) at 25 °C"
+DATA:
+  - type: tabulated nk
+    data: |
+        0.200 1.396 1.10E-7
+        0.225 1.373 4.90E-8
+        0.250 1.362 3.35E-8
+        0.275 1.354 2.35E-8
+        0.300 1.349 1.60E-8
+        0.325 1.346 1.08E-8
+        0.350 1.343 6.50E-9
+        0.375 1.341 3.50E-9
+        0.400 1.339 1.86E-9
+        0.425 1.338 1.30E-9
+        0.450 1.337 1.02E-9
+        0.475 1.336 9.35E-10
+        0.500 1.335 1.00E-9
+        0.525 1.334 1.32E-9
+        0.550 1.333 1.96E-9
+        0.575 1.333 3.60E-9
+        0.600 1.332 1.09E-8
+        0.625 1.332 1.39E-8
+        0.650 1.331 1.64E-8
+        0.675 1.331 2.23E-8
+        0.700 1.331 3.35E-8
+        0.725 1.330 9.15E-8
+        0.750 1.330 1.56E-7
+        0.775 1.330 1.48E-7
+        0.800 1.329 1.25E-7
+        0.825 1.329 1.82E-7
+        0.850 1.329 2.93E-7
+        0.875 1.328 3.91E-7
+        0.900 1.328 4.86E-7
+        0.925 1.328 1.06E-6
+        0.950 1.327 2.93E-6
+        0.975 1.327 3.48E-6
+        1.0 1.327 2.89E-6
+        1.2 1.324 9.89E-6
+        1.4 1.321 1.38E-4
+        1.6 1.317 8.55E-5
+        1.8 1.312 1.15E-4
+        2.0 1.306 1.10E-3
+        2.2 1.296 2.89E-4
+        2.4 1.279 9.56E-4
+        2.6 1.242 3.17E-3
+        2.65 1.219 6.70E-3
+        2.70 1.188 0.019
+        2.75 1.157 0.059
+        2.80 1.142 0.115
+        2.85 1.149 0.185
+        2.90 1.201 0.268
+        2.95 1.292 0.298
+        3.00 1.371 0.272
+        3.05 1.426 0.240
+        3.10 1.467 0.192
+        3.15 1.483 0.135
+        3.20 1.478 0.0924
+        3.25 1.467 0.0610
+        3.30 1.450 0.0368
+        3.35 1.432 0.0261
+        3.40 1.420 0.0195
+        3.45 1.410 0.0132
+        3.50 1.400 0.0094
+        3.6 1.385 0.00515
+        3.7 1.374 0.00360
+        3.8 1.364 0.00340
+        3.9 1.357 0.00380
+        4.0 1.351 0.00460
+        4.1 1.346 0.00562
+        4.2 1.342 0.00688
+        4.3 1.338 0.00845
+        4.4 1.334 0.0103
+        4.5 1.332 0.0134
+        4.6 1.330 0.0147
+        4.7 1.330 0.0157
+        4.8 1.330 0.0150
+        4.9 1.328 0.0137
+        5.0 1.325 0.0124
+        5.1 1.322 0.0111
+        5.2 1.317 0.0101
+        5.3 1.312 0.0098
+        5.4 1.305 0.0103
+        5.5 1.298 0.0116
+        5.6 1.289 0.0142
+        5.7 1.277 0.0203
+        5.8 1.262 0.0330
+        5.9 1.248 0.0622
+        6.0 1.265 0.107
+        6.1 1.319 0.131
+        6.2 1.363 0.0880
+        6.3 1.357 0.0570
+        6.4 1.347 0.0449
+        6.5 1.339 0.0392
+        6.6 1.334 0.0356
+        6.7 1.329 0.0337
+        6.8 1.324 0.0327
+        6.9 1.321 0.0322
+        7.0 1.317 0.0320
+        7.1 1.314 0.0320
+        7.2 1.312 0.0321
+        7.3 1.309 0.0322
+        7.4 1.307 0.0324
+        7.5 1.304 0.0326
+        7.6 1.302 0.0328
+        7.7 1.299 0.0331
+        7.8 1.297 0.0335
+        7.9 1.294 0.0339
+        8.0 1.291 0.0343
+        8.2 1.286 0.0351
+        8.4 1.281 0.0361
+        8.6 1.275 0.0372
+        8.8 1.269 0.0385
+        9.0 1.262 0.0399
+        9.2 1.255 0.0415
+        9.4 1.247 0.0433
+        9.6 1.239 0.0454
+        9.8 1.229 0.0479
+        10.0 1.218 0.0508
+        10.5 1.185 0.0662
+        11.0 1.153 0.0968
+        11.5 1.126 0.142
+        12.0 1.111 0.199
+        12.5 1.123 0.259
+        13.0 1.146 0.305
+        13.5 1.177 0.343
+        14.0 1.210 0.370
+        14.5 1.241 0.388
+        15.0 1.270 0.402
+        15.5 1.297 0.414
+        16.0 1.325 0.422
+        16.5 1.351 0.428
+        17.0 1.376 0.429
+        17.5 1.401 0.429
+        18.0 1.423 0.426
+        18.5 1.443 0.421
+        19.0 1.461 0.414
+        19.5 1.476 0.404
+        20.0 1.480 0.393
+        21.0 1.487 0.382
+        22 1.500 0.373
+        23 1.511 0.367
+        24 1.521 0.361
+        25 1.531 0.356
+        26 1.539 0.350
+        27 1.545 0.344
+        28 1.549 0.338
+        29 1.551 0.333
+        30 1.551 0.328
+        32 1.546 0.324
+        34 1.536 0.329
+        36 1.527 0.343
+        38 1.522 0.361
+        40 1.519 0.385
+        42 1.522 0.409
+        44 1.530 0.436
+        46 1.541 0.462
+        48 1.555 0.488
+        50 1.587 0.514
+        60 1.703 0.587
+        70 1.821 0.576
+        80 1.886 0.547
+        90 1.924 0.536
+        100 1.957 0.532
+        110 1.966 0.531
+        120 2.004 0.526
+        130 2.036 0.514
+        140 2.056 0.500
+        150 2.069 0.495
+        160 2.081 0.496
+        170 2.094 0.497
+        180 2.107 0.499
+        190 2.119 0.501
+        200 2.130 0.504
+SPECS:
+    n_absolute: true
+    wavelength_vacuum: true
+    temperature: 25 °C

examples/H2O/calc-spectra.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2019  Konstantin Ladutenko <kostyfisik@gmail.com>
+#
+#    This file is part of python-scattnlay
+#
+#    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 3 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.
+#
+#    The only additional remark is that we expect that all publications
+#    describing work using this software, or all commercial products
+#    using it, cite the following reference:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+import sys
+sys.path.insert(0,'../..')  # to be able to import scattnlay from the upper dir
+
+from scattnlay import scattnlay,scattcoeffs,fieldnlay
+
+import matplotlib.pyplot as plt
+import numpy as np
+import cmath
+
+from optical_constants import read_refractive_index_from_yaml as get_index
+
+from_rWL = 0.01
+to_rWL = 5  # limit from H2O-Hale.yml data
+step_rWL = 0.01
+rWLs = np.arange(from_rWL, to_rWL+step_rWL/2., step_rWL);
+WLs = 1/rWLs #mkm
+
+index_H2O = get_index('H2O-Hale.yml', WLs, "mkm")
+
+print(index_H2O)
+
+x = np.ones((1), dtype = np.float64)
+m = np.ones((1), dtype = np.complex128)
+
+
+core_r = 1 #mkm
+
+
+Qext_vec = []
+
+for i in range(len(WLs)):
+    WL = WLs[i]
+    x[0] = 2.0*np.pi*core_r/WL#/4.0*3.0
+    m[0] = index_H2O[:,1][i]
+    
+    terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+        np.array(x), np.array(m),
+        mp=True
+        # mp=False
+    )
+    print(np.array([Qext]))
+    Qext_vec.append(Qext)
+
+fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+axs.plot(rWLs, Qext_vec, color="black")
+plt.ylim(0, 4.3) 
+axs.set_xlabel("$1/\lambda, \mu m^{-1}$")
+axs.set_ylabel("$Q_{ext}$")
+plt.title("Scattnlay")
+plt.savefig("spectra.pdf",pad_inches=0.02, bbox_inches='tight')
+plt.show()
+plt.clf()
+plt.close()

examples/H2O/optical_constants.py

@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+# GPL license from Smuthi project
+"""Provide functionality to read optical constants in format provided by `refractiveindex.info <https://refractiveindex.info/>`_ website"""
+from scipy.interpolate import interp1d
+import io
+import numpy as np
+import yaml
+
+
+def read_refractive_index_from_yaml(filename, vacuum_wavelength, units="mkm", kind=1):
+    """Read optical constants in format provided by refractiveindex.info website.
+
+    Args:
+            filename (str): path and file name for yaml data
+                            downloaded from refractiveindex.info
+            vacuum_wavelength (float or np.array): wavelengths where refractive
+                                                   index data is needed
+            units (str): units for wavelength. currently, microns ('mkm' or 'um')
+                         and nanometers ('nm') can be selected
+            kind (int): order of interpolation
+    
+    Returns:
+            A pair (or np.array of pairs) of wavelength and
+            corresponding refractive index (complex)
+    """
+    if units == "nm": 
+        factor = 1000
+    elif units in ("mkm", "um"):
+        factor = 1
+    else:
+        raise NotImplementedError("Converting wavelength into '"+units
+                                  +"' units for refractive index data"
+                                  +" was not implemented.")
+    
+    the_file = yaml.load(open(filename))['DATA'][0]
+    data_type = the_file['type']
+    if data_type != 'tabulated nk':
+        raise NotImplementedError("Input data type '"+data_type
+                                  +"' available in file "+filename
+                                  +" was not implemented.")
+    data = the_file['data'].splitlines()
+    data_split = []
+    for wl in data:
+        data_split.append(wl.split())
+    data_num = []
+    for wl in data_split:
+        record = []
+        for val in wl:
+            record.append(float(val))
+        data_num.append(record)
+    data_np = np.array(data_num)
+    data_wl = data_np[:,0]*factor
+    eps_re = data_np[:,1]
+    eps_im = data_np[:,2]
+    f_re = interp1d(data_wl, eps_re, kind=kind)
+    f_im = interp1d(data_wl, eps_im, kind=kind)
+    data_out = np.transpose(np.vstack((vacuum_wavelength, f_re(vacuum_wavelength)+f_im(vacuum_wavelength)*1j)))
+    if len(data_out) == 1:
+        return data_out[0]
+    return data_out