# Hyperboloidoid with two sheets point A at Inf

Hyperboloid Sigma: x^2/a^2 + y^2/b^2 - z^2)/c^2 = -1
Fixed point A[r*cos(u0)*sin(v0), r*sin(u0)*sin(v0), r*cos(v0)]
when r -> inf

## Intersection ("antipodal") point Dinf when A --> Inf
in direction (cos(u0)sin(v0), sin(u0)sin(v0), cos(v0))

kill(all)$

assume(0<a, 0<b, 0<c, 1<s, 0<v0, v0<%pi); 

assume(0<u0, u0<%pi/2);

/* assume(%pi/2<u0, u0<%pi); */

/* assume(%pi<u0, u0<3*%pi/2); */

/* assume(3*%pi/2<u0, u0<2*%pi); */

Generic point A

x0: r*cos(u0)*sin(v0)$ y0: r*sin(u0)*sin(v0)$ z0: r*cos(v0)$

Generic point C on Sigma

C: matrix([xhat], [yhat], [zhat])$

Line l determined  by A and C

ec1: x - x0= lambda*(xhat - x0)$ ec2: y - y0 = lambda*(yhat - y0)$ ec3: z - z0 = lambda*(zhat - z0)$ 

Auxiliary functions

ec2/ec1$ k: rhs(%)$ y: y0 + k*(x - x0)$

ec3/ec1$ m: rhs(%)$ z: z0 + m*(x - x0)$

Quadratic polynomial p(x) = a2*x^2 + a1*x + a0

p: trigsimp(x^2/a^2 + y^2/b^2 - z^2/c^2 + 1)$ a2: ratcoef(p,x,2)$ a1: ratcoef(p,x,1)$

Appliyng Vieta's formula to get point D(x1,y1,z1):

x1: trigsimp(-a1/a2 - xhat)$ y1: y0 + k*(x1-x0)$ z1: z0 + m*(x1-x0)$

Limit when r -> infty

x1inf: trigsimp(limit(x1,r,inf ))$ y1inf: trigsimp(limit(y1,r,inf))$ z1inf: trigsimp(limit(z1,r,inf))$

Matrix LD such that LD*C = Dinf

LD: 1/denom(x1inf)*matrix(
 [(((b^2+a^2)*c^2*sin(u0)^2-b^2*c^2)*sin(v0)^2-a^2*b^2*cos(v0)^2),
   -2*a^2*c^2*cos(u0)*sin(u0)*sin(v0)^2,
  2*a^2*b^2*cos(u0)*cos(v0)*sin(v0)], 
 [-2*b^2*c^2*cos(u0)*sin(u0)*sin(v0)^2,
  (((b^2+a^2)*c^2*cos(u0)^2-a^2*c^2)*sin(v0)^2-a^2*b^2*cos(v0)^2),
  2*a^2*b^2*sin(u0)*cos(v0)*sin(v0)], 
 [-2*b^2*c^2*cos(u0)*cos(v0)*sin(v0),
  -2*a^2*c^2*sin(u0)*cos(v0)*sin(v0),
  (((b^2-a^2)*c^2*cos(u0)^2+a^2*c^2)*sin(v0)^2+a^2*b^2*cos(v0)^2)]
);

Direct calculation to verify that LD*C = Dinf

Dinf: matrix(
 [x1inf], 
 [y1inf], 
 [z1inf]
);

trigsimp(LD.C - Dinf) ;

## Locus point Einf = sC + (1-s)Dinf

xEinf: ratsimp(s*xhat + (1-s)*x1inf)$ yEinf: ratsimp(s*yhat + (1-s)*y1inf)$ zEinf: ratsimp(s*zhat + (1-s)*z1inf)$

Einf: matrix(
 [xEinf], 
 [yEinf], 
 [zEinf]
);

Matrix LE such that LE*C = Einf 

I: matrix([1,0,0],  [0,1,0],  [0,0,1])$ 

LE: s*I + (1-s)*LD;

Direct calculation to verify that LE*C = Einf

trigsimp(LE.C - Einf);

## Implicit equation of Delta_inf (using quadratic form representations)

Hyperboloid Sigma: x^2/a^2 + y^2/b^2 - z^2/c^2 = -1

kill(x,y,z)$

xs: matrix( [x],  [y],  [z],  [1]);

QSigma: matrix(
 [b^2*c^2,0,0,0], 
 [0,a^2*c^2,0,0], 
 [0,0,-a^2*b^2,0], 
 [0,0,0,a^2*b^2*c^2]
);

ec_Sigma: expand(1/(a^2*b^2*c^2)*transpose(xs).QSigma.xs) = 0;

Equation of the hyperbploid Delta_inf (image of Sigma under matrix LE)

LE1: matrix(
 [LE[1,1],LE[1,2],LE[1,3],0], 
 [LE[2,1],LE[2,2],LE[2,3],0], 
 [LE[3,1],LE[3,2],LE[3,3],0], 
 [0,0,0,1]
);  /* Matrix of the quadratic form that defines LE*/ 

QDelta: transpose(invert(LE1)).QSigma.invert(LE1)$ /* Matrix of the quadratic form that defines Delta_inf */

Implicit equation of Delta_inf: Ax^2 + Bxy + Cxz + Dy^2 + Eyz + Fz^2 + Gx + Hy + Iz + J = 0
eq_gen: transpose(xs).QDelta.xs = 0; 

QE1: trigsimp(invert(LE1))$ 

QE1t: transpose(QE1)$

temp1: trigsimp(QE1.xs)$

temp2: trigsimp(QSigma.temp1)$

temp3: trigsimp(QE1t.temp2)$

eq_Delta_temp: trigsimp(transpose(xs).temp3=0)$

For s >0 such that s <> 1/2,

d: denom(lhs(eq_Delta_temp)) ;

subst([a=5, b=4, c=3, s=20, u0=%pi/3, v0=%pi/6],d);

subst([a=5, b=4, c=3, s=20, u0=%pi/6, v0=%pi/3],d);

ratsimp(d/(s-1/2)^2); /* d/(s-1/2)^2 */

is not zero, so we can write the equation of Delta as

eq_Delta: expand(-num(lhs(eq_Delta_temp))) = 0;

subst([a=5,b=4,c=3,s=20,u0=%pi/6,v0=%pi/3], eq_Delta), numer;

### The quadric Delta: Ax^2 + Bxy + Cxz + Dy^2 + Eyz + Fz^2 + Gx + Hy + Iz + J= 0
is an hyperboloid with two sheets

A: ratsimp(ratcoef(num(lhs(eq_Delta)),x,2));

B: ratsimp(ratcoef(num(lhs(eq_Delta)),x*y,1));

C: ratsimp(ratcoef(num(lhs(eq_Delta)),x*z,1));

D: ratsimp(ratcoef(num(lhs(eq_Delta)),y,2));

E: ratsimp(ratcoef(num(lhs(eq_Delta)),y*z,1));

F: ratsimp(ratcoef(num(lhs(eq_Delta)),z,2));

G: 0;

H: 0;

I: 0;

J: (((4*a^4*b^2-4*a^2*b^4)*c^4*s^2 + (4*a^2*b^4-4*a^4*b^2)*c^4*s + (a^4*b^2-a^2*b^4)*c^4)*cos(u0)^2 -
      4*a^4*b^2*c^4*s^2 + 4*a^4*b^2*c^4*s - a^4*b^2*c^4)*sin(v0)^2+
(4*a^4*b^4*c^2*s^2-4*a^4*b^4*c^2*s+a^4*b^4*c^2)*cos(v0)^2;

M1: matrix(
 [A,B/2,C/2], 
 [B/2,D,E/2], 
 [C/2,E/2,F] 
);

M2: matrix(
 [A,B/2,C/2,G/2], 
 [B/2,D,E/2,H/2], 
 [C/2,E/2,F,I/2], 
 [G/2,H/2,I/2,J]
);

rho3: rank(M1);

rho4: rank(M2);

The properties of the quadrics reduces the kind of the locus surface Delta to be an
(a) ellipsoid
(b) hyperboloid with one sheet
(c) hyperboloid with two sheets

IF THE MATRIX LE IS NOT SINGULAR, (a) and (b) are excluded since LE^(-1)(Delta) = Sigma would be a connected set 

## Calculation of intersection of the hyperboloid Sigma with locus surface Delta_inf (see Maple worksheet)

/* trigsimp(eigenvalues(LD)); /* -1 with multiplicity 1 and 1 with multiplicity 2 */ */

/* v1: eigenvectors(LD)[2][1][1]; /* eigenvector corresponding to eigenvalue mu1 = -1 */ */

/* v2: eigenvectors(LD)[2][2][1]; /* first eigenvector corresponding to eigenvalue mu2 = 1 */ */

/* v3: eigenvectors(LD)[2][2][2]; /* second eigenvector corresponding to eigenvalue mu2 = 1 */ */

### Fixed points for LD

v1: [cos(u0)*sin(v0)/cos(v0), sin(u0)*sin(v0)/cos(v0), 1] ;

v2: [a^2*cos(v0)/(c^2*cos(u0)*sin(v0)), 0, 1];

v3: [(-(a^2*sin(v0)^2*sin(u0))) / (b^2*cos(u0)*sin(v0)^2), 1, 0];

### Hyperbolical curve gamma(a,b,c,u0,v0)

load(vect)$

w: express(v2~v3) /* cross product of v2 and v3 */;

Equation of the subspace (plane) generated by v2 and v3

kill(x,y,z)$

eq3: z = expand(-(w[1]*x + w[2]*y)/ w[3]);

Equation of Sigma

eq4: z = sqrt(c^2*(1 + x^2/a^2 + y^2/b^2 ));

eq5: rhs(eq3)^2 = rhs(eq4)^2;

eq6_1: subst([x=t], solve(eq5,y)[1]);

eq6_2: subst([x=t], solve(eq5,y)[2]);

eq7_1: subst([x=t, eq6_1], eq3);

eq7_2: subst([x=t, eq6_2], eq3);

Example (just to check the intersection curve found using GeoGebra):

subst([a=5, b=4, c=3, u0=%pi/6, v0=%pi/3],[x=t, eq6_1, eq7_1]); /* One half of the intersection curve */

subst([a=5, b=4, c=3, u0=%pi/6, v0=%pi/3],[x=t, eq6_2, eq7_2]); /* Another half of the intersection curve */

Gradient of Sigma at point D=(x,y,z)

[diff(x^2/a^2 + y^2/b^2 - z^2/c^2 +1,x), diff(x^2/a^2 + y^2/b^2 - z^2/c^2 +1,y), diff(x^2/a^2 + y^2/b^2 - z^2/c^2 +1,z)];

Gradient of Delta_inf at point D=(x,y,z)

[diff(lhs(eq_Delta),x), diff(lhs(eq_Delta),y), diff(lhs(eq_Delta),z)];

## Eigenvalues and eigenvectors of Qdelta (principal axes theorem)

aa: 5; bb: 4; cc: 3;

uu0: %pi/6; vv0: %pi/3; ss: 20;

QQ: subst([a=aa,b=bb,c=cc,u0=uu0,v0=vv0,s=ss],QDelta),numer; 

eigenpairs: eigens_by_jacobi(QQ);

lambda1: eigenpairs[1][1]; lambda2: eigenpairs[1][2]; lambda3: eigenpairs[1][3]; lambda4: eigenpairs[1][4];

w1: transpose(eigenpairs[2])[1]; w2: transpose(eigenpairs[2])[2]; w3: transpose(eigenpairs[2])[3]; w4: transpose(eigenpairs[2])[4];

semiaxis1Length:sqrt(abs(lambda4/lambda1)); semiaxis2Length:sqrt(abs(lambda4/lambda2)); semiaxis3Length:sqrt(abs(lambda4/lambda3));

vv1: subst([u0=uu0,v0=vv0],v1),numer;

anglev1w1: acos((w1[1]*vv1[1] + w1[2]*vv1[2] + w1[3]*vv1[3])/ sqrt(vv1[1]*vv1[1] + vv1[2]*vv1[2] + vv1[3]*vv1[3]));

## Locus surface when s -> Inf: hyperbolical cylinder

QDeltaInf: limit(QDelta,s,inf)$

eigenvalues(QDeltaInf)$

u1: eigenvectors(QDeltaInf)[2][1][1];

u2: eigenvectors(QDeltaInf)[2][2][1];

u3: eigenvectors(QDeltaInf)[2][3][1];

u4: eigenvectors(QDeltaInf)[2][4][1];

ec_OmegaInf: trigsimp(transpose(xs).QDeltaInf.xs = 0);

rank(submatrix(4,QDeltaInf,4)); /* For an hyperbolical cylinder should be ?*/

rank(QDeltaInf);  /* an hyperbolical cylinder should be ?*/