use Matrix.frink

symbolicMode[true]
showApproximations[false]
mx = new Matrix[noEval[[[cos[phi], sin[phi]], [-sin[phi], cos[phi]]]]]
m1 = (new Matrix[[(x1-x2)/2, (y1-y2)/2]]).transpose[]

mt = mx.multiply[m1]
println[mt.formatMatrix[]]

x1p = mt.get[1,1]
y1p = mt.get[2,1]
println[x1p]
println[y1p]

part1 = sqrt[(rx^2 ry^2 - rx^2 y1p^2 - ry^2 x1p^2) / (rx^2 y1p^2 - ry^2 x1p^2)]
if fa == fs
   part1 = -part1

cxp =   part1 rx y1p / ry
cyp = - part1 ry x1p / rx

mc = new Matrix[noEval[[[cos[phi], -sin[phi]], [sin[phi], cos[phi]]]]]
m2 = (new Matrix[noEval[[cxp, cyp]]]).transpose[]
m3 = (new Matrix[noEval[[(x1+x2)/2, (y1+y2)/2]]]).transpose[]
println[m3.formatMatrix[]]

mr = mc.multiply[m2].add[m3]
println[mr.formatMatrix[]]
cx = mr.get[1,1]
cy = mr.get[2,1]
println[cx]
println[cy]

dotProduct[u, v] := sum[mul[u, v]]

dotProductNormalized[u, v] :=
{
   clamp[dotProduct[u, v] / (hypotenuse[u] * hypotenuse[v]), -1, 1]
}

// Eq 5.4
angleBetween[u, v] :=
{
   sign = realSignum[u@0 * v@1 - u@1 * v@0]
   println["sign is $sign"]
   return sign arccos[dotProductNormalized[u,v]]
}

vx = [1,0]
v1 = [(x1p - cxp) / rx, (y1p - cyp)/ry]

theta1 = angleBetween[vx, v1]
println["theta1 is $theta1"]

v2 = [(-x1p - cxp) / rx, (-y1p - cyp)/ry]
println["v1 is $v1"]
println["v2 is $v2"]
dtheta = angleBetween[v1, v2] mod circle
println["dtheta is $dtheta"]