I am a bit confused on how the calculation of the expansion coefficients work for the spatial legendre filter. I ran the slab geometry example from the functional expansion jupyter notebook with an additional flux tally on a 1x1x40 mesh for comparison. The flux calculated from using the legendre polynomial expansion is the same shape as the flux calculated on the mesh but is about twice as large.

I have repeated this on a couple more geometries. The shape always matches but magnitude is always off by some (problem dependent) factor. Is there something that I am missing with regard to calculating the coefficients?

Here is the python code I ran to generate the plot from above:

```
import openmc
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 500
plt.rcParams['savefig.dpi'] = 500
fuel = openmc.Material()
fuel.add_element('U', 1.0, enrichment=4.5)
fuel.add_nuclide('O16', 2.0)
fuel.set_density('g/cm3', 10.0)
b4c = openmc.Material()
b4c.add_element('B', 4.0)
b4c.add_element('C', 1.0)
b4c.set_density('g/cm3', 2.5)
zmin, zmax = -10., 10.
box = openmc.model.RectangularPrism(10., 10., boundary_type='reflective')
bottom = openmc.ZPlane(z0=zmin, boundary_type='vacuum')
boron_lower = openmc.ZPlane(z0=-0.5)
boron_upper = openmc.ZPlane(z0=0.5)
top = openmc.ZPlane(z0=zmax, boundary_type='vacuum')
# Create three cells and add them to geometry
fuel1 = openmc.Cell(fill=fuel, region=-box & +bottom & -boron_lower)
absorber = openmc.Cell(fill=b4c, region=-box & +boron_lower & -boron_upper)
fuel2 = openmc.Cell(fill=fuel, region=-box & +boron_upper & -top)
geom = openmc.Geometry([fuel1, absorber, fuel2])
settings = openmc.Settings()
spatial_dist = openmc.stats.Box(*geom.bounding_box)
settings.source = openmc.IndependentSource(space=spatial_dist)
settings.batches = 210
settings.inactive = 10
settings.particles = 1000
flux_tally = openmc.Tally()
flux_tally.scores = ['flux']
# Create a Legendre polynomial expansion filter and add to tally
order = 16
expand_filter = openmc.SpatialLegendreFilter(order, 'z', zmin, zmax)
flux_tally.filters.append(expand_filter)
flux_mesh_tally = openmc.Tally()
flux_mesh_tally.scores = ['flux']
mesh = openmc.RegularMesh()
mesh.dimension = (1,1,40)
mesh.lower_left = geom.bounding_box[0]
mesh.upper_right = geom.bounding_box[1]
meshfilter = openmc.MeshFilter(mesh)
flux_mesh_tally.filters.append(meshfilter)
tallies = openmc.Tallies([flux_tally,flux_mesh_tally])
model = openmc.model.Model(geometry=geom, settings=settings, tallies=tallies)
sp_file = model.run(output=False)
#%%
with openmc.StatePoint(sp_file) as sp:
df1 = sp.tallies[flux_tally.id].get_pandas_dataframe()
df2 = sp.tallies[flux_mesh_tally.id].get_pandas_dataframe()
n = np.arange(order + 1)
a_n = (2*n + 1)/2 * df1['mean']
# divided by 10 for change of variables
phi10 = np.polynomial.Legendre(a_n/10, domain=(zmin, zmax))
# divided by 20 to match mesh tally
phi20 = np.polynomial.Legendre(a_n/20, domain=(zmin, zmax))
z = np.linspace(zmin, zmax, 1000)
fig,ax = plt.subplots()
ax.plot(z,phi10(z),label='expansion tally',alpha=0.75)
ax.plot(z,phi20(z),label='expansion tally halved',alpha=0.75)
z2 = np.linspace(zmin, zmax,41)
z2 = z2[:-1] + np.diff(z2)/2
ax.plot(z2,df2['mean'],label='mesh tally',alpha=0.75)
ax.set_xlabel('z(cm)')
ax.set_ylabel('flux')
plt.legend()
```