From 8881f9c2a5a1aac9a40c5de65da536174e4439a2 Mon Sep 17 00:00:00 2001 From: yurkin Date: Mon, 30 Sep 2013 06:56:15 +0000 Subject: [PATCH] Implemented calculations of decay-rate enhancement for point-dipole incident beam. Total part is calculated with the help of a new function DecayCross() in crosssec.c/h (similar to Cext), non-radiative part - by normalizing Cabs. In case of surface, all parts of decay-rate enhancement (total, rad, nonrad) are normalized to that in free space, but additionally enhancement due to surface alone is shown. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Still, it is not clear whether the non-radiative part should be adjusted to account for absorption in the substrate (when it is absorbing). Also, incident polarizations (Y and X) are not shown in the log for point-dipole incident beam. This is a bit problematic for standard LDR polarizability, since its factor S depends on both prop and incPol. But it probably makes little sense to use LDR polarizability for such incident field anyway (recommendation should be added to the manual). One test was run for dipole near a sphere was run and matched exactly to first row of Table 1 in S. D'Agostino, F. Della Sala, and L.C. Andreani, “Dipole-excited surface plasmons in metallic nanoparticles: Engineering decay dynamics within the discrete-dipole approximation,” Phys. Rev. B 87, 205413 (2013). However, a few hacks were used to exactly match the simulation conditions (based on raw data provided by Stefania D'Agostino. --- src/CalculateE.c | 18 ++++++++++++++++++ src/GenerateB.c | 28 ++++++++++++++++++++++++++++ src/crosssec.c | 23 +++++++++++++++++++++++ src/crosssec.h | 1 + src/param.c | 21 ++++++++++++++------- 5 files changed, 84 insertions(+), 7 deletions(-) diff --git a/src/CalculateE.c b/src/CalculateE.c index 99ddbe16..cd7c0e06 100644 --- a/src/CalculateE.c +++ b/src/CalculateE.c @@ -46,6 +46,8 @@ extern double * restrict muel_alpha; // defined and initialized in crosssec.c extern const Parms_1D phi_sg; extern const double ezLab[3],exSP[3]; +// defined and initialized in GenerateB.c +extern const double C0dipole,C0dipole_refl; // defined and initialized in param.c extern const bool store_int_field,store_dip_pol,store_beam,store_scat_grid,calc_Cext,calc_Cabs, calc_Csca,calc_vec,calc_asym,calc_mat_force,store_force,store_ampl; @@ -716,6 +718,22 @@ static void CalcIntegralScatQuantities(const enum incpol which) CCfile=FOpenErr(fname_cs,"w",ONE_POS); if (calc_Cext) PrintBoth(CCfile,"Cext\t= "GFORM"\nQext\t= "GFORM"\n",Cext,Cext*inv_G); if (calc_Cabs) PrintBoth(CCfile,"Cabs\t= "GFORM"\nQabs\t= "GFORM"\n",Cabs,Cabs*inv_G); + if (beamtype==B_DIPOLE) { + double self=1; + if (surface) self+=C0dipole_refl/C0dipole; + double tot=self+DecayCross()/C0dipole; + fprintf(CCfile,"\nDecay-rate enhancement\n\n"); + printf("\nDecay-rate enhancement:\n"); + PrintBoth(CCfile,"Total\t= "GFORM"\n",tot); + if (calc_Cabs) { // for simplicity we keep a single condition here + double nonrad=Cabs/C0dipole; + double rad=tot-nonrad; + // TODO: !!! how nonrad should account for absorption inside the surface??? + PrintBoth(CCfile,"Radiat\t= "GFORM"\n",rad); + PrintBoth(CCfile,"Nonrad\t= "GFORM"\n",nonrad); + } + if (surface) PrintBoth(CCfile,"Surface\t= "GFORM"\n",self); + } if (all_dir) fprintf(CCfile,"\nIntegration\n\n"); if (calc_Csca) { Csca=ScaCross(f_suf); diff --git a/src/GenerateB.c b/src/GenerateB.c index 16aa6c0b..daae0fd7 100644 --- a/src/GenerateB.c +++ b/src/GenerateB.c @@ -46,6 +46,9 @@ extern const char *beam_fnameY; extern const char *beam_fnameX; extern const opt_index opt_beam; +// used in CalculateE.c +double C0dipole,C0dipole_refl; // inherent cross sections of exciting dipole (in free space and addition due to surface) + // used in crosssec.c double beam_center_0[3]; // position of the beam center in laboratory reference frame /* complex wave amplitudes of secondary waves (with phase relative to particle center); @@ -327,6 +330,31 @@ void GenerateB (const enum incpol which, // x - or y polarized incident light cvAdd(v1,b+j,b+j); } } + /* calculate dipole inherent cross sections (for decay rate enhancements) + * General formula is C0=4pi*k*Im(p0*.G(r0,r0).p0) and it is used for reflection; but for direct + * interaction it is expected to be singular, so we use an explicit formula for point dipole: + * C0=(8pi/3)*k^4*|p0|^2. Thus we discard the choice of "-int ...", but still the used value should be + * correct for most formulations, e.g. poi, fcd, fcd_st, igt_so. Moreover, it is also logical, since the + * exciting dipole is really point one, in contrast to the dipoles composing the particle. + */ + C0dipole=2*FOUR_PI_OVER_THREE*WaveNum*WaveNum*WaveNum*WaveNum; + printf("C0=%g\n",C0dipole); + if (surface) { + r1[0]=r1[1]=0; + r1[2]=2*(beam_center[2]+hsub); + (*ReflTerm_real)(r1,gt); + double tmp; + /* the following expression uses that prop is real and a specific (anti-)symmetry of the gt + * a general expression is commented out below + */ + tmp=prop[0]*prop[0]*cimag(gt[0])+2*prop[0]*prop[1]*cimag(gt[1])+prop[1]*prop[1]*cimag(gt[3]) + +prop[2]*prop[2]*cimag(gt[5]); +// v1[0]=prop[0]; v1[1]=prop[1]; v1[2]=prop[2]; +// cReflMatrVec(gt,v1,v2); +// tmp=cDotProd_Im(v2,v1); + C0dipole_refl=FOUR_PI*WaveNum*tmp; + printf("C0refl=%g\n",C0dipole_refl); + } return; case B_LMINUS: case B_DAVIS3: diff --git a/src/crosssec.c b/src/crosssec.c index 5b8d2eec..0807f1fe 100644 --- a/src/crosssec.c +++ b/src/crosssec.c @@ -948,6 +948,29 @@ double AbsCross(void) //====================================================================================================================== +double DecayCross(void) +// computes total cross section for the dipole incident field; similar to Cext +// 4pi*k*Im[p0(*).Escat(r0)] +{ + double sum; + size_t i; + + /* This is a correct expression only _if_ exciting p0 is real, then + * (using G(r1,r2) = G(r2,r1)^T, valid also for surface) + * p0(*).Escat_i(r0) = p0(*).G_0i.p_i = p_i.G_i0.p0(*) = p_i.G_i0.p0 = p_i.Einc_i + * => Im(p0(*).Escat(r0)) = sum{Im(P.E_inc)} + * + * For complex p0 an efficient calculation strategy (not to waste evaluations of interaction) is to compute an array + * of G_i0.p0(*) together with Einc and use it here afterwards. + */ + sum=0; + for (i=0;i ","Sets propagation direction of incident radiation, float. Normalization (to the unity " - "vector) is performed automatically.\n" + "vector) is performed automatically. For point-dipole incident beam this determines its direction.\n" "Default: 0 0 1",3,NULL}, {PAR(recalc_resid),"","Recalculate residual at the end of iterative solver.",0,NULL}, #ifndef SPARSE @@ -2198,9 +2198,12 @@ void PrintInfo(void) // log incident beam and polarization fprintf(logfile,"\n---In laboratory reference frame:---\nIncident beam: %s\n",beam_descr); fprintf(logfile,"Incident propagation vector: "GFORMDEF3V"\n",COMP3V(prop_0)); - fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY_0)); - fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n",COMP3V(incPolX_0)); - if (surface) { // include surface-specific vectors + if (beamtype==B_DIPOLE) fprintf(logfile,"(dipole orientation)\n"); + else { // polarizations are not shown for dipole incident field + fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY_0)); + fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n",COMP3V(incPolX_0)); + } + if (surface && beamtype==B_DIPOLE) { // include surface-specific vectors fprintf(logfile,"Reflected propagation vector: "GFORMDEF3V"\n",COMP3V(prIncRefl)); fprintf(logfile,"Transmitted propagation vector: "); if (msubInf) fprintf (logfile,"none\n"); @@ -2221,10 +2224,14 @@ void PrintInfo(void) if (beam_asym) fprintf(logfile,"Incident Beam center position: "GFORMDEF3V"\n", beam_center[0],beam_center[1],beam_center[2]); fprintf(logfile,"Incident propagation vector: "GFORMDEF3V"\n",COMP3V(prop)); - fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY)); - fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n\n",COMP3V(incPolX)); + if (beamtype==B_DIPOLE) fprintf(logfile,"(dipole orientation)\n"); + else { // polarizations are not shown for dipole incident field + fprintf(logfile,"Incident polarization Y(par): "GFORMDEF3V"\n",COMP3V(incPolY)); + fprintf(logfile,"Incident polarization X(per): "GFORMDEF3V"\n",COMP3V(incPolX)); + } } - else fprintf(logfile,"Particle orientation: default\n\n"); + else fprintf(logfile,"Particle orientation: default\n"); + fprintf(logfile,"\n"); } // log type of scattering matrices if (store_mueller) {