We obtain a theoretical description of the height (z) distribution of flare hard X-rays in the collisional thick-target model as a function of photon energy epsilon. This depends on the target atmosphere density structure n(z) and on the beam spectral index delta. We show that by representing the data in terms of the 1-D function z(max) (epsilon) defining where the emission peaks as a function of epsilon it is possible to derive n(z) from data on zmax(E). This is done first on the basis of a simple stopping depth argument then refined to allow for the dependence on spectral index delta. The latter is worked out in detail for the case of a parameterization n(z) = n(o) (z/z(o))(-b) which yields numerical results for zmax(epsilon) well fit by z(max)(epsilon) similar to epsilon(-alpha), with alpha dependent on delta, which is also found to fit well to actual observations. This enables derivation of flare loop n(z) in terms of n(o), b from RHESSI data in an entirely novel way, independent of other density diagnostic methods, and also of how n(z) varies with time in flares such as by evaporation, as detailed in companion Paper II.