(\sin \theta _1 \sin \phi _1 \cos \theta _2 - \sin \theta _2 \sin \phi _2 \cos \theta _1) \sin \theta \cos \phi +
   + (\cos \theta _1 \sin \theta _2 \cos \phi _2 - \cos \theta _2 \sin \theta _1 \cos \phi _1) \sin \theta \sin \phi +
   + (\sin \theta _1 \cos \phi _1 \sin \theta _2 \sin \phi _2 - \sin \theta _2 \cos \phi _2 \sin \theta _1 \sin \phi _1) \cos \theta =
   = 0
твоя формула даст тебе тот же результат, только тебе придётся проверить на сколько сдвинуть арксинус или арккосинус чтобы попасть в область значений теты. просто в данном случае, арктангенс как раз покрывает ровно то что нужно, а деление на ноль означает арктангенс бесконечности то есть \pm \frac {\pi}2, там уже по непрерывности, иначе никак (в твоём методе тоже)