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()