Energies don't add up with UsePolarizationDecomposition

[sbyrnes321]

With UsePolarizationDecomposition(), I often get the situation where reflected plus transmitted power is significantly different (higher or lower) than incident power, even though the structure is purely dielectric so there shouldn't be any loss or gain. Here is an example:

S = S4.NewSimulation()

S:SetLattice({0.7,0}, {0, 0.3})
S:SetNumG(50)
S:AddMaterial("Glass", {1.5^2,0})
S:AddMaterial("Air", {1,0})

S:AddLayer('Substrate', 0, 'Glass')
S:AddLayer('PatternLayer', 0.3, 'Air')
S:AddLayer('Air', 0, 'Air')

--S:SetLayerPatternEllipse(layer, material, center, angle, halfwidths)
S:SetLayerPatternEllipse('PatternLayer', 'Glass', {.45,0}, 0, {0.2,0.2})

--S:SetLayerPatternPolygon(layer, material, center, angle, vertices)
S:SetLayerPatternPolygon('PatternLayer', 'Glass', {0,0}, 0, {0,0,  0.2,0,  0.2,0.2,  0,0.2})

S:SetFrequency(2)

S:SetExcitationPlanewave(
    {0,0},  -- incidence angles (spherical coordinates: phi in [0,180], theta in [0,360])
    {1,0},  -- s-polarization amplitude and phase (in degrees)
    {1,0}) -- p-polarization amplitude and phase

S:UsePolarizationDecomposition()

incident_power, total_reflected_power = S:GetPowerFlux('Substrate', 0)
total_transmitted_power, this_should_be_zero = S:GetPowerFlux('Air', 0)
assert(this_should_be_zero == 0) -- this part works fine

print('total reflected power:', total_reflected_power)
print('total transmitted power', total_transmitted_power)
print('total reflected plus transmitted power:', math.abs(total_reflected_power) + math.abs(total_transmitted_power))
print('...whereas the incident power is:', incident_power)

Then the output (using the Windows .exe) is:

total reflected power:  -0.11690432308553
total transmitted power 1.4025400903532
total reflected plus transmitted power: 1.5194444134387
...whereas the incident power is:       1.3333333333333

If I understand correctly, the last two lines should be equal. Indeed, if I remove the UsePolarizationDecomposition line, then they are truly equal down to the last decimal point:

.......UsePolarizationDecomposition turned off.......
total reflected power:  -0.059337095586464
total transmitted power 1.2739962377469
total reflected plus transmitted power: 1.3333333333333
...whereas the incident power is:       1.3333333333333

Do you agree? Thanks so much for your hard and generous work, and thanks in advance for your help here. :-D

[sbyrnes321]

oops, that was a bad example: The ellipse intersects its periodic replicas, and the documentation states that this is not allowed. And if I remove the intersection then the energies almost add up (1e-3 accuracy or better). I think I have seen bigger discrepancies without self-intersections, but I can't find a simple example right now. I'll close this issue until I do. Sorry about that!