close all
clc

Angles   = 90.00001:0.00001:269.99999;
SinAngle = sin(Angles * unitsratio('rad', 'deg'));
CosAngle = cos(Angles * unitsratio('rad', 'deg'));

% AngularRadius = [87.25, 95.2, 139.0]';
AngularRadius = [45.0, 85.25, 89.25, 93.2, 97.2, 137.0, 141.0]';

SphereRadius = 1;

CircleX0 = SphereRadius * cos(AngularRadius * unitsratio('rad', 'deg'));
CircleR  = SphereRadius * sin(AngularRadius * unitsratio('rad', 'deg'));

for i = 1:length(AngularRadius)

    CircleX(i, :) = ones(1, length(CosAngle)) * CircleX0(i);
    CircleY(i, :) = CircleR(i) * CosAngle;
    CircleZ(i, :) = CircleR(i) * SinAngle;

end

[SphereX, SphereY, SphereZ] = sphere(90);

figure
hold on
grid off
axis equal

mesh(SphereX, SphereY, SphereZ, 'edgecolor', 'b');

for i = 1:length(AngularRadius)

    plot3(CircleX(i, :), CircleY(i, :), CircleZ(i, :), 'r', 'linewidth', 2)

end

view(3)

[CircleTheta, CirclePhi, CircleR] = cart2sph(CircleX, CircleY, CircleZ);

FigureSize = 1000;

CircleTheta = CircleTheta * unitsratio('deg', 'rad') * (FigureSize / 180);
CirclePhi   = CirclePhi   * unitsratio('deg', 'rad') * (FigureSize / 180);

figure
hold on
grid off
axis off
axis tight
axis fill
axis image

for i = 1:length(AngularRadius)

    plot( CircleTheta(i, :), CirclePhi(i, :), 'k')
    plot(-CircleTheta(i, :), CirclePhi(i, :), 'k')

end

xlim([-FigureSize FigureSize]);
ylim([-FigureSize / 2 FigureSize / 2]);

set(gcf, 'color', 'white')

set(gca, 'position', [0 0 1 1])

print('-dpng', '-r150', 'light_zone_borders.png')
