#procedures
#
# Frank Sottile
# 25 November 2004
# Mexico City
#
######################################################################
with(plots):
with(plottools):
plotsetup(x11,plotoptions=`width=3in,height=2.5in`):
######################################################################
#
#    This file contains all the procedures for the construction and 
# counting of the 62 lines.
#
######################################################################
#
# ParameterizeLine(Line,t)  Parameterizes the line with t=0 the first point
#                             and t=1 the second.   Works in 2D and 3D.
# Meet(Line1,Line2)         Returns the equation for Line1 to meet Line2
# MakeRegulus(Line1, Line2, Line3, t)   Gives the linked parametrizations of
#                           Lines 1 and 2 for the regulus defined by the 3 lines
#                           This parametrizes a hyperboloid
# PlaneEquation(Line, pt)   Returns the equation of the plane through
#                           Line and pt.  
# traceHyperboloid(Line1, Line2, Line3, Eqn)    Finds the 3-dimensional trace of
#                           the hyperboloid through Line1--3 in the plane
#                           defined by Eqn. 
# TangentLine (C,t0)        Returns the tangent line to a space curve C which 
#                           is parametrized by t.
# ParameterizePlane(Line, pt, u, v)        This is a plane through Line  
#                           and pt, with parameters u and v.  pt=(0,-1)
# CurveInPlane(Plane,C)     Writes a point C in local coordinates Plane
# LineInPlane(Plane,Line)   Writes Line in local coordinates Plane
# PerturbLine(Line,pt,pivot,e)       Returns a perturbed version of Line, 
#                           which meets it at Line(pivot), but with second point
#                           Line(1-pivot)+e*pt. 
# perturbTan(tan,pt,e,f)    Moves the first point of tan the fraction e
#                           further from pt, and the secont point f closer
# shear(tan,pt)             Applies a shear to pt in local coordinates for the plane
#                           so that tan is vertical.
# EXpand(P,P0,t0,ex,ey,et)  Scales and reparametrizes a point P in a plane to enable
#                           better viewing:  Basically a workaround for an apparent
#                           MAPLE bug that prevents enlarging a picture by too much.
# Triangle(l,k,m)           returns the three vertices of the triangle formes by the 
#                           lines l,k,m 
# numberOfTangents(T)       Computes the number of tangents of a given 4 triangles T.  
# CheckMeet(L)              Checks to see if the lines meeting the four lines stored in L
#                           meet them in the four intervals stored in L
#
################################################################################
#
ParameterizeLine := proc(Line,t)
#
#  Given a line by two points [ [x1,y1,z1], [x2,y2,z2] ]
#  Parameterizes the line with t=0 the first point and t=1 the second
#  Returns [x(t),y(t),z(t)].  Also works in 2D.
#
 if nops(Line[1])=3 then
   return([Line[1][1]+t*(Line[2][1]-Line[1][1])
          ,Line[1][2]+t*(Line[2][2]-Line[1][2])
          ,Line[1][3]+t*(Line[2][3]-Line[1][3])]):
  else
   return([Line[1][1]+t*(Line[2][1]-Line[1][1])
          ,Line[1][2]+t*(Line[2][2]-Line[1][2])]):
 end if:
end proc:
################################################################################
#
Meet := proc(Line1,Line2)
#
#  Returns the equation for Line1 to meet Line2
#
 return(linalg[det](linalg[augment](linalg[matrix](4,1,1),
              linalg[matrix]([Line1[],Line2[]])))):
end proc:
################################################################################
#
MakeRegulus := proc(Line1,Line2, Line3,t)
#
#  This gives parametrized points on Line1 and Line2 that are collinear with
#  Line 3.  This defines the regulus of lines meeting all three.  It also
#  serves to parametrize the hyperboloid, via the command
#    ParameterizeLine(MakeRegulus(Line1, Line2, Line3,t), s):
#  Returns []
#
#
 local s:
 return( [ParameterizeLine(Line1,t),
       subs(s=solve(Meet([ParameterizeLine(Line1,t),ParameterizeLine(Line2,s)],
          Line3)=0,s),ParameterizeLine(Line2,s))])
end proc:
###################################################################################
#
PlaneEquation:=proc(Line, pt)
#
#  Returns the 1-form defining the plane containing the Line and the pt
#
 local A,a,b,c,d,eq,var:
 A:=solve({a*Line[1][1]+b*Line[1][2]+c*Line[1][3]+d=0, 
          a*Line[2][1]+b*Line[2][2]+c*Line[2][3]+d=0, 
          a*pt[1]+b*pt[2]+c*pt[3]+d=0},{a,b,c,d}):
 eq:=subs(op(A),a*x+b*y+c*z+d):
 var:=op(indets(eq) minus ({x,y,z} union indets(pt))):
 return(numer(subs(var=1,eq))):
end proc:
################################################################################
#
traceHyperboloid:=proc(Line1, Line2, Line3, Eqn)
#
#  Find the 3-dimensional trace of the hyperboloid through Line1--3
#  in the plane defined by Eqn.
#
 local R, S:
 R:=ParameterizeLine(MakeRegulus(Line1, Line2, Line3, t), s);
 S:=solve(subs(x=R[1], y=R[2], z=R[3], Eqn)=0,s);
 return(simplify( subs(s=S,R) )):
end proc:
######################################################################
TangentLine := proc(C,t0)
#
# Returns the tangent line to a space curve C, which is a function of t
#
 local A, B:
 A:=subs(t=t0,C):
 B:=subs(t=t0,diff(C,t)):
 return([A,[A[1]+B[1], A[2]+B[2], A[3]+B[3]]]):
end proc:
################################################################################
#
ParameterizePlane := proc(Line, pt, u, v)
#
#  This is a plane through Line and pt, 
#  the parameters are u and v.
#
 local U, V:
 U:=ParameterizeLine(Line,u);
 V:=[Line[1][1]-pt[1],Line[1][2]-pt[2],Line[1][3]-pt[3]]:
 return([U[1]+v*V[1],U[2]+v*V[2],U[3]+v*V[3]]):
end proc:
################################################################################
#
CurveInPlane:=proc(Plane,C)
#
#  Intersect C with Plane and writes in local coordinates
#
 local A, U, V:
 A:=solve({Plane[1]=C[1], Plane[2]=C[2], Plane[3]=C[3]},{u,v}):
 if member(u,indets(op(1,A))) then 
   U:=op(1,A): V:=op(2,A):
   else
  V:=op(1,A): U:=op(2,A):
 end if:
 return([op(2,U),op(2,V)]):
end proc:
###############################################################################
LineInPlane:=proc(Plane,Line)
#
# Returns a line in a plane
#
 return([CurveInPlane(Plane,Line[1]), CurveInPlane(Plane,Line[2])])
end proc:
###############################################################################
#
PerturbLine:=proc(Line,pt,pivot,e)
#
#    Returns a perturbed version of Line, which meets it at Line(pivot),
# but with second point Line(1)+e*pt.
#
 local A,B:
 A:=ParameterizeLine(Line,pivot):
 B:=ParameterizeLine([Line[1],pt],e):  #Line(1)+e*(pt-Line(1))
 return([B, ParameterizeLine([B,A],1/pivot)]):
end proc:
###############################################################################
#
perturbTan:=proc(tan,pt,e,f)
#
#   perturbs the tangent line to (hopefully) a secant line
#
 return([[tan[1][1]+e*(-pt[1]+tan[1][1]),
          tan[1][2]+e*(-pt[2]+tan[1][2]),
          tan[1][3]+e*(-pt[3]+tan[1][3])],
         [tan[2][1]+f*(pt[1]-tan[2][1]),
          tan[2][2]+f*(pt[2]-tan[2][2]),
          tan[2][3]+f*(pt[3]-tan[2][3])]]):
end proc:
##############################################################################
#
shear:=proc(tan,pt)
#
#   This applies a shear in local coordinates for a plane in order to make the
# short side of the triangle vertical.  (Makes it easier to view.)
#
 local v:
 v:=diff(tan,t):
 return([pt[1]-v[1]*pt[2]/v[2], pt[2]]):
end proc:
##############################################################################

EXpand:=proc(P,P0,t0,ex,ey,et)
#
# This scales a point in the plane P about the point P0 by the factors ex and ey 
# in the x- and y-directions, while also scaling t about the point t0 by the 
# factor et.   It returns P centered at the origin.
 if type(t0,symbol) then
   return([P0[1]+(P[1]-P0[1])*ex,P0[2]+(P[2]-P0[2])*ey]):
  else
   return(subs(t=t0 + (t-t0)/et, [P0[1]+(P[1]-P0[1])*ex,P0[2]+(P[2]-P0[2])*ey])):
 end if:
end proc:

##############################################################################
Triangle:=proc(l,k,m)
#
#   returns the three vertices of the triangle formes by the lines l,k,m
#
 local A,B,C,s,t,u,P,Q,R:
 A:= ParameterizeLine(l,t):
 B:= ParameterizeLine(k,s):
 C:= ParameterizeLine(m,u):
 P:=solve({A[1]=B[1],A[2]=B[2],A[3]=B[3]},{s,t}):
 Q:=solve({A[1]=C[1],A[2]=C[2],A[3]=C[3]},{u,t}):
 R:=solve({C[1]=B[1],C[2]=B[2],C[3]=B[3]},{u,s}):
 return([subs(op(P),A), subs(op(Q),A), subs(op(R),B)]):
end proc:
##################################################################################
numberOfTangents := proc(T)
#
#  Procedure to compute the number of transversals of a given
#  triangle T.  
#
#################################################################
 local file, NumberOfTransversals, Line, L, Regulus, 
       Equation, Points, j, P, M, N, SOL:
 ################################################################
 NumberOfTransversals:=0:
 ################################################################
 for Line[1] in combinat[choose](T[1],2) do 
  for  Line[2] in combinat[choose](T[2],2) do 
   L[2]:=ParameterizeLine(Line[2],t):
   for  Line[3] in combinat[choose](T[3],2) do 
    L[3]:=ParameterizeLine(Line[3],t):
    for  Line[4] in combinat[choose](T[4],2) do 
     L[4]:=ParameterizeLine(Line[4],t):

     Regulus:=MakeRegulus(Line[1], Line[2], Line[3], t):
     Equation:=numer(Meet(Line[4],Regulus)):
     Points:=[solve(Equation=0,t)]:
#
     if not type(Points[1],nonreal) then 
      for j from 1 to 2 do 
       P:=Points[j]:
       if 0<=evalf(P) and evalf(P)<=1 then 
        M[j]:=subs(t=P,Regulus):
#       
# By the parametrization of the regulus, the second coordinate is on Line[2]
#       
        P:=solve(L[2][1]=M[j][2][1],t):
        if 0<=evalf(P) and evalf(P)<=1 then
#        
# The lines 3 and 4 are slightly different
         N[j]:=ParameterizeLine(M[j],s):
         SOL:=solve({L[3][1]=N[j][1],L[3][2]=N[j][2],L[3][3]=N[j][3]},{s,t}):
         P:=subs(op(SOL),t):
         if 0<=evalf(P) and evalf(P)<=1 then
          SOL:=solve({L[4][1]=N[j][1],L[4][2]=N[j][2],L[4][3]=N[j][3]},{s,t}):
          P:=subs(op(SOL),t):
          if 0<=evalf(P) and evalf(P)<=1 then
           NumberOfTransversals:=NumberOfTransversals+1: 
          end if:
         end if:
        end if:
       end if:
      end do:
     end if:
##################################################################
    end do:
   end do:
  end do:
 end do:
 return(NumberOfTransversals):
end proc:
##################################################################
####################################################################################
CheckMeet:=proc(L,N)
#
# Checks to see if the lines meeting the four lines stored in L
# meet them in the four intervals stored in L
#
 local P, S, St, Ls, DI:
 DI:=30:
 # Find the points where the common lines meet Line a
 #
 P:=[solve(Meet(MakeRegulus(L[1], L[2], L[3], t), L[4])=0,t)]:
 #
 # Check to see if these points are real
 #
 if type(P[1][1],nonreal) then
   lprint("Only Complex Solutions");
  else
   #
   #   Computes the two common lines
   #
   S:=[MakeRegulus(L[1], L[2], L[3], P[1]), MakeRegulus(L[1], L[2], L[3], P[2])]:
   #
   #   Computes the points where the two common lines meet the other
   # 3 lines.  First it parametrizes all of them, and then it solves
   # the equations.
   #
   St:=[ParameterizeLine(S[1],t),ParameterizeLine(S[2],t)]:
   Ls:=[ParameterizeLine(L[2],s),ParameterizeLine(L[3],s),ParameterizeLine(L[4],s)]:
   #
   #  First line
   #
   P:=[[P[1],
     subs(solve({St[1][1]=Ls[1][1],St[1][2]=Ls[1][2],St[1][3]=Ls[1][3]},{s,t}),s),
     subs(solve({St[1][1]=Ls[2][1],St[1][2]=Ls[2][2],St[1][3]=Ls[2][3]},{s,t}),s),
     subs(solve({St[1][1]=Ls[3][1],St[1][2]=Ls[3][2],St[1][3]=Ls[3][3]},{s,t}),s)],
    [P[2],
     subs(solve({St[2][1]=Ls[1][1],St[2][2]=Ls[1][2],St[2][3]=Ls[1][3]},{s,t}),s),
     subs(solve({St[2][1]=Ls[2][1],St[2][2]=Ls[2][2],St[2][3]=Ls[2][3]},{s,t}),s),
     subs(solve({St[2][1]=Ls[3][1],St[2][2]=Ls[3][2],St[2][3]=Ls[3][3]},{s,t}),s)]]:

   if N=1 then 
    lprint(evalf(P[1],20));
    lprint(evalf(P[2],20));
   end if:

   if (evalb( evalf(P[1][1],DI) >=evalf(P[1][1]^2,DI) ) and 
       evalb( evalf(P[1][2],DI) >=evalf(P[1][2]^2,DI) ) and 
       evalb( evalf(P[1][3],DI) >=evalf(P[1][3]^2,DI) ) and 
       evalb( evalf(P[1][4],DI) >=evalf(P[1][4]^2,DI) )) then 
    lprint("First line meets 4 triangles");
   end if:
   if (evalb( evalf(P[2][1],DI) >=evalf(P[2][1]^2,DI) ) and 
       evalb( evalf(P[2][2],DI) >=evalf(P[2][2]^2,DI) ) and 
       evalb( evalf(P[2][3],DI) >=evalf(P[2][3]^2,DI) ) and 
       evalb( evalf(P[2][4],DI) >=evalf(P[2][4]^2,DI) )) then 
    lprint("Seecond line meets 4 triangles");
   end if:
  end if:      #   Ends the loop for reality of solutions
end proc: