Hello!
I am comparing the flux of the whole simulation cell using meep.get_fluxes with
the far-field flux out from a sphere with a radius much larger than the length
of the sim. cell. However, I get a complete mismatch of the flux values where I
believe they should be close to same and I can't figure out the error.
Here is a scaled down version of the code:
sx=sy=sz=4
dpml=1.6
symmetry=[mp.Mirror(mp.X),mp.Mirror(mp.Y,phase=-1)]
pml_layer=[mp.PML(dpml)]
padding=0.1
fcen=2
df=1.2
nfreq=4
source2=[mp.Source(mp.ContinuousSource(frequency=fcen,start_time=0,end_time=1),
component=mp.Ey,
center=mp.Vector3(0,0,0))]
sim=mp.Simulation(cell_size=
mp.Vector3(sx+2*dpml,sy+2*dpml,sz+2*dpml),
symmetries=symmetry,
sources=source2,
dimensions=3,
boundary_layers=pml_layer,
resolution=10)
fluxregion1=mp.FluxRegion( #region -y to calculate flux from
center=mp.Vector3(0,-sy/2+padding,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y,
weight=-1)
fluxregion2=mp.FluxRegion( #region y to calculate flux from
center=mp.Vector3(0,sy/2-padding,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y)
fluxregion3=mp.FluxRegion( # region -x to calculate flux from
center=mp.Vector3(-sx/2+padding,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X,
weight=-1)
fluxregion4=mp.FluxRegion( #region x to calculate flux from
center=mp.Vector3(sx/2-padding,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X)
fluxregion5=mp.FluxRegion( #z-bottom region to calculate
flux from
center=mp.Vector3(0,0,sz/2-padding),
size=mp.Vector3(sx-padding*2,sz-padding*2,0),
direction=mp.Z)
fluxregion6=mp.FluxRegion( #z-top region to calculate flux
from
center=mp.Vector3(0,0,-sz/2+padding),
size=mp.Vector3(sx-padding*2,sy-padding*2,0),
direction=mp.Z,
weight=-1)
flux_total=sim.add_flux(fcen,df,nfreq,
fluxregion1,fluxregion2,fluxregion3,fluxregion4,fluxregion5,fluxregion6)
--------------
nearfieldregion1=mp.Near2FarRegion( #nearfield -z.
above pyramid.
center=mp.Vector3(0,-sz/2+padding,0),
size=mp.Vector3(sx-padding*2,sy-padding*2,0),
direction=mp.Z,
weight=-1)
nearfieldregion6=mp.Near2FarRegion( #nearfield z.
below pyramid.
center=mp.Vector3(0,sz/2+padding,0),
size=mp.Vector3(sx-padding*2,sy-padding*2,0),
direction=mp.Z)
nearfieldregion2=mp.Near2FarRegion(
center=mp.Vector3(sx/2,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X)
nearfieldregion3=mp.Near2FarRegion(
center=mp.Vector3(-sx/2,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X,
weight=-1)
nearfieldregion4=mp.Near2FarRegion(
center=mp.Vector3(0,sy/2,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y)
nearfieldregion5=mp.Near2FarRegion(
center=mp.Vector3(0,-sy/2,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y,
weight=-1)
nearfield=sim.add_near2far(fcen,df,nfreq,nearfieldregion1,nearfieldregion2,nearfieldregion3,
nearfieldregion4,nearfieldregion5,nearfieldregion6)
-------------
sim.run(until=10)
P_tot_ff = np.zeros(nfreq)
r=100
npts=10
Px=0
Py=0
Pz=0
for n in range(npts):
angleN=(n/npts)*math.pi*2
for m in range(npts):
angleM=(m/npts)*math.pi
ff=sim.get_farfield(nearfield,
mp.Vector3(r*math.sin(angleM)*math.cos(angleN),r*math.sin(angleM)*math.sin(angleN),-r*math.cos(angleM)))
i=0
for k in range(nfreq):
Px=(ff[i+1]*np.conjugate(ff[i+5])-ff[i+2]*np.conjugate(ff[i+4]))
Py=(ff[i+2]*np.conjugate(ff[i+3])-ff[i]*np.conjugate(ff[i+5]))
Pz=(ff[i]*np.conjugate(ff[i+4])-ff[i+1]*np.conjugate(ff[i+3]))
Pr=(Px)*math.cos(angleN)*math.sin(angleM)+(Py)*math.sin(angleN)*math.sin(angleM)-(Pz)*math.cos(angleM)
surface_Element=math.pow(r,2)*math.pow(angle,2)*math.pow((1/npts),2)
P_tot_ff[k] += surface_Element*(1/2)*np.real(Pr)
print(n,m,k,P_tot_ff[k])
i=i+6
flux_tot_out = mp.get_fluxes(flux_total)
P_tot_ff= [0.00021628 0.00140266 0.00369171 0.00201077]
flux_tot_out = [0.031127320105939483, 0.19604245927203143, 0.3616725444974181,
0.20233403694380303]
P_tot_ff should be the radial component of the flux summed over the whole
sphere. So in my mind, P_tot_ff and flux_tot_out should be very close to the
same, but they obviously are not?
I would be very grateful for some help.
Warm Greetings,
Karl-Johan
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
