Thanks for the reply. In my particular instance, there are a lot of
constants, and the problem looks a bit difficult to automate. Here are
the specifics:
var('x1,x2,x3,x4,k,x,y,z', domain=RR)
# Definition of P = Im[(x1 + i x2)^k]
P(k, x1, x2, x3, x4) = (x1^2 + x2^2)^(k/2)*sin(k*arctan2(x2,x1))
#
Example; intersection of an elipsoid and sphere:
var('x,y,z')
solve([x^2 +y^2+z^2 ==1,x^2+y^2+2*z^2 ==1],[x,y,z])
#[[x == r1, y == -sqrt(-r1^2 + 1), z == 0], [x == r2, y == sqrt(-r2^2
+ 1), z == 0]]
var('r1')
p1=parametric_plot3d([r1,-sqrt(-r1^2 + 1),0],
(-1,1),thickness=10,color='red')
p2=paramet
