A method for generating long periodic variations in satellite orbital elements when perturbed by discrete gravity anomalies is presented. The method consists of developing a disturbing potential as a function of orbital and gravity anomaly parameters, and generating partial derivatives of the potential with respect to the orbital elements. The partials are averaged over the period of the satellite to eliminate the short periodic variations. The averaged partials are substituted into the variation of parameter equations to give the mean orbital rates. Classically orbital elements are used in generating gravity field and thus the method is dynamic in nature. The problem is extremely cumbersome and complex when multi-state parameters have to be estimated from a considerably large data set. However, when mean orbital rates are used, the problem reduces to a simple linear static case, where only the gravity parameters have to be estimated, and it is a simple matrix inversion problem. Thus the method developed here was utilized in reducing Appolo 15 and 16 subsatellite radio tracking data to produce a lunar gravity field represented by point masses.