Area of intersection between two circles in the hyperbolic plane

code.sage

# See http://math.stackexchange.com/q/684834/35416 for details
var('a b x y r d')
var('R1', latex_name='R')
var('R2', latex_name="R'")
# The following two antiderivatives were found by Wolfram Alpha
upperArc(x) = -(b*arctan((b*x)/(sqrt(b^2-r^2)*sqrt(r^2-x^2))))/sqrt(b^2-r^2)+(b*arctan(x/sqrt(b^2-r^2)))/sqrt(b^2-r^2)+arctan(x/sqrt(r^2-x^2))
lowerArc(x) = (-sqrt(b^2-r^2)*arctan(x/sqrt(r^2-x^2))+b*arctan((b*x)/(sqrt(b^2-r^2)*sqrt(r^2-x^2)))+b*arctan(x/sqrt(b^2-r^2)))/sqrt(b^2-r^2)
y1 = cosh(R1)
r1 = sinh(R1)
y2 = (cosh(d)+sinh(d))*cosh(R2)
r2 = (cosh(d)+sinh(d))*sinh(R2)
s1 = solve([x^2+(y-y1)^2==r1^2, x^2+(y-y2)^2==r2^2], [x,y], solution_dict=True)
x1 = s1[0][x]
x2 = s1[1][x]
res = (lowerArc(x2)-lowerArc(x1)).subs(b=y2,r=r2) - (upperArc(x2)-upperArc(x1)).subs(b=y1,r=r1)

result.txt

2*((cosh(d) + sinh(d))*arctan(1/2*sqrt(-8*cosh(d)*sinh(d)^3 - 8*sinh(d)^4 - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 8*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + cosh(d) + 2*sinh(d))*cosh(R1)*cosh(R2) - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - 4*sinh(d)^2)*(cosh(d) + sinh(d))*cosh(R2)/((sqrt(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*(cosh(d) + sinh(d))*cosh(R2) - sqrt(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1))*sqrt((2*(4*cosh(d)*sinh(d)^2 + 4*sinh(d)^3 + cosh(d) + 3*sinh(d))*cosh(R1)*cosh(R2)^3 - (8*cosh(d)*sinh(d)^3 + 8*sinh(d)^4 + 4*cosh(d)*sinh(d) + 8*sinh(d)^2 + 1)*cosh(R2)^4 - 2*cosh(d)*sinh(d)^3 - 2*sinh(d)^4 - 2*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + sinh(d))*cosh(R1)*cosh(R2) + (8*cosh(d)*sinh(d)^3 + 8*sinh(d)^4 - (2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 2*cosh(d)*sinh(d) + 6*sinh(d)^2)*cosh(R2)^2 - sinh(d)^2)/(2*(cosh(d) + sinh(d))*cosh(R1)*cosh(R2) - (2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - cosh(R1)^2))))*cosh(R2) + (cosh(d) + sinh(d))*arctan(1/2*sqrt(-8*cosh(d)*sinh(d)^3 - 8*sinh(d)^4 - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 8*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + cosh(d) + 2*sinh(d))*cosh(R1)*cosh(R2) - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - 4*sinh(d)^2)/(((cosh(d) + sinh(d))*cosh(R2) - cosh(R1))*sqrt(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)))*cosh(R2) + (arctan(1/2*sqrt(-8*cosh(d)*sinh(d)^3 - 8*sinh(d)^4 - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 8*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + cosh(d) + 2*sinh(d))*cosh(R1)*cosh(R2) - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - 4*sinh(d)^2)*cosh(R1)/(((cosh(d) + sinh(d))*cosh(R2) - cosh(R1))*sqrt(-((2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2*cosh(R2)^2 + cosh(R1)^4 + 2*cosh(d)*sinh(d)^3 + 2*sinh(d)^4 + 2*(cosh(d)*sinh(d) + sinh(d)^2)*cosh(R1)^2 - 2*((cosh(d) + sinh(d))*cosh(R1)^3 + (2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + sinh(d))*cosh(R1))*cosh(R2) + sinh(d)^2)/(2*(cosh(d) + sinh(d))*cosh(R1)*cosh(R2) - (2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - cosh(R1)^2))))*cosh(R1) - arctan(1/2*sqrt(-8*cosh(d)*sinh(d)^3 - 8*sinh(d)^4 - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 8*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + cosh(d) + 2*sinh(d))*cosh(R1)*cosh(R2) - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - 4*sinh(d)^2)/((cosh(d) + sinh(d))*cosh(R2) - cosh(R1)))*cosh(R1) - arctan(1/2*sqrt(-8*cosh(d)*sinh(d)^3 - 8*sinh(d)^4 - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 8*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + cosh(d) + 2*sinh(d))*cosh(R1)*cosh(R2) - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - 4*sinh(d)^2)/(((cosh(d) + sinh(d))*cosh(R2) - cosh(R1))*sqrt(-((2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2*cosh(R2)^2 + cosh(R1)^4 + 2*cosh(d)*sinh(d)^3 + 2*sinh(d)^4 + 2*(cosh(d)*sinh(d) + sinh(d)^2)*cosh(R1)^2 - 2*((cosh(d) + sinh(d))*cosh(R1)^3 + (2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + sinh(d))*cosh(R1))*cosh(R2) + sinh(d)^2)/(2*(cosh(d) + sinh(d))*cosh(R1)*cosh(R2) - (2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - cosh(R1)^2)))) - arctan(1/2*sqrt(-8*cosh(d)*sinh(d)^3 - 8*sinh(d)^4 - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 8*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + cosh(d) + 2*sinh(d))*cosh(R1)*cosh(R2) - 4*(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - 4*sinh(d)^2)/(((cosh(d) + sinh(d))*cosh(R2) - cosh(R1))*sqrt((2*(4*cosh(d)*sinh(d)^2 + 4*sinh(d)^3 + cosh(d) + 3*sinh(d))*cosh(R1)*cosh(R2)^3 - (8*cosh(d)*sinh(d)^3 + 8*sinh(d)^4 + 4*cosh(d)*sinh(d) + 8*sinh(d)^2 + 1)*cosh(R2)^4 - 2*cosh(d)*sinh(d)^3 - 2*sinh(d)^4 - 2*(2*cosh(d)*sinh(d)^2 + 2*sinh(d)^3 + sinh(d))*cosh(R1)*cosh(R2) + (8*cosh(d)*sinh(d)^3 + 8*sinh(d)^4 - (2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R1)^2 + 2*cosh(d)*sinh(d) + 6*sinh(d)^2)*cosh(R2)^2 - sinh(d)^2)/(2*(cosh(d) + sinh(d))*cosh(R1)*cosh(R2) - (2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)*cosh(R2)^2 - cosh(R1)^2)))))*sqrt(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1))/sqrt(2*cosh(d)*sinh(d) + 2*sinh(d)^2 + 1)