subroutine refresh_periodic_zero_mode_state(plan, charge, e_bottom, z_gauge, phi_gauge, state)
type(periodic_zero_mode_plan_type), intent(in) :: plan
real(dp), intent(in) :: charge(:)
real(dp), intent(in) :: e_bottom, z_gauge, phi_gauge
type(periodic_zero_mode_state_type), intent(out) :: state
real(dp), allocatable :: difference(:, :)
real(dp) :: z0, z1, z2, denominator, coefficient(3)
integer(i32) :: elem, interval, ninterval, i0, i1, i2
if (size(charge) /= plan%nelem) error stop 'periodic zero-mode charge size mismatch.'
if (any(.not. ieee_is_finite(charge))) error stop 'periodic zero-mode charges must be finite.'
ninterval = plan%nbreak + 1_i32
allocate (difference(3, ninterval + 1_i32))
allocate (state%cumulative_charge_coeff(3, ninterval))
allocate (state%primitive_at_break(plan%nbreak), state%sheet_charge(plan%nbreak))
difference = 0.0_dp
state%sheet_charge = 0.0_dp
state%e_bottom = e_bottom
state%z_gauge = z_gauge
state%phi_gauge = phi_gauge
state%total_charge = sum(charge)
do elem = 1, plan%nelem
z0 = plan%panel(elem)%z(1)
z1 = plan%panel(elem)%z(2)
z2 = plan%panel(elem)%z(3)
i0 = plan%panel(elem)%break_index(1)
i1 = plan%panel(elem)%break_index(2)
i2 = plan%panel(elem)%break_index(3)
if (plan%panel(elem)%horizontal) then
state%sheet_charge(i0) = state%sheet_charge(i0) + charge(elem)
else
if (i1 > i0) then
denominator = (z1 - z0)*(z2 - z0)
coefficient = charge(elem)*[z0*z0/denominator, -2.0_dp*z0/denominator, 1.0_dp/denominator]
call add_interval_range(difference, i0 + 1_i32, i1, coefficient)
end if
if (i2 > i1) then
denominator = (z2 - z0)*(z2 - z1)
coefficient = charge(elem)*[ &
1.0_dp - z2*z2/denominator, 2.0_dp*z2/denominator, -1.0_dp/denominator &
]
call add_interval_range(difference, i1 + 1_i32, i2, coefficient)
end if
end if
call add_interval_range(difference, i2 + 1_i32, ninterval, [charge(elem), 0.0_dp, 0.0_dp])
end do
state%cumulative_charge_coeff(:, 1) = difference(:, 1)
do interval = 2, ninterval
state%cumulative_charge_coeff(:, interval) = &
state%cumulative_charge_coeff(:, interval - 1) + difference(:, interval)
end do
state%primitive_at_break(1) = 0.0_dp
do interval = 1, plan%nbreak - 1
state%primitive_at_break(interval + 1) = state%primitive_at_break(interval) + &
integrate_polynomial( &
state%cumulative_charge_coeff(:, interval + 1), &
plan%break_z(interval), plan%break_z(interval + 1) &
)
end do
end subroutine refresh_periodic_zero_mode_state