kill(all)$
load("lrats");
isolate_wrt_times : true$

/* the three equations */
eq1: (x1-x)^2+(y1-y)^2+(z1-z)^2=l1^2$
eq2: (x2-x)^2+(y2-y)^2+(z2-z)^2=l2^2$
eq3: (x3-x)^2+(y3-y)^2+(z3-z)^2=l3^2$

/* expand */
eq1_expanded : fullratsimp(expand(eq1))$
eq2_expanded : fullratsimp(expand(eq2))$
eq3_expanded : fullratsimp(expand(eq3))$

/* subtract eq2 and eq3 from eq1 */
eq1_m_eq2 : rat(eq1_expanded - eq2_expanded)$
eq1_m_eq3 : rat(eq1_expanded - eq3_expanded)$

result : solve( [eq1_m_eq2, eq1_m_eq3],[x,y])$

/* substitute some 'patterned and beautiful' shared subexpressions to simplify the expressions */
result : lratsubst( [x1^2+y1^2+z1^2=w1, x2^2 + y2^2 + z2^2 = w2, x3^2 + y3^2 + z3^2 = w3 ], result )$
result : lratsubst( [(y1-y2)*z3 + (y3-y1)*z2 + (y2-y3)*z1=v1, (x2-x1)*z3 + (x1-x3)*z2 + (x3-x2)*z1=v2, (x2-x1)*y3+(x1-x3)*y2+(x3-x2)*y1=d], result)$

disp("Result after substitutions:")$
disp(result)$

x_result : expand(isolate(rhs(result[1][1]),z));
y_result : expand(isolate(rhs(result[1][2]),z));

/* insert resulting x and y into eq1 */
eq1_substituted : subst(x_result, x, eq1_expanded)$
eq1_substituted : subst(y_result, y, eq1_substituted)$

z_final : solve(eq1_substituted,z)[1]$

/* boooom */
z_final : rat( z_final )$
disp( "z independent:")$
disp( z_final )$

/* trying to improve on the form, but doesn't work :( */
/* this one should simplify the first subexpression inside the sqrt as well as the denominator of the whole expression and a whole lot of other places, but it doesn't have any effect :( */
z_final : lratsubst( [%t19^2+%t16^2=c1], z_final);

/* also I can't find out how to factor out the obvious common factors and constansts like '2' */
/* this doesn't work: */
z_final : factorsum( z_final );

disp("finished")$