from math import pi,asin,sin,cos,acos
from PADS.SVG import *
import sys

# Circular segment formulas from https://en.wikipedia.org/wiki/Circular_segment

def areaFromRadius(r): return pi*r**2

def twoCircleArea(r1,r2): return areaFromRadius(r1)+areaFromRadius(r2)

def angleFromChord(r,c):
    return 2*asin(c/(2*r))

def areaFromAngle(r,theta):
    return r**2*(theta-sin(theta))/2

def areaFromChord(r,c):
    theta = angleFromChord(r,c)
    return areaFromAngle(r,theta)

def distanceToChord(r,c):
    theta = angleFromChord(r,c)
    return r*cos(theta/2)

def separationFromChord(r1,r2,c):
    return distanceToChord(r1,c)+distanceToChord(r2,c)

# Bisection to solve Mrs. Miniver's problem within floating point accuracy

def miniverSeparation(r1,r2):
    if r1 > r2: r1,r2 = r2,r1   # make sure small circle first
    totalArea = twoCircleArea(r1,r2)
    targetArea = totalArea/3

    def miniverTest(x):
        # x = cos(chord on small circle), -1: outside, 1: inside
        chord = 2*r1*sin(acos(x))
        area1 = areaFromChord(r1,chord)
        if x > 0:
            area1 = areaFromRadius(r1) - area1
        area2 = areaFromChord(r2,chord)
        return area1 + area2 < targetArea

    lo,hi = -1.0,1.0
    for i in range(100):
        mid = (hi+lo)/2
        if miniverTest(mid):
            lo = mid
        else:
            hi = mid
    chord = 2*r1*sin(acos(mid))
    return distanceToChord(r2,chord) - r1*mid

# Make it into a nice picture
# Monochromatic and with a loose frame for now; we'll fix it up later

gridUnit = 150
targetArea = 100000.0
boundingBox = (6+4j)*gridUnit
output = SVG(boundingBox,sys.stdout)
output.group(fill="none",stroke="#000")

for i in (0,1,2):
    ratio = 2.0**(i*0.25)
    areaWhenSmall = twoCircleArea(1.0,ratio)
    factor = (targetArea/areaWhenSmall)**0.5
    r1,r2 = factor,factor*ratio
    s = miniverSeparation(r1,r2)
    o = (1+2j+2*i)*gridUnit
    output.circle(o + (s+r2-r1)*0.5j, r1)
    output.circle(o - (s+r1-r2)*0.5j, r2)
  
output.ungroup()
output.close()