A non‐singular analytical theory for the motion of near‐Earth satellite orbits with air drag is developed in terms of KS elements, utilizing a model of the atmosphere that allows for oblateness and in which the density behaviour approximates to the observed diurnal variation. The series expansions include up to cubic terms in e (eccentricity) and c (a small parameter dependent on the flattening of the atmosphere). Only two of the nine equations are solved analytically to compute the state vector at the end of each revolution due to symmetry in the KS element equations. Numerical comparisons between the analytical and numerically integrated values of the drag perturbed orbital parameters (semi‐major axis, eccentricity and argument of perigee) are made up to 100 revolutions for the test cases with different inclinations, over a range of perigee and apogee heights, which experience significant diurnal bulge effects. The comparisons are found to be good. The accuracies of the numerical computations are checked with the bilinear relation in KS elements, and are found to be very satisfactory.