x,y,z = var('x, y, z')
Z(x,y) = sqrt(x^2 + y^2 + 1)
x1, y1 = [3, 4]
x2, y2 = [5, 7]
p1 = [x1,y1,Z(x1,y2)]
p2 = [x2,y2,Z(x2,y2)]
xhat(x,y) = x/(Z(x,y)+1)
yhat(x,y) = y/(Z(x,y)+1)
q1 = [xhat(x1,y1),yhat(x1,y1)]
q2 = [xhat(x2,y2),yhat(x2,y2)]
v = vector([x, y, z])
v1 = vector(p1)
v2 = vector(p2)
n = v1.cross_product(v2)
a = n[0]
b = n[1]
c = n[2]
p(x,y) = c*x^2 + c*y^2 + 2*a*x + 2*b*y + c
r(x,y) = x^2 + y^2 - 1