|
|
@@ -0,0 +1,121 @@
|
|
|
+#!/usr/bin/env python
|
|
|
+# -*- coding: UTF-8 -*-
|
|
|
+#
|
|
|
+# Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.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/>.
|
|
|
+
|
|
|
+# This test case calculates the electric field in the
|
|
|
+# E-k plane, for an spherical Si-Ag-Si nanoparticle. Core radius is 17.74 nm,
|
|
|
+# inner layer 23.31nm, outer layer 22.95nm. Working wavelength is 800nm, we use
|
|
|
+# silicon epsilon=13.64+i0.047, silver epsilon= -28.05+i1.525
|
|
|
+
|
|
|
+import os, cmath
|
|
|
+import numpy as np
|
|
|
+from scattnlay import fieldnlay
|
|
|
+from fieldplot import fieldplot
|
|
|
+
|
|
|
+
|
|
|
+if __name__ == '__main__':
|
|
|
+ import argparse
|
|
|
+
|
|
|
+ parser = argparse.ArgumentParser()
|
|
|
+
|
|
|
+ parser.add_argument("dirnames", nargs='*', default='.', help="read all data from DIR(S)")
|
|
|
+ parser.add_argument("-f", "--filename", dest="fname", nargs='?', default=None,
|
|
|
+ help="name of 'n' file")
|
|
|
+ parser.add_argument("-w", "--wavelength", dest="wl", default=3.75, type=float,
|
|
|
+ help="wavelength of electromagnetic wave")
|
|
|
+ parser.add_argument("-r", "--radius", dest="rad", default=None, type=float,
|
|
|
+ help="radius of PEC sphere")
|
|
|
+ parser.add_argument("-t", "--thickness", dest="tc", default=0.8, type=float,
|
|
|
+ help="thickness of cloaking layer")
|
|
|
+ parser.add_argument("-n", "--npoints", dest="npts", default=101, type=int,
|
|
|
+ help="number of points for the grid")
|
|
|
+
|
|
|
+ args = parser.parse_args()
|
|
|
+
|
|
|
+ for dirname in args.dirnames:
|
|
|
+ print "Calculating spectra for data file(s) in dir '%s'..." % (dirname)
|
|
|
+
|
|
|
+ wl = args.wl # cm
|
|
|
+ if (args.rad is None):
|
|
|
+ Rs = 0.75*wl # cm
|
|
|
+ else:
|
|
|
+ Rs = args.rad # cm
|
|
|
+ tc = args.tc # cm
|
|
|
+
|
|
|
+ if (args.fname is None):
|
|
|
+ files = [x for x in os.listdir('%s/' % (dirname)) if x.endswith('.dat')]
|
|
|
+ files.sort()
|
|
|
+ else:
|
|
|
+ files = [args.fname]
|
|
|
+
|
|
|
+ npts = args.npts # cm
|
|
|
+
|
|
|
+ if not os.path.exists('%s/flow-results/' % (dirname)):
|
|
|
+ os.makedirs('%s/flow-results/' % (dirname))
|
|
|
+
|
|
|
+ Rt = Rs + tc # cm
|
|
|
+
|
|
|
+ print "Wl = %.2f, Rs = %.2f, tc = %.2f, Rt = %.2f" % (wl, Rs, tc, Rt)
|
|
|
+
|
|
|
+ ms = 1.0 + 40.0j
|
|
|
+ for i, fname in enumerate(files):
|
|
|
+ print "Calculating spectra for file '%s'..." % (fname)
|
|
|
+
|
|
|
+ basename = os.path.splitext(fname)[0]
|
|
|
+
|
|
|
+ nvalues = np.loadtxt('%s/%s' % (dirname, fname))*1.0 + 1e-11j
|
|
|
+
|
|
|
+ tl = tc/len(nvalues) # cm
|
|
|
+ r = [Rs]
|
|
|
+
|
|
|
+ for i in range(len(nvalues)):
|
|
|
+ r += [r[i] + tl]
|
|
|
+
|
|
|
+ x = np.ones((1, len(nvalues) + 1), dtype = np.float64)
|
|
|
+ m = np.ones((1, len(nvalues) + 1), dtype = np.complex128)
|
|
|
+
|
|
|
+ x[0] = 2.0*np.pi*np.array(r, dtype = np.float64)/wl
|
|
|
+ m[0] = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
|
|
|
+
|
|
|
+ factor = 2
|
|
|
+ comment='PEC-'+basename
|
|
|
+ WL_units='cm'
|
|
|
+ flow_total = 39
|
|
|
+ #flow_total = 0
|
|
|
+ #crossplane='XZ'
|
|
|
+ crossplane='YZ'
|
|
|
+ #crossplane='XY'
|
|
|
+
|
|
|
+ # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
|
|
|
+ field_to_plot='Pabs'
|
|
|
+ #field_to_plot='Eabs'
|
|
|
+ #field_to_plot='angleEx'
|
|
|
+ fieldplot(x[0],m[0], wl, comment, WL_units, crossplane, field_to_plot, npts,
|
|
|
+ factor, flow_total, pl=0, outline_width=0.1)
|
|
|
+
|
|
|
+
|
|
|
+ print "Done!!"
|
|
|
+
|
