AOC 2024 Day 24 Part 2 via matrix and vector math

day24_via_vectors.m

% xyz = readcsv('day24.csv');

n1 = 100;
n2 = 200;
n3 = 300;

r1 = xyz(1:3,n1);
r2 = xyz(1:3,n2);
r3 = xyz(1:3,n3);
v1 = xyz(4:6,n1);
v2 = xyz(4:6,n2);
v3 = xyz(4:6,n3);

% vector and matrix methods of determining the initial position (r0 = [x0]) and velocity (v0 = [vx0])
%                                                                     [y0]                     [vy0]
%                                                                     [z0]                     [vz0]

% here is the fast track to the answer, derivation below
% use a function cross_ from below to create the 4 (3x3) cross product matrix components of (C) and 2 (3x1) cross product components of (D)
C11 =  cross_(v2-v1);               % 3x3 cross product matrix associated with  (v2-v1) x r0 
C12 = -cross_(r2-r1);               % 3x3 cross product matrix associated with -(r2-r1) x v0
C21 =  cross_(v3-v1);               % 3x3 cross product matrix associated with  (v3-v1) x r0
C22 = -cross_(r3-r1);               % 3x3 cross product matrix associated with -(r3-r1) x v0
D1 = cross_(r2)*v2 - cross_(r1)*v1; % 3x1 vector associated with (r2 x v2) - (r1 x v1)
D2 = cross_(r3)*v3 - cross_(r1)*v1; % 3x1 vector associated with (r3 x v3) - (r1 x v1)

% form the (C) and (D) matrices
C = [C11 C12
     C21 C22];

D = [D1 
     D2];

% use any method to compute the inverse of (C) = (C^-1). Here is the answer
rv0 = C^-1*D

% here are the details of the derivation

% some notes on 3D-vector cross product (x) properties and their matrix form
% a x b = [ a3*b2-a2*b3]
%         [-a3*b1+a1*b3]
%         [ a2*b1-a1*b2]
% a x b = [ 0   a3 -a2] [b1] = A*b  cross product matrix, A=cross_(a), (based on vector a) is (3x3), (b) is (3x1) 
%         [-a3  0   a1] [b2]        
%         [ a2 -a1  0 ] [b3]
% a x a = [0]
%         [0] 
%         [0]
% b x a = -a x b

% trajectory of any 3 of the known initial postion & velocity vectors
% (r1-r0) = -(v1-v0)*t1
% (r2-r0) = -(v2-v0)*t2
% (r3-r0) = -(v3-v0)*t3

% take cross product of each side to eliminate time from the equations
% (v1-v0) x (r1-r0) = v1 x r1 - v1 x r0 - v0 x r1 + v0 x r0 = v1 x r1 - v1 x r0 + r1 x v0 + v0 x r0 = -(v1-v0) x (v1-v0)*t1 = [0]
%                                                                                                                             [0]
%                                                                                                                             [0]
% (v2-v0) x (r2-r0) = v2 x r2 - v2 x r0 - v0 x r2 + v0 x r0 = v2 x r2 - v2 x r0 + r2 x v0 + v0 x r0 = -(v2-v0)x(v2-v0)*t2 = [0]
%                                                                                                                           [0]
%                                                                                                                           [0]
% (v3-v0) x (r3-r0) = v3 x r3 - v3 x r0 - v0 x r3 + v0 x r0 = v3 x r3 - v3 x r0 + r3 x v0 + v0 x r0 = -(v3-v0)x(v3-v0)*t3 = [0]
%                                                                                                                           [0]
%                                                                                                                           [0]

% subtract pairs of equations to eliminate v0 x r0 leaving linear equations in r0 and v0
% (v2 x r2)-(v1 x r1) - (v2-v1) x r0 + (r2-r1) x v0 = 0
% (v3 x r3)-(v1 x r1) - (v3-v1) x r0 + (r3-r1) x v0 = 0
% rearranging
% (v2-v1) x r0 - (r2-r1) x v0 = (v2 x r2) - (v1 x r1)
% (v3-v1) x r0 - (r3-r1) x v0 = (v3 x r3) - (v1 x r1)
% or in matrix form
% C11*r0 + C12*v0 = D1
% C21*r0 + C22*v0 = D2
% [C11 C12] [r0] = [D1]
% [C21 C22] [v0]   [D2]
% C11 = cross product matrix of  (v2-v1)
% C12 = cross product matrix of -(r2-r1)
% C21 = cross product matrix of  (v3-v1)
% C22 = cross product matrix of -(r3-r1)
% D1 = [(v2 x r2)-(v1 x r1)]
% D2 = [(v3 x r3)-(v1 x r1)]
% reform into linear matrix equation
% C * rv0 = D
% rv0 = C^-1 * D
% rv0 = [x0 ]
%       [y0 ]
%       [z0 ]
%       [vx0]
%       [vy0]
%       [vz0]
% C = [C11 C12]       (6x6) matrix formed from 4 (3x3) cross product matrices
%     [C21 C22]
% D = [D1]      (6x1) matrix formed from 2 (3x1) vectors
%     [D2]
% rv0 = [C11 C12]^-1 * [D1]
%       [C21 C22]      [D2]   

% (C) and (D) matrix expansions in vector form to visualize 

x = [1 0 0];
y = [0 1 0];
z = [0 0 1];

C  = [ 0          z*(v2-v1) -y*(v2-v1)  0         -z*(r2-r1)  y*(r2-r1) 
      -z*(v2-v1)  0          x*(v2-v1)  z*(r2-r1)  0         -x*(r2-r1) 
       y*(v2-v1) -x*(v2-v1)  0         -y*(r2-r1)  x*(r2-r1)  0 
       0          z*(v3-v1) -y*(v3-v1)  0         -z*(r3-r1)  y*(r3-r1)  
      -z*(v3-v1)  0          x*(v3-v1)  z*(r3-r1)  0         -x*(r3-r1)
       y*(v3-v1) -x*(v3-v1)  0         -y*(r3-r1)  x*(r3-r1)  0];

D  = [cross(r2,v2)-cross(r1,v1)
      cross(r3,v3)-cross(r1,v1)];

rv0 = C^-1*D

% both the calculations for rv0 yield
% 
% rv0 =
% 
%            404422374079783
%            199182431001928
%            166235642339249
%                       -228
%                        166
%                        245