#10.r.maple
#
#   Draws the hyperbolae from the second half of the seventh ideal
#
interface(quiet=true):
with(plots):
plotsetup(x11,plotoptions=`width=3in,height=2.5in`):

############################################
#
Sphere := (a,b,c,r) ->
  linalg[matrix]([
   [a^2+b^2+c^2-r^2, -a, -b, -c],
   [      -a     , 1 , 0 , 0 ],
   [      -b     , 0 , 1 , 0 ],
   [      -c     , 0 , 0 , 1 ]]):

Wedge := proc(Sph)
local Pairs, i, j, wedge:
 wedge := linalg[matrix](6, 6, 0):
 Pairs := combinat[choose]([1,2,3,4],2):
 for i from 1 to 6 do
  for j from 1 to 6 do
   wedge[i,j] := linalg[det](linalg[submatrix](Sph,Pairs[i],Pairs[j])):
  od:
 od:
 evalm(wedge)
end:

Equation := proc(a,b,c,r)
local  plucker:
 plucker:=linalg[vector]([p01,p02,p03,p12,p13,p23]):
 linalg[innerprod](plucker,Wedge(Sphere(a,b,c,r)),plucker)
end:

Quadric := (a,b,c,d,e,f,g,h,x,y) ->
  linalg[matrix]([
   [ a , b , c , d ],
   [ b , e , f , g ],
   [ c , f , h , x ],
   [ d , g , x , y ]]):

##########################################
#  We assume that A = b = 1.
#  a is the function a(B) stored in aSOLS
#
TanPoint := proc(Q,B,a)
local Point, plugin, touch:
 Point:=linalg[vector]([1,B,t*a,t]):
 plugin:=linalg[innerprod](Point,Q,Point):
 touch:= expand(-coeff(plugin,t)/2/coeff(plugin,t^2)):
# print(plugin, touch):
 [B, simplify(touch*a), touch]
end:
# [1,B,1,b]
#  B is the function B(b) stored in aSOLS
TanPointb := proc(Q,B,b)
local Point, plugin, touch:
 Point:=linalg[vector]([1,B,t,t*b]):
 plugin:=linalg[innerprod](Point,Q,Point):
 touch:= expand(-coeff(plugin,t)/2/coeff(plugin,t^2)):
 [B, simplify(touch), simplify(touch*b)]
end:

Coords:=linalg[vector]([1,X,Y,Z]):
plucker:=linalg[vector]([p01,p02,p03,p12,p13,p23]):

F:=subs(p02=aA*Aa,p03=bB*Aa,p12=aA*Bb,p13=bB*Bb,p01=0,p23=0,
    linalg[innerprod](plucker,Wedge(Quadric(-1,0,-1,0,1,0,0,1,0,-1)),plucker)):

aSOLS:=[solve(subs(Bb=B,bB=1,Aa=1,F)=0,aA)]:
bSOLS:=[solve(subs(Bb=B,aA=1,Aa=1,F)=0,bB)]:
######################################################################
#
s:=-2^(1/2):   # The `sister quadric to this has s:=2^(1/2):
s2:=-s:
Q10:=Quadric(s*g^2, (1+s)*g^2/2, s*g, g, 
            g^2, g, g, 4-2*s, 1, 2*(s-1)):

##################################################################
#     [ r ,  0 ,  0 ,  0 ],
#     [ a ,  d ,  0 ,  0 ],     The list Eqs is for this 
#     [ b ,  e ,  x ,  0 ],     matrix (See parameter.maple
#     [ c ,  f ,  y ,  z ]])):
Eqs:=[r^2*g^2+8*s+12, 2*a-r*s-r, 2*b+r*s*g+r*g, 2*c+r*s*g+r*g, 
      2*z^2-s-1, 2*y^2+29*s+41, x+2*y*s-2*y, f^2-12*s-17, 
      e+f, 4*d-12*f*r^2*s*g+17*f*r^2*g]:
r:=simplify([solve(Eqs[1]=0,r)][1]):   #  The solutions have same square
a:=solve(Eqs[2]=0,a):
b:=solve(Eqs[3]=0,b):
c:=solve(Eqs[4]=0,c):
z:=[solve(Eqs[5]=0,z)][1]:    #  The solutions have same square
y:=[solve(Eqs[6]=0,y)][1]:    #  The solutions have same square
x:=solve(Eqs[7]=0,x):
f:=[solve(Eqs[8]=0,f)][1]:    #  The solutions have same square
e:=solve(Eqs[9]=0,e):
d:=solve(Eqs[10]=0,d):

A:=linalg[inverse](linalg[matrix]([
     [ r ,  0 ,  0 ,  0 ],
     [ a ,  d ,  0 ,  0 ], 
     [ b ,  e ,  x ,  0 ], 
     [ c ,  f ,  y ,  z ]])):

for i from 1 to 4 by 3 do for j from 1 to 4 do A[i,j]:=simplify(A[i,j]/I): od: od:
r:= A[1,1]:
Xp:=solve(t*cos(th) + r*sin(th) = A[2,1] + A[2,2]*Xp, Xp):
Yp:=solve(t*sin(th) - r*cos(th) = A[3,1] + A[3,2]*Xp + A[3,3]*Yp, Yp):
Zp:=solve(t = A[4,1] + A[4,2]*Xp + A[4,3]*Yp + A[4,4]*Zp, Zp):

g:=2:

Hyp10:=plot3d([Xp, Yp, Zp ],t=-10*abs(g)..10*abs(g),th=-Pi..Pi,color=tan,grid=[40,40]):
tC1 :=TanPoint(Q10,B,aSOLS[1]):
tC2 :=TanPoint(Q10,B,aSOLS[2]):
tC3 :=TanPointb(Q10,B,bSOLS[1]):
tC4 :=TanPointb(Q10,B,bSOLS[2]):
###################################################################
NFR:=180:

Eps:=0.08:
#Eps:=0.28:

tC1[1]:=tC1[1]+Eps*sin(th):
tC2[1]:=tC2[1]+Eps*sin(th):
tC3[1]:=tC3[1]+Eps*sin(th):
tC4[1]:=tC4[1]+Eps*sin(th):
tC1[2]:=tC1[2]+Eps*cos(th):
tC2[2]:=tC2[2]+Eps*cos(th):
tC3[2]:=tC3[2]+Eps*cos(th):
tC4[2]:=tC4[2]+Eps*cos(th):
TC1:=[tC1[1]*cos(T) + tC1[2]*sin(T), 
      tC1[2]*cos(T) - tC1[1]*sin(T), tC1[3]]:
TC2:=[tC2[1]*cos(T) + tC2[2]*sin(T), 
      tC2[2]*cos(T) - tC2[1]*sin(T), tC2[3]]:
TC3:=[tC3[1]*cos(T) + tC3[2]*sin(T), 
      tC3[2]*cos(T) - tC3[1]*sin(T), tC3[3]]:
TC4:=[tC4[1]*cos(T) + tC4[2]*sin(T), 
      tC4[2]*cos(T) - tC4[1]*sin(T), tC4[3]]:
######################################################
eps:=0.08: 
Axis:=[X*cos(ph) + eps*cos(th)*sin(ph),
       eps*cos(th)*cos(ph) - X*sin(ph), eps*sin(th)]:
s2:=sqrt(2):

plotsetup(gif,plotoutput=`10.r.gif`,plotoptions=`height=240,width=240`):
display([
animate3d([TC1[]],B=-1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[30,5],frames=NFR),
animate3d([TC1[]],B=-5..-s2,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[28,5],frames=NFR),
animate3d([TC1[]],B=s2..5.6,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[20,5],frames=NFR),
animate3d([TC2[]],B=-1..-0.2,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[28,5],frames=NFR),
animate3d([TC2[]],B=0.2..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[28,5],frames=NFR),
animate3d([TC2[]],B=-3..-s2,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[18,5],frames=NFR),
animate3d([TC2[]],B=s2..2.81,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[18,5],frames=NFR),
#  X-axis
animate3d(Axis,X=-8..8,th=-Pi..Pi,ph=-Pi..Pi,color=magenta,grid=[2,5],frames=NFR)
,animate3d([(Xp)*cos(T)+sin(T)*(Yp),
             -(Xp)*sin(T)+cos(T)*(Yp),Zp],
              t=-8*abs(g)..8*abs(g),th=-Pi..Pi,T=-Pi..Pi,
              color=tan,grid=[29,35],frames=NFR)
], insequence=false, view=[-9..10,-8..7,-5..5],
#], insequence=false, view=[-19..20,-18..17,-15..15],
  orientation=[0,60],scaling = CONSTRAINED);