Distance between two circles with known size and intersection area Caution - biologist at work

I am trying to plot some circles and want to work out how far apart my circles should be with a target in mind. The ditance between the centre points of the two circles is $d$. The first circle has a radius of $r_1$ = 5 cm and the second has a radius of $r_2$ = 4 cm. The area of these circles is 78.54 cm$^2$ and 50.27 cm$^2$ respectively, making the joint area $a_{1,2}$ = 128.81 cm$^2$. I want 30% of that area, $0.3 \times 128.81$ cm$^2$ to be intersection area (area shared by the two circles). Therefore I have a known pair of radii and a known intersection area, and an unknown distance. By trial and error with various values of $d$ in the equation
$$ a_{1,2} = r_1^2 cos^{-1} \frac{d^2 + r_1^2 - r_2^2}{2dr_1} + r_2^2 cos^{-1} \frac{d^2 + r_2^2 - r_1^2}{2dr_2} - \frac{1}{2} \sqrt{(-d+r_1+r_2)(d+r_1-r_2)(d-r_1+r_2)(d+r_1+r_2)} $$
I find that, to get 30% of the total area into the intersection area, when the radii are 5 and 4, $d$ needs to be approximately 2.71. However, is there a way to solve $d$ from known radii and intersection area? I've been searching for hours and come up short - the internet is dominated by solutions to find the intersection are from known radii and distance.
FYI, I'm programming this in R - bonus kudos for answers providing accompanying script!
r1 <- 5
r2 <- 4
a1 <- pi*r1^2
a2 <- pi*r2^2
aJ <- 0.3 * (a1+a2)

intersectionArea <- function(r1, r2, d){
  (r1^2 * acos((d^2 + r1^2 - r2^2)/(2*d*r1))) +
  (r2^2 * acos((d^2 + r2^2 - r1^2)/(2*d*r2))) -
  (0.5 * sqrt((-d+r1+r2)*(d+r1-r2)*(d-r1+r2)*(d+r1+r2)))
}

aJ/intersectionArea(r1,r2, 2.709)

 A: Here is the Python script that effectively solves your problem (sorry, I don't do R). It's simple enough than it can be translated into any other programming language. The solution is obtained iteratively but it converges quite quickly (range of possible solutions is cut by half in each iteration).
import math

def area(r1, r2, d):
    s1 = r1**2 * math.acos((d**2 - r2**2 + r1**2)/(2 * d * r1))
    s2 = r2**2 * math.acos((d**2 - r1**2 + r2**2)/(2 * d * r2))
    s3 = math.sqrt((-d + r1 + r2) * (d + r1 - r2) * (d + r2 - r1) * (d + r1 + r2)) / 2
    return s1 + s2 - s3

def solve(r1, r2, targetArea, error):
    d1 = r1 + r2
    a1 = area(r1, r2, d1)
    d2 = abs(r1 - r2)
    a2 = area(r1, r2, d2)
    while True:
        d = (d1 + d2) / 2
        a = area(r1, r2, d)
        if abs(a - targetArea) <= error:
            return d
        if a > targetArea:
            d2 = d
            a2 = a
        else:
            d1 = d
            a1 = a

# input data
r1 = 4
r2 = 5
a = 0.3 * 128.81
error = 0.00001

# find d for the given input 
d = solve(r1, r2, a, error)
print("The result is d=%.5f" % d)

# calculate area for the calculated value of d
a2 = area(r1, r2, d)

# print for verification
print("For r1=%.5f, r2=%.5f, d=%.5f the area is %.5f and the target area is %.5f"
      % (r1, r2, d, a2, a))

In the hart of this script is the following function:
def solve(r1, r2, targetArea, error)        

Simply, it takes 4 parameters: two radii, intersection area and the accuracy that you want to achieve (0.00001 or 0.01 or any other small value).
The output of the script is:
The result is d=2.70874
For r1=4.00000, r2=5.00000, d=2.70874 the area is 38.64300 and the target area is 38.64300

