#7.r.maple
#
#   Draws the hyperbolae from the first half of the sixth 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)]:

#########################################################
# Initial hyperbola
#
Hyp:=plot3d([2^(1/2)*cos(th)+t*sin(th),
            (-2^(1/2)*sin(th)+t*cos(th)-1),t]
            ,t=-8..8,th=-Pi..Pi,color=coral):

######################################################################
#
s:=-2^(1/2):   # The `sister quadric to this has s:=2^(1/2):
s2:=2^(1/2):   # The `sister quadric to this has s:=2^(1/2):
Family7:=Quadric(g^2*(2+s)/(1+s), 1/2*g^2/(1+s), g*(2+s)/(1+s), 
           g, -g^2, -g, g, -4-2*s, 1, 2+2*s):
##################################################################
#     [ 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 2 to 4 by 2 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 = A[3,1] + A[3,2]*Xp + A[3,3]*Yp, Yp):
Zp:=solve(t*sin(th) - r*cos(th) =  A[4,1] + A[4,2]*Xp + A[4,3]*Yp + A[4,4]*Zp, Zp):
                     
g:=2:

###################################################################
NFR:=180:
Eps:=0.08:
tC1:=TanPoint(Family7,B,aSOLS[1]):
tC2:=TanPoint(Family7,B,aSOLS[2]):

tC1[1]:=tC1[1]+Eps*sin(th):
tC2[1]:=tC2[1]+Eps*sin(th):
tC1[2]:=tC1[2]+Eps*cos(th):
tC2[2]:=tC2[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]]:
######################################################
eps:=0.08: 
Axis:=[X*cos(ph) + eps*cos(th)*sin(ph),
       eps*cos(th)*cos(ph) - X*sin(ph), eps*sin(th)]:

plotsetup(gif,plotoutput=`7.r.gif`,plotoptions=`height=220,width=300`):
display([
animate3d([TC1[]],B=-1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC1[]],B=-5..-s2,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC1[]],B=s2..5,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC2[]],B=-1..-0.1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC2[]],B=0.1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC2[]],B=-5..-s2,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC2[]],B=s2..5,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,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=[30,40],frames=NFR)
], insequence=false, view=[-10..10,-10..10,-5..5],
  orientation=[0,75],scaling = CONSTRAINED);