  We find universal structure and scaling of BEC statistics and thermodynamics
for mesoscopic canonical-ensemble ideal gas in a trap for any parameters,
including critical region. We identify universal constraint-cut-off mechanism
that makes BEC fluctuations non-Gaussian and is responsible for critical
phenomena. Main result is analytical solution to problem of critical phenomena.
It is derived by calculating universal distribution of noncondensate occupation
(Landau function) and then universal functions for physical quantities. We find
asymptotics of that solution and its approximations which describe universal
structure of critical region in terms of parabolic cylinder or confluent
hypergeometric functions. Results for order parameter, statistics, and
thermodynamics match known asymptotics outside critical region. We suggest
2-level and 3-level trap models and find their exact solutions in terms of
cut-off negative binomial distribution (that tends to cut-off gamma
distribution in continuous limit) and confluent hypergeometric distribution. We
introduce a regular refinement scheme for condensate statistics approximations
on the basis of infrared universality of higher-order cumulants and method of
superposition and show how to model BEC statistics in actual traps. We find
that 3-level trap model with matching the first 4 or 5 cumulants is enough to
yield remarkably accurate results in whole critical region. We derive exact
multinomial expansion for noncondensate occupation distribution and find its
high temperature asymptotics (Poisson distribution). We demonstrate that
critical exponents and a few known terms of Taylor expansion of universal
functions, calculated previously from fitting finite-size simulations within
renorm-group theory, can be obtained from presented solutions.
