Multipole decomposition and scattering cross section using COMSOL

Below are the files you should save as notepad texts. Four files in total, please name them "Parameters", "Exact formulae", "Integral", "Results" and save them on your desktop. You will call them in subsequent events as I will show in the video.

Please notice that you can use the last three files for whatever nanoparticle you may have. The only changes you may need to update are, the area of your particle, so instead of (pi*rdisk^2) in the "Results" file, put your area. Also, please make sure you always name the component coupling as "intop1" for the nanoparticle in case you may have more than one component coupling, or in case COMSOL default name is changed in other versions.

The link for the paper is:
https://doi.org/10.1016/j.optcom.2017...
----------------------------------------------------------------------------------------------------------
Multipole calculations files:

Parameters:

rdisk 0.225[um] ""
hdisk 0.4[um] ""
R 0.75[um] ""

Exact formulae:

k ewfd.k0 ""
r sqrt(x^2+y^2+z^2) ""
kr k*r ""
Px ewfd.Px ""
Py ewfd.Py ""
Pz ewfd.Pz ""
iw ewfd.iomega ""
rdP x*Px+y*Py+z*Pz ""
EDx Px*jn(0,kr)+k^2/2*(3*rdP*x-Px*r^2)*jn(2,kr)/kr^2 ""
EDy Py*jn(0,kr)+k^2/2*(3*rdP*y-Py*r^2)*jn(2,kr)/kr^2 ""
EDz Pz*jn(0,kr)+k^2/2*(3*rdP*z-Pz*r^2)*jn(2,kr)/kr^2 ""
rcp_x y*Pz-z*Py ""
rcp_y z*Px-x*Pz ""
rcp_z x*Py-y*Px ""
MDx -3/2*iw*rcp_x*jn(1,kr)/kr ""
MDy -3/2*iw*rcp_y*jn(1,kr)/kr ""
MDz -3/2*iw*rcp_z*jn(1,kr)/kr ""
EQxx "3*( (3*(x*Px+x*Px)-2*rdP)*jn(1,kr)/kr + 2*k^2*jn(3,kr)/kr^3*(5*x*x*rdP-r^2*(x*Px+x*Px)-r^2*rdP) )" ""
EQyy "3*( (3*(y*Py+y*Py)-2*rdP)*jn(1,kr)/kr + 2*k^2*jn(3,kr)/kr^3*(5*y*y*rdP-r^2*(y*Py+y*Py)-r^2*rdP) )" ""
EQzz "3*( (3*(z*Pz+z*Pz)-2*rdP)*jn(1,kr)/kr + 2*k^2*jn(3,kr)/kr^3*(5*z*z*rdP-r^2*(z*Pz+z*Pz)-r^2*rdP) )" ""
EQyz "3*( (3*(y*Pz+z*Py))*jn(1,kr)/kr + 2*k^2*jn(3,kr)/kr^3*(5*y*z*rdP-r^2*(y*Pz+z*Py)) )" ""
EQxz "3*( (3*(x*Pz+z*Px))*jn(1,kr)/kr + 2*k^2*jn(3,kr)/kr^3*(5*x*z*rdP-r^2*(x*Pz+z*Px)) )" ""
EQxy "3*( (3*(x*Py+y*Px))*jn(1,kr)/kr + 2*k^2*jn(3,kr)/kr^3*(5*x*y*rdP-r^2*(x*Py+y*Px)) )" ""
MQxx -iw*15*(2*x*(y*Pz-z*Py))*jn(2,kr)/kr^2 ""
MQyy -iw*15*(2*y*(z*Px-x*Pz))*jn(2,kr)/kr^2 ""
MQzz -iw*15*(2*z*(x*Py-y*Px))*jn(2,kr)/kr^2 ""
MQyz -iw*15*(y*(x*Py-y*Px)+z*(z*Px-x*Pz))*jn(2,kr)/kr^2 ""
MQxz -iw*15*(x*(x*Py-y*Px)+z*(y*Pz-z*Py))*jn(2,kr)/kr^2 ""
MQxy -iw*15*(x*(z*Px-x*Pz)+y*(y*Pz-z*Py))*jn(2,kr)/kr^2 "“


Integral:

p_sq abs(intop1(EDx))^2+abs(intop1(EDy))^2+abs(intop1(EDz))^2 ""
m_sq (abs(intop1(MDx))^2+abs(intop1(MDy))^2+abs(intop1(MDz))^2)/c_const^2 ""
eq_sq (abs(intop1(EQxx))^2+abs(intop1(EQyy))^2+abs(intop1(EQzz))^2+2*(abs(intop1(EQyz))^2+abs(intop1(EQxz))^2+abs(intop1(EQxy))^2))*k^2 ""
mq_sq (abs(intop1(MQxx))^2+abs(intop1(MQyy))^2+abs(intop1(MQzz))^2+2*(abs(intop1(MQyz))^2+abs(intop1(MQxz))^2+abs(intop1(MQxy))^2))*k^2/c_const^2 "“


Results:

ewfd.k0^4/6/pi/epsilon0_const^2*p_sq/(pi*rdisk^2)/(1[V/m])^2 rad ED
ewfd.k0^4/6/pi/epsilon0_const^2*m_sq/(pi*rdisk^2)/(1[V/m])^2 rad MD
ewfd.k0^4/6/pi/epsilon0_const^2*eq_sq/(pi*rdisk^2)/120/(1[V/m])^2 rad EQ
ewfd.k0^4/6/pi/epsilon0_const^2*mq_sq/(pi*rdisk^2)/120/(1[V/m])^2 rad MQ
ewfd.k0^4/6/pi/epsilon0_const^2*(p_sq+m_sq+eq_sq/120+mq_sq/120)/(pi*rdisk^2)/(1[V/m])^2 rad TOT