  The standard variational derivation of stellar matter structure in the
Wigner-Seitz approximation is generalized to the finite temperature situation
where a wide distribution of different nuclear species can coexist in the same
density and proton fraction condition, possibly out of $\beta$-equilibrium. The
same theoretical formalism is shown to describe on one side the single-nucleus
approximation (SNA), currently used in most core collapse supernova
simulations, and on the other side the nuclear statistical equilibrium (NSE)
approach, routinely employed in r- and p-process explosive nucleosynthesis
problems. In particular we show that in-medium effects have to be accounted for
in NSE to have a theoretical consistency between the zero and finite
temperature modeling. The bulk part of these in-medium effects is analytically
calculated and shown to be different from a van der Waals excluded volume term.
This unified formalism allows controlling quantitatively the deviations from
the SNA in the different thermodynamic conditions, as well as having a NSE
model which is reliable at any arbitrarily low value of the temperature, with
potential applications for neutron star cooling and accretion problems. We
present different illustrative results with several mass models and effective
interactions, showing the importance of accounting for the nuclear species
distribution even at temperatures lower than 1 MeV.
