нарисовать спутниковую зону на равносторонней проекции

Мне нужно нарисовать границы зоны наблюдения спутника на равносторонней проекции. Я нашел эту формулу (1) и рисунок:

sin(fi) = cos(alpha) * sin(fiSat) – sin(alpha) * sin (Beta) * cos (fiSat);

sin(lambda) = (cos(alpha) * cos(fiSat) * sin(lambdaSat)) / cos(asin(sin(fi))) +
(sin(alpha) * sin(Beta) * sin(fiSat) * sin(lambdaSat)) / cos(asin(sin(fi))) -
(sin(alpha) * cos(Beta) * cos(lambdaSat))/cos(asin(sin(fi)));

cos(lambda) = (cos(alpha) * cos(fiSat) * cos(lambdaSat)) / cos(asin(sin(fi))) +
(sin(alpha) * sin(Beta) * sin(fiSat) * cos(lambdaSat)) / cos(asin(sin(fi))) -
(sin(alpha) * cos(Beta) * sin(lambdaSat)) / cos(asin(sin(fi)));

Сечения Земли в различных плоскостях:

И система уравнений (2) с рисунком:

if sin(lambda) > 0, cos(lambda) > 0 then lambda = asin(sin(lambda));
if sin(lambda) > 0, cos(lambda) < 0 then lambda = 180 - asin(sin(lambda));
if sin(lambda) < 0, cos(lambda) < 0 then lambda = 180 - asin(sin(lambda));
if sin(lambda) < 0, cos(lambda) > 0 then lambda = asin(sin(lambda));

Схема опорных углов для долготы Земли:

Where:  alpha – polar angle;
fiSat, lambdaSat – latitude, longitude of satellite;
Beta – angle which change from 0 to 2*Pi and help to draw the observation zone;
fi, lambda – latitude, longitude of point B on the border of observation zone;

Я повторяю обе формулы (1) и (2) в цикле от 0 до 2 * Pi, чтобы нарисовать границу зоны наблюдения. Но я не совсем уверен в (2) системе уравнений.

Внутри интервалов [-180; -90], [-90; 90], [90; 180] зона рисуется правильно.

Центр в -35; 45:

Центр в 120; 60:

Центр в -120; -25

Но на границе -90 и 90 градусов это становится грязно:

Центр в -95; -50

Центр в 95; 30

Можете ли вы помочь мне с формулами (1) и (2) или написать другие?

double deltaB = 1.0*M_PI/180;
observerZone.clear();
for (double Beta = 0.0; Beta <= (M_PI * 2) ; Beta += deltaB){
double sinFi = cos(alpha) * sin(fiSat) - sin(alpha) * sin(Beta) * cos(fiSat);
double sinLambda = (cos(alpha) * cos(fiSat) * sin(lambdaSat))/cos(asin(sinFi)) +
(sin(alpha) * sin(Beta) * sin(fiSat) * sin(lambdaSat))/cos(asin(sinFi)) -
(sin(alpha) * cos(Beta) * cos(lambdaSat))/cos(asin(sinFi));
double cosLambda = (cos(alpha) * cos(fiSat) * cos(lambdaSat))/cos(asin(sinFi)) +
(sin(alpha) * sin(Beta) * sin(fiSat) * cos(lambdaSat))/cos(asin(sinFi)) -
(sin(alpha) * cos(Beta) * sin(lambdaSat))/cos(asin(sinFi));
if (sinLambda > 0) {
if (cosLambda > 0 ){
sinLambda = asin(sinLambda);
sinFi = asin(sinFi);
}
else {
sinLambda = M_PI - asin(sinLambda);
sinFi = asin(sinFi);
}
}
else if (cosLambda > 0) {
sinLambda = asin(sinLambda);
sinFi = asin(sinFi);
}
else {
sinLambda = -M_PI - asin(sinLambda);
sinFi = asin(sinFi);
}
Point point;
point.latitude = qRadiansToDegrees(sinFi);
point.longitude = qRadiansToDegrees(sinLambda);
observerZone.push_back(point);
}