% sphere.inc 

/Pi 3.14159265358979 def

% matrix A => axis of A
% assumes (Ae0 - e0)x(Ae1 - e1)

/axis {
4 dict begin
/A exch def

[
A 0 get 0 get 1 sub
A 1 get 0 get
A 2 get 0 get
]

[
A 0 get 1 get
A 1 get 1 get 1 sub
A 2 get 1 get
]

cross-product

end
} def

% [w0 w1] t
% [w0 w1] of same length, w0 perp to w1 - the arc from w0 through w1
% cos(t).w0 + sin(t).w1 
% velocity -sin(t).w0 + cos(t).w1

/great-circle {
8 dict begin
/t exch def
/pars exch def
/w0 pars 0 get def
/w1 pars 1 get def
/t t 180 mul Pi div def
[
[ 
t cos w0 0 get mul
t sin w1 0 get mul
add
t cos w0 1 get mul
t sin w1 1 get mul
add
t cos w0 2 get mul
t sin w1 2 get mul
add
]
[
t sin neg w0 0 get mul
t cos w1 0 get mul
add
t sin neg w0 1 get mul
t cos w1 1 get mul
add
t sin neg w0 2 get mul
t cos w1 2 get mul
add
]
]
end
} def

% axis
% assuming axis x [0 0 -1] not [0 0 0]
% where axis-circle goes from visible to invisible in projection

/horizon {
[0 0 -1]
cross-product
normalized
end
} def

% locmat w0 w1 N

/sphere-arc {
8 dict begin
/N exch def
/w1 exch def
/w0 exch def
/locmat exch def

/theta w0 w1 angle-between degrees-to-radians def
/alpha w0 w1 cross-product normalized def
/wperp alpha w0 cross-product def
locmat [w0 wperp] /great-circle 0 theta N mkpath3d 
end
} def

% [R phi theta]
% phi = latitude 

/spherical-to-rect {
4 dict begin
aload pop
/theta exch def
/phi exch def
/R exch def
[
R theta cos mul phi cos mul 
R theta sin mul phi cos mul 
R phi sin mul 
]
end
} def

% w0 w1 => h = axis x [0 0 1] = last visible point on path w0 -> w1
% if w0 is visible

/horizon {
cross-product normalized
[0 0 1] cross-product normalized
} def

% R w0 w1

/true-horizon {
4 dict begin
/w1 exch def
/w0 exch def
/R exch def
R w0 matrix-vector
R w1 matrix-vector
horizon
R transpose exch matrix-vector
end
} def

% - testing ---------------------------------------------------

false {

/RUN {} def
RUN (../ps/mkpath.inc) run
RUN (../ps/mkpath3d.inc) run
RUN (../ps/matrix.inc) run

72 72 scale
4.25 5.5 translate
0.012 setlinewidth

gsave
newpath
-2 0 moveto
2 0 lineto
0 -2 moveto
0 2 lineto
% stroke
grestore

newpath
0 0 2 0 360 arc
closepath
stroke

/R 
[1 0 1] 20 rotation-matrix
def

use-projection

/tau [0 0 5] def

/A [2 0 0] def
/B [0 2 0] def
/C [0 0 2] def

gsave
0.8 setgray
newpath

[R tau]
A B
16 sphere-arc 

[R tau]
B C
16 sphere-arc 

[R tau]
C A
16 sphere-arc 

closepath
fill

0 setgray
newpath

[R tau]
A B
16 sphere-arc 

[R tau]
B C
16 sphere-arc 

[R tau]
C A
16 sphere-arc 

stroke
grestore

} if