A fast iterative method for gravity field determination from low Earth satellite orbit coordinates has been developed and implemented successfully. The method is based on energy conservation and avoids problems related to orbit dynamics and initial state. In addition, the particular geometry of a repeat orbit is exploited by using a very efficient iterative estimation scheme, in which a set of normal equations is approximated by a sparse block-diagonal equivalent. Recovery experiments for spherical harmonic gravity field models up to degree and order 80 and 120 were conducted based on a 29-day simulated data set of orbit coordinates. The method was found to be very flexible and could be easily adapted to include observations of non-conservative accelerations, such as (to be) provided by satellites like CHAMP, GRACE, and GOCE. A serious drawback of the method is its large sensitivity to satellite velocity errors. Existing orbit determination strategies need to be altered or augmented to include algorithms that focus on optimizing the accuracy of estimated velocities.