add large sphere example

examples/Opt.Commun.-2010-Geints-Nanojets_of_dielectric_microspheres/fig1.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 -*-
+#
+#    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  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 at least one of the following references:
+#    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
+#        a multilayered sphere," Computer Physics Communications,
+#        vol. 180, Nov. 2009, pp. 2348-2354.
+#    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
+#        calculation of electromagnetic near-field for a multilayered
+#        sphere," Computer Physics Communications, vol. 214, May 2017,
+#        pp. 225-230.
+#
+#    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 scattnlay
+from scattnlay import fieldnlay
+from scattnlay import scattnlay
+import numpy as np
+from matplotlib import pyplot as plt
+
+npts = 151
+factor = 3. # plot extent compared to sphere radius
+index_H2O = 1.33+0.j
+
+WL = 0.532 #mkm
+total_r = 3 #mkm
+
+nm = 1.0 # host medium
+x = 2.0 * np.pi * np.array([total_r / 4.0 * 3.0, total_r], dtype=np.float64) / WL
+m = np.array((index_H2O, index_H2O), dtype=np.complex128) / nm
+
+print("x =", x)
+print("m =", m)
+terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+    np.array([x]), np.array([m]))
+print("Qsca = " + str(Qsca))
+terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+    np.array([x]), np.array([m]), mp=True)
+print("mp Qsca = " + str(Qsca))
+
+scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
+zero = np.zeros(npts*npts, dtype = np.float64)
+
+coordX, coordZ = np.meshgrid(scan, scan)
+coordX.resize(npts * npts)
+coordZ.resize(npts * npts)
+coordY = zero
+
+terms, E, H = fieldnlay(
+    np.array([x]), np.array([m]),
+    coordX, coordY, coordZ,
+    mp=True
+)
+Ec = E[0, :, :]
+Er = np.absolute(Ec)
+Eabs2 = (Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
+Eabs_data = np.resize(Eabs2, (npts, npts))
+label = r'$|E|^2$'
+plt.imshow(Eabs_data,
+           cmap='jet',
+           vmin=0., vmax=14
+
+           )
+print(np.min(Eabs_data), np.max(Eabs_data))
+plt.savefig("R"+str(total_r)+"mkm_mp.jpg")
+# plt.show()

examples/field-Ag-flow.py

@@ -34,25 +34,26 @@ from fieldplot import fieldplot
 
 import numpy as np
 import cmath
+
 # # a)
-#WL=400 #nm
-#core_r = WL/20.0
-#epsilon_Ag = -2.0 + 10.0j
+# WL=400 #nm
+# core_r = WL/20.0
+# epsilon_Ag = -2.0 + 10.0j
 
 # # b)
-#WL=400 #nm
-#core_r = WL/20.0
-#epsilon_Ag = -2.0 + 1.0j
+# WL=400 #nm
+# core_r = WL/20.0
+# epsilon_Ag = -2.0 + 1.0j
 
 # c)
-WL=354 #nm
-core_r = WL/20.0
+WL = 354  # nm
+core_r = WL / 20.0
 epsilon_Ag = -2.0 + 0.28j
 
 # d)
-#WL=367 #nm
-#core_r = WL/20.0
-#epsilon_Ag = -2.71 + 0.25j
+# WL=367 #nm
+# core_r = WL/20.0
+# epsilon_Ag = -2.71 + 0.25j
 
 # WL=500 #nm
 # core_r = 615.0
@@ -64,36 +65,37 @@ index_Ag = np.sqrt(epsilon_Ag)
 # n2 = 0.565838 + 7.23262j
 nm = 1.0
 
-x = 2.0*np.pi*np.array([core_r/4.0*3.0, core_r], dtype = np.float64)/WL
+x = 2.0 * np.pi * np.array([core_r / 4.0 * 3.0, core_r], dtype=np.float64) / WL
 
-m = np.array((index_Ag, index_Ag), dtype = np.complex128)/nm
+m = np.array((index_Ag, index_Ag), dtype=np.complex128) / nm
 
-print( "x =", x)
-print( "m =", m)
+print("x =", x)
+print("m =", m)
 
-comment='bulk-WL'+str(WL)+'nm_r'+str(core_r)+'nm_epsilon'+str(epsilon_Ag)+'-flow'
-WL_units='nm'
+comment = 'bulk-WL' + str(WL) + 'nm_r' + str(core_r) + 'nm_epsilon' + str(epsilon_Ag) + '-flow'
+WL_units = 'nm'
 npts = 251
-factor=2.1
+factor = 2.1
 flow_total = 9
-#flow_total = 21
-#flow_total = 0
-crossplane='XZ'
-#crossplane='YZ'
+# flow_total = 21
+# flow_total = 0
+crossplane = 'XZ'
+# crossplane='YZ'
 # crossplane='XY'
 
 # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
-field_to_plot='Pabs'
-#field_to_plot='angleEx'
+field_to_plot = 'Pabs'
+# field_to_plot='angleEx'
 
 
 import matplotlib.pyplot as plt
-fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+
+fig, axs = plt.subplots(1, 1)  # , sharey=True, sharex=True)
 fig.tight_layout()
-fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
-          subplot_label=' ',is_flow_extend=False)
+fieldplot(fig, axs, x, m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+          subplot_label=' ', is_flow_extend=False)
 
-#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+# fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
 
 # for ax in axs:
 #     ax.locator_params(axis='x',nbins=5)
@@ -101,8 +103,8 @@ fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts,
 
 fig.subplots_adjust(hspace=0.3, wspace=-0.1)
 
-plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
-                    +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
+plt.savefig(comment + "-R" + str(int(round(x[-1] * WL / 2.0 / np.pi))) + "-" + crossplane + "-"
+            + field_to_plot + ".pdf", pad_inches=0.02, bbox_inches='tight')
 
 plt.draw()