k = 2*Pi
ph = ArcTan[x, y]
f = 1.2*10^-7
l = 5.32*10^-7
NA = 0.85
th1 = ArcSin[NA]
ybc = 8.854187*10^-12
c = 3*10^8
P = 50
R = 2.67*10^-1
E0 = P/(Pi*ybc*c*(R)^2)
arf = 0.5
r = Sqrt[x^2 + y^2]
ph0 = 0 Pi/180
ph1 = 90 Pi/180
ph2 = 180 Pi/180
Bm = BesselJ[n, k*r*Sin[th]]
EZ = (-2*f*E0)/l*
Sum[Sin[th]*Sin[th]*Sqrt[Cos[th]]*Cos[th]*
Exp[I*k*z*
Cos[th]]*((Cos[ph0]*
Sum[((arf - 1)*Cos[arf*Pi] + I*n*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^(n + 1)*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}]) - (Sin[ph0]*
Sum[(n*Cos[arf*Pi] + I*(arf - 1)*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^n*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}])), {th, 0, th1, 0.1}]
PSZ = Im[D[Conjugate[EZ]*EZ]]
POZ = Im[Conjugate[EZ]*D[EZ, ph]]
PZ = PSZ + POZ
EZ1 = (-2*f*E0)/l*
Sum[Sin[th]*Sin[th]*Sqrt[Cos[th]]*Cos[th]*
Exp[I*k*z*
Cos[th]]*((Cos[ph1]*
Sum[((arf - 1)*Cos[arf*Pi] + I*n*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^(n + 1)*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}]) - (Sin[ph1]*
Sum[(n*Cos[arf*Pi] + I*(arf - 1)*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^n*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}])), {th, 0, th1, 0.1}]
PSZ1 = Im[D[Conjugate[EZ1]*EZ1]]
POZ1 = Im[Conjugate[EZ1]*D[EZ1, ph]]
PZ1 = PSZ1 + POZ1
EZ2 = (-2*f*E0)/l*
Sum[Sin[th]*Sin[th]*Sqrt[Cos[th]]*Cos[th]*
Exp[I*k*z*
Cos[th]]*((Cos[ph2]*
Sum[((arf - 1)*Cos[arf*Pi] + I*n*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^(n + 1)*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}]) - (Sin[ph2]*
Sum[(n*Cos[arf*Pi] + I*(arf - 1)*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^n*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}])), {th, 0, th1, 0.1}]
PSZ2 = Im[D[Conjugate[EZ2]*EZ2]]
POZ2 = Im[Conjugate[EZ2]*D[EZ2, ph]]
PZ2 = PSZ2 + POZ2
z = 0
Plot[{PZ, PZ1, PZ2}, {y, -2, 2}, PlotRange -> All,
PlotLegends -> {"\!\(\*SubscriptBox[\(\[CurlyPhi]\), \(0\)]\)=0",
"\!\(\*SubscriptBox[\(\[CurlyPhi]\), \(0\)]\)=\!\(\*FractionBox[\(\
\[Pi]\), \(2\)]\)",
"\!\(\*SubscriptBox[\(\[CurlyPhi]\), \(0\)]\)=\[Pi]"},
Frame -> True]
ph = ArcTan[x, y]
f = 1.2*10^-7
l = 5.32*10^-7
NA = 0.85
th1 = ArcSin[NA]
ybc = 8.854187*10^-12
c = 3*10^8
P = 50
R = 2.67*10^-1
E0 = P/(Pi*ybc*c*(R)^2)
arf = 0.5
r = Sqrt[x^2 + y^2]
ph0 = 0 Pi/180
ph1 = 90 Pi/180
ph2 = 180 Pi/180
Bm = BesselJ[n, k*r*Sin[th]]
EZ = (-2*f*E0)/l*
Sum[Sin[th]*Sin[th]*Sqrt[Cos[th]]*Cos[th]*
Exp[I*k*z*
Cos[th]]*((Cos[ph0]*
Sum[((arf - 1)*Cos[arf*Pi] + I*n*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^(n + 1)*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}]) - (Sin[ph0]*
Sum[(n*Cos[arf*Pi] + I*(arf - 1)*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^n*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}])), {th, 0, th1, 0.1}]
PSZ = Im[D[Conjugate[EZ]*EZ]]
POZ = Im[Conjugate[EZ]*D[EZ, ph]]
PZ = PSZ + POZ
EZ1 = (-2*f*E0)/l*
Sum[Sin[th]*Sin[th]*Sqrt[Cos[th]]*Cos[th]*
Exp[I*k*z*
Cos[th]]*((Cos[ph1]*
Sum[((arf - 1)*Cos[arf*Pi] + I*n*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^(n + 1)*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}]) - (Sin[ph1]*
Sum[(n*Cos[arf*Pi] + I*(arf - 1)*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^n*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}])), {th, 0, th1, 0.1}]
PSZ1 = Im[D[Conjugate[EZ1]*EZ1]]
POZ1 = Im[Conjugate[EZ1]*D[EZ1, ph]]
PZ1 = PSZ1 + POZ1
EZ2 = (-2*f*E0)/l*
Sum[Sin[th]*Sin[th]*Sqrt[Cos[th]]*Cos[th]*
Exp[I*k*z*
Cos[th]]*((Cos[ph2]*
Sum[((arf - 1)*Cos[arf*Pi] + I*n*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^(n + 1)*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}]) - (Sin[ph2]*
Sum[(n*Cos[arf*Pi] + I*(arf - 1)*Sin[arf*Pi])/((arf - 1)^2 -
n^2)*I^n*Exp[I*n*ph]*Sin[arf*Pi]*Bm, {n, -10, 10,
1}])), {th, 0, th1, 0.1}]
PSZ2 = Im[D[Conjugate[EZ2]*EZ2]]
POZ2 = Im[Conjugate[EZ2]*D[EZ2, ph]]
PZ2 = PSZ2 + POZ2
z = 0
Plot[{PZ, PZ1, PZ2}, {y, -2, 2}, PlotRange -> All,
PlotLegends -> {"\!\(\*SubscriptBox[\(\[CurlyPhi]\), \(0\)]\)=0",
"\!\(\*SubscriptBox[\(\[CurlyPhi]\), \(0\)]\)=\!\(\*FractionBox[\(\
\[Pi]\), \(2\)]\)",
"\!\(\*SubscriptBox[\(\[CurlyPhi]\), \(0\)]\)=\[Pi]"},
Frame -> True]
