% SphereProblem.m
% Solution to intersecting sphere problem
% Data provided - each sphere is centred at (xc,yc,zc) and has radius r
xc1=0; yc1=0; zc1=1200; r1=646.12;
xc2=-957.0680502; yc2=1350.871118; zc2=100; r2=1548.28;
xc3=-1706.834233; yc3=712.6642874; zc3=800; r3= 1249.19;
xc4=-1464.077028; yc4=-532.89593; zc4=600; r4=1267.62;

% These each give rise to a sphere equation
% (x-xc)^2+(y-yc)^2+(z-zc)^2 = r^2

% We could solve the problem with only 3 data points, but using 4 makes it
% easier to solve -> the fourth point turns solving a set of quadratic
% problems into solving a system of linear equations

% Consider taking the first two sphere equations
% (x-xc1)^2+(y-yc1)^2+(z-zc1)^2 = r1^2   [eqn. 1]
% (x-xc2)^2+(y-yc2)^2+(z-zc2)^2 = r2^2   [eqn. 2]

% Now expand them:
% x^2-2x*xc1 + xc1^2 + y^2 - 2y*yc1 + yc1^2 + z^2 - 2z*zc1 + zc1^2 = r1^2 [eqn. 3]
% x^2-2x*xc2 + xc2^2 + y^2 - 2y*yc2 + yc2^2 + z^2 - 2z*zc2 + zc2^2 = r2^2 [eqn. 4]

% And take the difference of eqns [3] - [4]:
% x*(2xc2-2xc1) + y*(2yc2-2yc1) + z*(2zc2-2zc1) = r1^2 - xc1^2 - yc1^2 - zc1^2
%                                                -r2^2 + xc2^2 + yc2^2 + zc2^2

% Repeat by subtracting the sphere equations for spheres 3 and 4, yielding
% three linear equations in (x,y,z) - this can then be solved as a system of
% linear equations A*[x;y;z] = b, where the matrix A consists of the
% coefficients for (x,y,z) and b the constant terms.

% Solve using MATLAB
A=[2*(xc2-xc1) 2*(yc2-yc1) 2*(zc2-zc1);
   2*(xc3-xc1) 2*(yc3-yc1) 2*(zc3-zc1);
   2*(xc4-xc1) 2*(yc4-yc1) 2*(zc4-zc1)];

b=[r1^2-xc1^2-xc1^2-zc1^2-r2^2+xc2^2+yc2^2+zc2^2;
   r1^2-xc1^2-yc1^2-zc1^2-r3^2+xc3^2+yc3^2+zc3^2;
   r1^2-xc1^2-yc1^2-zc1^2-r4^2+xc4^2+yc4^2+zc4^2];

% To solve Ax=b, find solution = A^{-1}*b, since (we assume) the rows in A
% are linearly independent.

solution=inv(A)*b;
x_sol=solution(1)
y_sol=solution(2)
z_sol=solution(3)

% Output:
% x_sol = -595.0350
% y_sol =  235.5398
% z_sol = 1.1110e+03 -> I guess this is what you're looking for?  Roughly an elevation of 1111

% Now check that each of the four equations is satisfied - each of these should
% give zero as the answer
(x_sol-xc1)^2+(y_sol-yc1)^2+(z_sol-zc1)^2 - r1^2
(x_sol-xc2)^2+(y_sol-yc2)^2+(z_sol-zc2)^2 - r2^2
(x_sol-xc3)^2+(y_sol-yc3)^2+(z_sol-zc3)^2 - r3^2
(x_sol-xc4)^2+(y_sol-yc4)^2+(z_sol-zc4)^2 - r4^2

% Output for all four of the equations is -5.4964 - which is pretty close, since it's comparing the
% point on the sphere to the squared radius (which is on the order of 10^6 units^2).