add comment

Makefile

@@ -28,10 +28,10 @@ source:
 cython: scattnlay.pyx
     # create c++ code for double precision module
 	$(CYTHON) --cplus scattnlay.pyx -o $(SRCDIR)/scattnlay.cpp
-	# # create c++ code for MP module
-	# ln -s scattnlay.pyx scattnlay_mp.pyx
-	# $(CYTHON) --cplus scattnlay_mp.pyx -o $(SRCDIR)/scattnlay_mp.cpp
-	# rm scattnlay_mp.pyx
+	# create c++ code for MP module
+	ln -s scattnlay.pyx scattnlay_mp.pyx
+	$(CYTHON) --cplus scattnlay_mp.pyx -o $(SRCDIR)/scattnlay_mp.cpp
+	rm scattnlay_mp.pyx
 
 python_ext: $(SRCDIR)/nmie.cc $(SRCDIR)/py_nmie.cc $(SRCDIR)/scattnlay.cpp $(SRCDIR)/scattnlay_mp.cpp
 	$(PYTHON) setup.py build_ext --inplace

examples/calc-expansion-coeffs-spectra.py

@@ -25,53 +25,65 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from scattnlay import fieldnlay, scattnlay, expansioncoeffs, scattcoeffs
+# This script reproduces Fig.2a from "Nonradiating anapole modes in
+# dielectric nanoparticles" by Miroshnichenko, Andrey E. et al in  
+# Nature Communications DOI:10.1038/ncomms9069
 
-import numpy as np
-import cmath
 
-WL=550 #nm
-core_r = 102
-#core_r = 120
+from scattnlay import fieldnlay, scattnlay, expansioncoeffs, scattcoeffs
+import cmath
+import matplotlib.pyplot as plt
+import numpy as np
 
-index_NP = 4.0
 x = np.ones((1), dtype = np.float64)
 m = np.ones((1), dtype = np.complex128)
 
-import matplotlib.pyplot as plt
+WL=550 #nm
+core_r = 180
+index_NP = 4.0
+
+from_R = 120/2.0
+to_R = 240/2.0
 
-# dipole - 1 -- 2
-# quad   - 2 -- 4
-# octo   - 3 -- 8
-# hex    - 4 -- 16
-# 32     - 5 -- 32
 npts = 151
 ext = ".png"
-# npts = 351
+npts = 351
 # ext = ".pdf"
-x[0] = 2.0*np.pi*core_r/WL#/4.0*3.0
-m[0] = index_NP
-#    for mode_type, field_to_plot, WL, mode_n,  crossplane, isStream in plot_params :
-comment='bulk-NP-R'+str(core_r)+'nm-WL'+str(WL)
-# scattnlay function process several NP designs in one call
-x = np.array([x])
-m = np.array([m])
-#print(x.shape)
-terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
-nmax = terms
-#print(nmax)
-terms, an, bn = scattcoeffs(x, m, nmax)
-#print(terms, an)
-terms, aln, bln, cln, dln = expansioncoeffs(x, m, nmax)
-print(terms, dln)
-# fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
-# fig.tight_layout()
-# fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+comment='bulk-NP-WL'+str(WL)
+
+val_all = []
+
+all_R = np.linspace(from_R, to_R, npts)
+for core_r in all_R:
+    x[0] = 2.0*np.pi*core_r/WL
+    m[0] = index_NP
+    terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(np.array([x]), np.array([m]))
+    terms, an, bn = scattcoeffs(np.array([x]), np.array([m]),terms)
 
-# # plt.savefig("Egor3/"+"%02d"%(i)+
-# #             comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
-# #                 +field_to_plot+"-mode"+mode+mt+st+ext,pad_inches=0.02, bbox_inches='tight')
-# plt.clf()
-# plt.close()
+    terms, aln, bln, cln, dln = expansioncoeffs(np.array([x]), np.array([m]), terms)
+    order = 0
+    val_all.append([#Qsca,
+                    np.abs(aln[0][-1][order]),np.abs(dln[0][0][order])
+        #,np.abs(an[0][order])
+    ])
+#print( )
+val_all = np.array(val_all)
+#print()
+#print(terms, np.abs(dln))
+fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+fig.tight_layout()
+fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+all_R = all_R*2
+plt.plot(all_R, val_all[:,0],label="electric dipole, scatt.")
+plt.plot(all_R, val_all[:,1],label="electric dipole, internal")
+plt.legend()
+#plt.plot(all_R, val_all[:,2])
+plt.ylim(0,4)
+plt.xlabel("D, nm")
+plt.ylabel("Mie coefficients")
+plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+ext,
+            pad_inches=0.02, bbox_inches='tight')
+plt.clf()
+plt.close()
 
 #print("end")