refresh_periodic_zero_mode_state Subroutine

  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
    !> Arrays already sized for this plan are reused; every component is overwritten.
    type(periodic_zero_mode_state_type), intent(inout) :: state
    real(dp) :: z0, z1, z2, denominator, offset, shift
    integer(i32) :: elem, interval, ninterval, i0, i1, i2, gauge_interval, gauge_breakpoint

    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
    if (allocated(state%cumulative_charge_coeff)) then
      if (size(state%cumulative_charge_coeff, 2) /= ninterval) deallocate (state%cumulative_charge_coeff)
    end if
    if (.not. allocated(state%cumulative_charge_coeff)) allocate (state%cumulative_charge_coeff(3, ninterval))
    call size_break_array(state%primitive_at_break, plan%nbreak)
    call size_break_array(state%sheet_charge, plan%nbreak)
    ! The per-interval differences accumulate in cumulative_charge_coeff and become prefix sums in place.
    state%cumulative_charge_coeff = 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)

    associate (difference => state%cumulative_charge_coeff)
      ! 傾斜 panel の寄与は、加える区間と外す区間それぞれの局所原点で表す。係数は panel 自身の高さ幅だけで
      ! 決まり、絶対 z に依らない（平行移動で値が変わらない）。
      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
            ! [z0, z1]: q*(z-z0)^2/denominator。原点 a では z-z0 = t+offset, offset=a-z0。
            denominator = (z1 - z0)*(z2 - z0)
            offset = plan%break_z(i0) - z0
            difference(:, i0 + 1_i32) = difference(:, i0 + 1_i32) + &
                                        charge(elem)*[offset*offset, 2.0_dp*offset, 1.0_dp]/denominator
            offset = plan%break_z(i1) - z0
            difference(:, i1 + 1_i32) = difference(:, i1 + 1_i32) - &
                                        charge(elem)*[offset*offset, 2.0_dp*offset, 1.0_dp]/denominator
          end if
          if (i2 > i1) then
            ! [z1, z2]: q*(1-(z2-z)^2/denominator)。原点 a では z2-z = offset-t, offset=z2-a。
            denominator = (z2 - z0)*(z2 - z1)
            offset = z2 - plan%break_z(i1)
            difference(:, i1 + 1_i32) = difference(:, i1 + 1_i32) + charge(elem)*[ &
                                        1.0_dp - offset*offset/denominator, 2.0_dp*offset/denominator, &
                                        -1.0_dp/denominator]
            offset = z2 - plan%break_z(i2)
            difference(:, i2 + 1_i32) = difference(:, i2 + 1_i32) - charge(elem)*[ &
                                        1.0_dp - offset*offset/denominator, 2.0_dp*offset/denominator, &
                                        -1.0_dp/denominator]
          end if
        end if
        difference(1, i2 + 1_i32) = difference(1, i2 + 1_i32) + charge(elem)
      end do

      ! 前の区間の多項式を次の区間の原点へ移してから差分を足す。移す量は隣接 break の間隔だけ。
      do interval = 2, ninterval
        shift = zero_mode_interval_origin(plan%break_z, interval) - &
                zero_mode_interval_origin(plan%break_z, interval - 1_i32)
        difference(:, interval) = shift_polynomial(difference(:, interval - 1), shift) + difference(:, interval)
      end do
    end associate
    ! Outside the source slab the integrated charge is known exactly. Do not
    ! extrapolate cancellation residue from the inclined-panel polynomials.
    state%cumulative_charge_coeff(:, 1) = 0.0_dp
    state%cumulative_charge_coeff(:, ninterval) = [state%total_charge, 0.0_dp, 0.0_dp]
    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_local_polynomial( &
                                               state%cumulative_charge_coeff(:, interval + 1), &
                                               plan%break_z(interval + 1) - plan%break_z(interval) &
                                               )
    end do
    call locate_zero_mode_interval(plan%break_z, plan%height_tolerance, z_gauge, gauge_interval, gauge_breakpoint)
    state%gauge_primitive = zero_mode_charge_primitive(plan, state, z_gauge, gauge_interval)
  end subroutine refresh_periodic_zero_mode_state