#rotate.maple
#
################################################################
#
# Here, we investigate the families of quadrics simultaneously
# tangent to the envelope of lines perpendicular to the x-axis and
# tangent to the hyperbola of radius 1 centered at (0,2,0).
# (These lines are transversal to both the x-axis and the y-z line
# at infinity.)
#
interface(quiet=true):
with(plots):
plotsetup(x11,plotoptions=`width=3in,height=2.5in`):
########################################################
#
#  Computes quadratic form of sphere
#
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 ]]):
########################################################
#
#  Computes second exterior power of a 4 by 4 matrix
#
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 in Pl\"ucker coordinates for lines tangent to the sphere
#   with center (a,b,c) and radius r
#
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:

#
#   Matrix for a general quadratic form
#
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 ]]):

#######################################################################
#
#  This procedure computes the point of tangency of a given tangent line.
#  We are working in the local coordinates for lines perpendicular to 
#  x-axis, where A = b = 1.
#
#  a is the function a(B) stored in BSOLS
#
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:
##################################################################

Coords:=linalg[vector]([1,X,Y,Z]):

##################################################################
plucker:=linalg[vector]([p01,p02,p03,p12,p13,p23]):
Hyp:=linalg[innerprod](plucker,Wedge(Quadric(3,0,-2,0,1,0,0,1,0,-1)),plucker):
F:=subs(p02=W*Y,p03=W*Z,p12=X*Y,p13=X*Z,p01=0,p23=0,Hyp):

ZSOLS:=[solve(subs(Y=1,W=1,F)=0,Z)]:
YSOLS:=[solve(subs(Z=1,W=1,F)=0,Y)]:

############################################################
#
#  We now draw some from each of the four real families
#
###########################################################
#
#  Hyperboloid of one sheet
#
Ideal1:= [x, h+y, g, f, d, c^2+4*e*y, b, a-3*e]:
#        e:=1    a,b,c,d,e,f,g,  h  ,x,   y
Family1:=Quadric(3,0,c,0,1,0,0,c^2/4,0,-c^2/4):
linalg[innerprod](Coords,Family1,Coords)-
expand(X^2+(c*Y+4)^2/4 - Z^2*c^2/4-1):
##########################################################################
#
#   Begin Drawing
#
#########################################################################
#
#   Original Hyperboloid
#
hyp:=[cos(th)+t*sin(th), -sin(th)+t*cos(th)+2, t]:
Hyp:=[hyp[1]*cos(ph) + hyp[2]*sin(ph), 
      hyp[2]*cos(ph) - hyp[1]*sin(ph), hyp[3]]:

#  Since I cannot animate a space curve, the tangents will be little cylinders
#
#
Eps:=0.05:  bot:=1/sqrt(1+ZSOLS[1]^2):
s1:=[X+Eps*cos(th),t + Eps*sin(th)*bot*ZSOLS[1], 
                      t*ZSOLS[1] +  Eps*sin(th)*bot]:
S1:=[s1[1]*cos(ph) + s1[2]*sin(ph), 
     s1[2]*cos(ph) - s1[1]*sin(ph), s1[3]]:
s2:=[X+Eps*cos(th),t + Eps*sin(th)*bot*ZSOLS[2], 
                      t*ZSOLS[2] +  Eps*sin(th)*bot]:
S2:=[s2[1]*cos(ph) + s2[2]*sin(ph), 
     s2[2]*cos(ph) - s2[1]*sin(ph), s2[3]]:
#Horizontal tangent
Htp:=[(Eps*cos(th)+1)*cos(ph) + t*sin(ph), 
      t*cos(ph) - (Eps*cos(th)+1)*sin(ph), Eps*sin(th)]:
Htm:=[(Eps*cos(th)-1)*cos(ph) + t*sin(ph), 
      t*cos(ph) - (Eps*cos(th)-1)*sin(ph), Eps*sin(th)]:


# The curves of tangenct will be just curves
c:=-2:
tC1:=TanPoint(Family1,X,YSOLS[1]):
tC2:=TanPoint(Family1,X,YSOLS[2]):
tC1[1]:=tC1[1]+Eps*sin(th):
tC2[1]:=tC2[1]+Eps*sin(th):
tC1[2]:=tC1[2]+Eps*cos(th)+Eps/2:
tC2[2]:=tC2[2]+Eps*cos(th)+Eps/2:

Xs:=fsolve(tC1[3]=3.9,X):

TC1:=[tC1[1]*cos(ph) + tC1[2]*sin(ph), 
      tC1[2]*cos(ph) - tC1[1]*sin(ph), tC1[3]]:
TC2:=[tC2[1]*cos(ph) + tC2[2]*sin(ph), 
      tC2[2]*cos(ph) - tC2[1]*sin(ph), 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)]:

macro(pacr=COLOR(RGB, 1, .8, 0)):
macro(pagr=COLOR(RGB, .5, 1, .5)):

NFR:=250:

plotsetup(gif,plotoutput=`oneHyp.rot.gif`,plotoptions=`height=250,width=250`):

display([
#  X-axis
animate3d(Axis,X=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=red,grid=[2,5],frames=NFR),

animate3d(Hyp,t=-3.9..3.9,th=-Pi..Pi,ph=-Pi..Pi,color=coral,grid=[15,30],frames=NFR),
animate3d([TC1[]],X=-Xs..-1.00001,th=-Pi..Pi,ph=-Pi..Pi,color=blue,grid=[100,5],frames=NFR),
animate3d([TC2[]],X=-Xs..-1.00001,th=-Pi..Pi,ph=-Pi..Pi,color=blue,grid=[100,5],frames=NFR),
animate3d([TC1[]],X=1.0001..Xs,th=-Pi..Pi,ph=-Pi..Pi,color=blue,grid=[100,5],frames=NFR),
animate3d([TC2[]],X=1.0001..Xs,th=-Pi..Pi,ph=-Pi..Pi,color=blue,grid=[100,5],frames=NFR),
#
#  Common transversals and tangents
#
animate3d(Htp,t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.021,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.021,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.11,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.11,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.34,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.34,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.7,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.7,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=2.1,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-2.1,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=2.75,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-2.75,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=3.5,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-3.5,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=4.4,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-4.4,S1),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
#
animate3d(Htm,t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.021,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.021,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.11,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.11,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.34,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.34,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=1.7,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-1.7,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=2.1,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-2.1,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=2.75,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-2.75,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=3.5,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-3.5,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=4.4,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR),
animate3d(subs(X=-4.4,S2),t=-6..6,th=-Pi..Pi,ph=-Pi..Pi,color=pagr,grid=[2,5],frames=NFR)
#
], insequence=false, view=[-6..6,-6..6,-4..4],orientation=[0,60]);


time();