Исходное уравнение геодезической, даже на поверхности вращения -- это уравнение второго порядка:

Здесь

-- криволинейные координаты на поверхности.
В случае поверхности вращения справедлива теорема Клеро: вдоль геодезической сохраняется произведение:
(расстояние от точки до оси вращения)

(косинус угла между геодезической и параллелью),
или, что то же
(расстояние от точки до оси вращения)

(синус угла между геодезической и меридианом)
Если ось вращения

, то расстояние

от точки до оси через декартовы координаты выражается так:

. Декартовы компоненты единичного вектора касательной к геодезической равны

. Единичный вектор, направленный по параллели (координатной линии

в цилиндрических координатах), равен

. Значит, упомянутый косинус угла равен скалярному произведению этих двух векторов, т.е.

. Теорема Клеро гласит, что если умножить его на

, то произведение будет константой на геодезической:

А это -- практически Ваше уравнение.
Левая часть этого уравнения,

, и есть первый интеграл уравнения геодезических (того, первоначального, второго порядка) для случая поверхности вращения.
Доказательство теоремы Клеро можно посмотреть здесь:
http://higeom.math.msu.su/people/chernavski/lectures.pdfСтраницы 75, 76.