#animate
#
################################################################
#
# Here, we animate the families of quadrics simultaneously
# tangent to the envelope of lines perpendicular to the x-axis and
# tangent to the sphere 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 ]]):

Coords:=linalg[vector]([1,X,Y,Z]):
#######################################################################
#
#  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]):

##################################################################
F:=subs(p02=W*Y,p03=W*Z,p12=X*Y,p13=X*Z,p01=0,p23=0,Equation(0,2,0,1)):

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
#
###########################################################
#
#  First family with the ellipses
#
#                Quadric(a,b,c,d,e,f,g,  h  ,x,  y  ) 
Family1 := evalm(Quadric(3,0,c,0,1,0,0,c^2/4,0,c^2/4)):
Ideal1  := [x,h-y,g,f,d,c^2-4*e*y,b,a-3*e]:
linalg[innerprod](Coords,Family1,Coords)-
  expand(-1+X^2+c^2/4*(Y+4/c)^2 +c^2/4*Z^2):
##############################################################
#
#
#
#               Quadric(a,b,c,d, e,f,g,  h   ,x,   y  ) 
Family2 := eval(Quadric(1,0,0,d,-1,0,0,d^2/4,0,d^2/4)):
#Ideal2  := [x,h-y,g,f,d^2+4*e*y,c,b,a+e]:
#expand(-3-X^2+d^2/4*Y^2+d^2/4*(Z+4/d)^2):
###############################################################
#
#
#
#               Quadric(a,b,c,d,e,f,g,  h  ,x,   y    ) 
Family3 := eval(Quadric(3,0,0,0,1,f,0,f^2/4,0,-3*f^2/4)):
Ideal3  := [x,3*h+y,g,4*e*y+3*f^2,d,c,b,a-3*e]:
#linalg[innerprod](Coords,Family3,Coords):
#Eq3:=expand(3 + (X+f*Y)^2 - 3*f^2/4*Y^2 - 3*f^2/4*Z^2):
#simplify(subs(X=1+2*t-2/s3*(cos(th)*(2+t)+sin(th)*s3*t),
#              Y=Ga/s3*(cos(th)*(2+t)+sin(th)*s3*t),
#              Z=Ga/s3*(-sin(th)*(2+t)+cos(th)*s3*t),
#              Ga=2/f,s3=3^(1/2),Eq3)):
################################################################
#
#               Quadric(a,b,c,d, e,f,g, h  ,x,     y    ) 
Family4 := eval(Quadric(1,0,0,0,-1,0,g,g^2/4,0,-3*g^2/4)):
#Ideal4  := [x,3*h+y,f,4*e*y-3*g^2,d,c,b,a+e]:
#linalg[innerprod](Coords,Family4,Coords):
#expand(1-(X-g*Z)^2 + g^2*Y^2/4 + g^2*Z^2/4 ):
##########################################################################
#
#   Begin Drawing
#
#########################################################################
#

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

AA := display3d(
#                   The  x - axis
spacecurve([t,0,0],t=-8..8,color=red,thickness=3,numpoints=2),
#
#   Tangent lines
#
spacecurve([1. ,t, 0],t=-10..10,color=pagr,numpoints=2),
spacecurve([-1.,t, 0],t=-10..10,color=pagr,numpoints=2),
spacecurve([.97,t, subs(X=.97,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.86,t, subs(X=.86,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.66,t, subs(X=.66,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.42,t, subs(X=.42,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.18,t, subs(X=.18,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([0,t, subs(X=0,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.18,t, subs(X=-.18,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.42,t, subs(X=-.42,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.66,t, subs(X=-.66,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.86,t, subs(X=-.86,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.97,t, subs(X=-.97,t*ZSOLS[1])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.97,t, subs(X=.97,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.86,t, subs(X=.86,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.66,t, subs(X=.66,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.42,t, subs(X=.42,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([.18,t, subs(X=.18,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([0,t, subs(X=0,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.18,t, subs(X=-.18,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.42,t, subs(X=-.42,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.66,t, subs(X=-.66,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.86,t, subs(X=-.86,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2),
spacecurve([-.97,t, subs(X=-.97,t*ZSOLS[2])],t=-10..10,color=pagr,numpoints=2)):


###########################################################################
#
#   The Original Sphere
#
Sphere := plot3d([sin(ph)*cos(th),-2-sin(ph)*sin(th),-cos(ph)],
            ph=0..Pi,th=-Pi..Pi,grid=[15,15],color=tan):

##########################################################
#
#   Begin animation
#
#
#
NFR:=250:  # Number of Frames

############################################################
#   Elipsoids from the first family
#
Sc:=2.5:
C:=c -> -2/(Sc*(sin(c)+1/2) + sqrt(1+Sc^2*(sin(c)+1/2)^2)):

Ellipse1:=animate3d([sin(ph)*cos(th),-2*C(c)+C(c)*sin(ph)*sin(th),C(c)*cos(ph)],
             ph=0..Pi,th=-Pi..Pi,c=-Pi..Pi,grid=[30,30],color=magenta,frames=NFR):

plotsetup(gif,plotoutput=`I-movie.gif`,plotoptions=`height=500,width=500`):
display3d([Sphere,Ellipse1,AA],insequence=false, 
       view=[-9..9,-6..12,-5..5],orientation=[200,40]);

################################################################
#   Hyperboloids from the second family
#
#  The hyperboloid at 2/d = Be
Sd:=2:
Be:=d -> 2/(Sd*(sin(d)+1/2) + sqrt(1+Sd^2*(sin(d)+1/2)^2)):

plotsetup(gif,plotoutput=`II-movie.gif`,plotoptions=`height=500,width=500`):
Hyp2:=animate3d([1+2*t,Be(d)*(cos(th)*(2+t)+sin(th)*3^(1/2)*t),
               -2*Be(d)+Be(d)*(-sin(th)*(2+t)+cos(th)*3^(1/2)*t)]
               ,t=-3..3,th=-Pi..Pi,d=-Pi..Pi,grid=[30,30],color=magenta,frames=NFR):
display3d([Sphere,Hyp2,AA],insequence=false, 
       view=[-8..8,-6..12,-5..5],orientation=[220,40]);

################################################################
#   Hyperboloids from the third family
#
#  The hyperboloid at 2/f = Ga
Sf:=2:
s3:=3^(1/2):
Ga:=f -> 2/(Sf*(sin(f)+1/2) + sqrt(1+Sf^2*(sin(f)+1/2)^2)):

plotsetup(gif,plotoutput=`III-movie.gif`,plotoptions=`height=500,width=500`):
Hyp3:=animate3d([1+2*t-2/s3*(cos(th)*(2+t)+sin(th)*s3*t),
             Ga(f)/s3*(cos(th)*(2+t)+sin(th)*s3*t),
             Ga(f)/s3*(-sin(th)*(2+t)+cos(th)*s3*t)]
             ,t=-15..15,th=-Pi..Pi,f=-Pi..Pi,grid=[40,35],color=magenta,frames=NFR):

display3d([Sphere,Hyp3,AA],insequence=false, 
       view=[-8..8,-6..12,-5..5],orientation=[200,40]);

################################################################
#   Hyperboloids from the fourth family
#
#  The hyperboloid at 
Sg:=2.5:
G:=g -> Sg*(sin(g)+1/2) + sqrt(1+Sg^2*(sin(g)+1/2)^2):

Hyp4a:=animate3d([G(g)*R*sin(th)+sqrt(1+G(g)^2/4*R^2), R*cos(th), R*sin(th)]
                  ,R=0..2,th=-Pi..Pi,g=-Pi..Pi,grid=[30,40],color=magenta,frames=NFR):
Hyp4b:=animate3d([G(g)*R*sin(th)-sqrt(1+G(g)^2/4*R^2), R*cos(th), R*sin(th)]
                 ,R=0..8,th=-Pi..Pi,g=-Pi..Pi,grid=[20,30],color=magenta,frames=NFR):
#
plotsetup(gif,plotoutput=`IV-movie.gif`,plotoptions=`height=500,width=500`):
display3d([Sphere,Hyp4a,Hyp4b,AA],insequence=false, 
       view=[-8..8,-6..12,-5..5],orientation=[60,100]);