#2.r.maple
#
#  Draws hyperbolae from the second family
#
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)):
 [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]):

s2:=2^(1/2):

#
# This hyperboloid is very interesting:  All 12 components are real;
# the ideal factors over an extension of Z.
#
##################################################################

F:=subs(p02=a*A,p03=b*A,p12=a*B,p13=b*B,p01=0,p23=0,
    linalg[innerprod](plucker,Wedge(Quadric(-1,0,-1,0,1,0,0,1,0,-1)),plucker)):

aSOLS:=[solve(subs(b=1,A=1,F)=0,a)]:
bSOLS:=[solve(subs(a=1,A=1,F)=0,B)]: 
#########################################################
# Initial hyperbola
#
Hyp:=plot3d([s2*cos(th)+t*sin(th),
              (-s2*sin(th)+t*cos(th)-1),t]
              ,t=-8..8,th=-Pi..Pi,color=coral):
#########################################################
#Ideal2:= [x,h+y,g,f,d^2+e*y,c,b,a+2*e]:
Family2:=Quadric(-2,0,0,d,1,0,0,d^2,0,-d^2):
d:=1:
#linalg[innerprod](Coords,Family2,Coords);quit;
#F;d:=1:
#subs(p02=a*A,p03=b*A,p12=a*B,p13=b*B,p01=0,p23=0,
#    linalg[innerprod](plucker,Wedge(Family2),plucker));quit;

linalg[innerprod](Coords,Family2,Coords)-
expand(X^2+(d*Y)^2-(d*Z-1)^2-1):
simplify(subs(X=cos(th)+t*sin(th),Y=(-sin(th)+t*cos(th))/d,Z=(t+1)/d,
X^2+(d*Y)^2-(d*Z-1)^2-1)):
d:=1:
NFR:=180:
Eps:=0.08:

tC1:=TanPoint(Family2,B,aSOLS[1]):
tC2:=TanPoint(Family2,B,aSOLS[2]):
tC3:=TanPointb(Family2,bSOLS[1],b):
tC4:=TanPointb(Family2,bSOLS[2],b):
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)+Eps/2:
tC2[2]:=tC2[2]+Eps*cos(th)+Eps/2:
tC3[2]:=tC3[2]+Eps*cos(th)+Eps/2:
tC4[2]:=tC4[2]+Eps*cos(th)+Eps/2:

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


plotsetup(gif,plotoutput=`2.r.gif`,plotoptions=`height=250,width=250`):
display([
animate3d([TC1[]],B=-1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC2[]],B=-1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC3[]],b=-1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
animate3d([TC4[]],b=-1..1,th=-Pi..Pi,T=-Pi..Pi,color=blue,grid=[50,5],frames=NFR),
#  X-axis
animate3d(Axis,X=-5..5,th=-Pi..Pi,ph=-Pi..Pi,color=magenta,grid=[2,5],frames=NFR)
,animate3d([(cos(th)+t*sin(th))*cos(T)+sin(T)*((-sin(th)+t*cos(th))/d),
           -(cos(th)+t*sin(th))*sin(T)+cos(T)*((-sin(th)+t*cos(th))/d),
              (t+1)/d],t=-6.5*abs(d)..6.5*abs(d),th=-Pi..Pi,T=-Pi..Pi,
              color=tan,grid=[25,40],frames=NFR)
], insequence=false, view=[-5.1..5.1,-5.1..5.1,-4..4],
  orientation=[0,75],scaling = CONSTRAINED);


time();