# Ellipsoid point A at Inf

Ellipsoid 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, 0<s, 0<v0, v0<%pi);

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

/* assume(%pi/2<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: ratcoeff(p,x,2)$ a1: ratcoeff(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(
 [(a^2*c^2*sin(u0)^2 - b^2*c^2*cos(u0)^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*c^2*cos(u0)^2 - a^2*c^2*sin(u0)^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),
  (a^2*c^2*sin(u0)^2 + b^2*c^2*cos(u0)^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)

Ellipsoid 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]
);

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

Equation of the ellipsoid 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]
);  /* Extended matrix associate with LE*/ 

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

eq_Delta: trigsimp(transpose(xs).QDelta.xs = 0); /* Implicit equation od Delta_inf: Ax^2 + Bxy + Cxz + Dy^2 + Eyz + Fz^2 + Gx + Hy + Iz + J = 0 */

## Calculation of intersection of the ellipsoid Sigma with locus ellipsoid Delta_inf

### Fixed points for LD

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 */

Singular values of LD

sqrt(eigenvalues(transpose(LD).LD));

### Elliptical 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 = -(w[1]*x + w[2]*y)/ w[3];

eq4: z = sqrt(c^2*(1 - x^2/a^2 - y^2/b^2 )); /* Equation of Sigma */

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/4],[x=t, eq6_1, eq7_1]); /* One half of the intersection curve */

subst([a=5, b=4, c=3, u0=%pi/6, v0=%pi/4],[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/3; vv0: %pi/4; 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(-lambda4/lambda1); semiaxis2Length:sqrt(-lambda4/lambda2); semiaxis3Length:sqrt(-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: elliptical cylinder

QDeltaInf: limit(QDelta,s,inf);

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

rank(submatrix(4,QDeltaInf,4)); /* For an elliptical cylinder should be 2*/

rank(QDeltaInf);  /* For an elliptical cylinder should be 3*/