#7.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)]:

#F;g:=1:s:=-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):
#primpart(subs(p02=aA*Aa,p03=bB*Aa,p12=aA*Bb,p13=bB*Bb,p01=0,p23=0,
#    linalg[innerprod](plucker,Wedge(Family7),plucker)));quit;


#########################################################
# 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):
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):
                     
#   Yet another check
#Co:=linalg[vector]([1,Xp,Yp,Zp]):
#evalf(linalg[innerprod](Co,Family7,Co));quit;

g:=2:

Hyp7:=plot3d([Xp, Yp, Zp ],t=-8*abs(g)..8*abs(g),th=-Pi..Pi,color=tan,grid=[40,50]):
TC1 :=TanPoint(Family7,B,aSOLS[1]):
TC2 :=TanPoint(Family7,B,aSOLS[2]):
###################################################################
macro(pagr=COLOR(RGB, .5, 1, .5)):

AA:=display3d(
spacecurve([t,0,0],t=-8..8,color=magenta,thickness=5),
Hyp,  Hyp7,
spacecurve([TC1[]],B=-1..1,color=blue,thickness=3,numpoints=100),
spacecurve([TC1[]],B=sqrt(2)..5.5,color=blue,thickness=3,numpoints=100),
spacecurve([TC1[]],B=-5.5..-sqrt(2),color=blue,thickness=3,numpoints=100),
spacecurve([TC2[]],B=-1..-0.1,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC2[]],B=0.1..1,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC2[]],B=sqrt(2)..5.5,color=blue,thickness=3,numpoints=100),
spacecurve([TC2[]],B=-6..-sqrt(2),color=blue,thickness=3,numpoints=100),
#
#  The tangent lines
#
spacecurve([1.  ,t*subs(B=1   ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.99,t*subs(B=0.99,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.96,t*subs(B=0.96,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.90,t*subs(B=0.90,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.79,t*subs(B=0.79,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.64,t*subs(B=0.64,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.46,t*subs(B=0.46,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.27,t*subs(B=0.27,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.09,t*subs(B=0.09,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.09,t*subs(B=-0.09,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.27,t*subs(B=-0.27,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.46,t*subs(B=-0.46,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.64,t*subs(B=-0.64,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.79,t*subs(B=-0.79,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.90,t*subs(B=-0.90,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.96,t*subs(B=-0.96,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.99,t*subs(B=-0.99,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.  ,t*subs(B=-1   ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
#
spacecurve([1.  ,t*subs(B=1   ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.99,t*subs(B=0.99,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.96,t*subs(B=0.96,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.90,t*subs(B=0.90,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.79,t*subs(B=0.79,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.64,t*subs(B=0.64,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.46,t*subs(B=0.46,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.27,t*subs(B=0.27,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([0.09,t*subs(B=0.09,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.09,t*subs(B=-0.09,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.27,t*subs(B=-0.27,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.46,t*subs(B=-0.46,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.64,t*subs(B=-0.64,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.79,t*subs(B=-0.79,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.90,t*subs(B=-0.90,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.96,t*subs(B=-0.96,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-0.99,t*subs(B=-0.99,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.  ,t*subs(B=-1   ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
#
#  Second set of tangent lines
#
spacecurve([4 ,t*subs(B=4 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([3.7 ,t*subs(B=3.7 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([3.4 ,t*subs(B=3.4 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([3.1 ,t*subs(B=3.1 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([2.8 ,t*subs(B=2.8 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([2.5 ,t*subs(B=2.5 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([2.21,t*subs(B=2.21,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.99,t*subs(B=1.99,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.77,t*subs(B=1.77,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.58,t*subs(B=1.58,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.48,t*subs(B=1.48,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.435,t*subs(B=1.435,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.4165,t*subs(B=1.4165,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.4165,t*subs(B=1.4165,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.435,t*subs(B=1.435,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.48,t*subs(B=1.48,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.58,t*subs(B=1.58,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.77,t*subs(B=1.77,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([1.99,t*subs(B=1.99,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([2.21,t*subs(B=2.21,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([2.5 ,t*subs(B=2.5 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([2.8 ,t*subs(B=2.8 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([3.1 ,t*subs(B=3.1 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([3.4 ,t*subs(B=3.4 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([3.7 ,t*subs(B=3.7 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([4 ,t*subs(B=4 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
#
spacecurve([-4 ,t*subs(B=4 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-3.7 ,t*subs(B=3.7 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-3.4 ,t*subs(B=3.4 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-3.1 ,t*subs(B=3.1 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-2.8 ,t*subs(B=2.8 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-2.5 ,t*subs(B=2.5 ,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-2.21,t*subs(B=2.21,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.99,t*subs(B=1.99,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.77,t*subs(B=1.77,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.58,t*subs(B=1.58,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.48,t*subs(B=1.48,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.435,t*subs(B=1.435,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.4165,t*subs(B=1.4165,aSOLS[2]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.4165,t*subs(B=1.4165,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.435,t*subs(B=1.435,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.48,t*subs(B=1.48,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.58,t*subs(B=1.58,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.77,t*subs(B=1.77,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.99,t*subs(B=1.99,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-2.21,t*subs(B=2.21,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-2.5 ,t*subs(B=2.5 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-2.8 ,t*subs(B=2.8 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-3.1 ,t*subs(B=3.1 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-3.4 ,t*subs(B=3.4 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-3.7 ,t*subs(B=3.7 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
spacecurve([-4 ,t*subs(B=4 ,aSOLS[1]),t],t=-10..10,color=pagr,numpoints=2),
#
view=[-6..6,-9..5,-5..5]):

plotsetup(gif,plotoutput=`7.gif`,plotoptions=`height=600,width=600`):
display(AA,orientation=[75,60],scaling=CONSTRAINED);


AA:=display3d(
spacecurve([t,0,0],t=-8..8,color=magenta,thickness=5),
Hyp7,
spacecurve([TC1[]],B=-1..1,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC1[]],B=sqrt(2)..5,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC1[]],B=-5..-sqrt(2),color=BLUE,thickness=3,numpoints=100),
spacecurve([TC2[]],B=-1..-0.1,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC2[]],B=0.1..1,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC2[]],B=sqrt(2)..5,color=BLUE,thickness=3,numpoints=100),
spacecurve([TC2[]],B=-6..-sqrt(2),color=BLUE,thickness=3,numpoints=100),
view=[-7..7,-9..5,-5..5]):
plotsetup(gif,plotoutput=`sm.gif`,plotoptions=`height=130,width=160`):
display(AA,orientation=[75,60],scaling=CONSTRAINED);