In [1]:

    

%pylab inline

from scipy.optimize import basinhopping
from scipy.interpolate import UnivariateSpline

import pandas as pd

import csv

import shapely.affinity as af
import shapely.geometry as sh

from equivalent_ellipse import *
from scaled_figures import *


    

    




Populating the interactive namespace from numpy and matplotlib

In [46]:

    

x_list = list()
y_list = list()

with open('../data/cutout_x.csv', 'r') as x_csvfile:
    with open('../data/cutout_y.csv', 'r') as y_csvfile:
        
        x_reader = csv.reader(x_csvfile, delimiter=',', lineterminator='\n')
        y_reader = csv.reader(y_csvfile, delimiter=',', lineterminator='\n')
        
        for row in x_reader:
            x_list += [row]
  
        for row in y_reader:
            y_list += [row]

            
num_cutouts = len(x_list)


x_array = [0,]*num_cutouts
y_array = [0,]*num_cutouts


for i in range(num_cutouts):

    x_array[i] = array(x_list[i], dtype='float')
    y_array[i] = array(y_list[i], dtype='float')

    
cutout = [0,]*num_cutouts

for i in range(num_cutouts):
    cutout[i] = shapely_cutout(x_array[i],y_array[i])


    

In [47]:

    

with open('cutout_factors','r') as f:
    
    loaded_factor = eval(f.read())
    
dictArray = array(list(loaded_factor.items()),float)
cutoutRef = dictArray[:,0].astype(int)
factor = dictArray[:,1]


    

In [48]:

    

len(cutout)


    

    
Out[48]:





116

In [49]:

    

with open('custom_cutout_factors','r') as f:
    
    loaded_custom_factor = eval(f.read())

# customDimensionsDict = {'5x13':[5,13],
#                    '5x10':[5,10],
#                    '5x8':[5,8],
#                    '4x13':[4,13],
#                    '4x10':[4,10],
#                    '4x8':[4,8],
#                    '4x6.5':[4,6.5],
#                    '3x13':[3,13],
#                    '3x9':[3,9],
#                    '3x6.5':[3,6.5],
#                    '3x5':[3,5]}

customDimensionsDict = {'5x8':[5,8],
                   '4x8':[4,8],
                   '4x6.5':[4,6.5],
                   '3x9':[3,9],
                   '3x6.5':[3,6.5],
                   '3x5':[3,5]}

customDataLabelList = list(customDimensionsDict.keys())

custom_factors = [loaded_custom_factor[x] for x in customDataLabelList]
custom_widths = [customDimensionsDict[x][0] for x in customDataLabelList]
custom_lengths = [customDimensionsDict[x][1] for x in customDataLabelList]

ellipse_cutouts = [0,]*len(custom_widths)

for i in range(len(custom_widths)):
    
    ellipse_cutouts[i] = create_ellipse([0,0,custom_widths[i],custom_lengths[i],-45])
    

cutoutRef = append(cutoutRef,arange(len(cutout),len(cutout)+len(ellipse_cutouts)))
cutout.extend(ellipse_cutouts)
factor = append(factor,custom_factors)


    

In [50]:

    

# testCutout = cutout[70]
# testCutout


    

In [7]:

    

with open('circle_cutout_factors','r') as f:
    
    loaded_circle_factor = eval(f.read())
    
circleDictArray = array(list(loaded_circle_factor.items()),float)

circleDiameter = circleDictArray[:,0]
circleFactor = circleDictArray[:,1]


    

In [8]:

    

scatter(circleDiameter,1/circleFactor)

circleFit = UnivariateSpline(circleDiameter,1/circleFactor)

diameter_i = linspace(circleDiameter.min(),circleDiameter.max())
plot(diameter_i,circleFit(diameter_i),'r--')


    

    
Out[8]:





[<matplotlib.lines.Line2D at 0x90635c0>]






    

In [9]:

    

minVal = 3
maxVal = 9
numZones = 100


    

In [10]:

    

# Create Zones

zoneWidth = linspace(minVal,maxVal,numZones)

zoneWeight = 70*circleFit.derivative()(zoneWidth)

circle = [0,]*(numZones)
circle[0] = sh.Point(0,0).buffer(minVal/2)

zone = [0,]*numZones
zone[0] = circle[0]

for i in range(1,numZones-1):

    circle[i] = sh.Point(0,0).buffer(zoneWidth[i]/2)
    zone[i] = circle[i].difference(circle[i-1])

zone[-1] = sh.Point(0,0).buffer(20).difference(circle[-2])

weightedArea = 0.
for i in range(numZones):
    
    weightedArea += cutout[0].intersection(zone[i]).area * zoneWeight[i]
    
weightedArea


    

    
Out[10]:





46.563494077548761

In [11]:

    

cutout[0].area


    

    
Out[11]:





40.994999999999976

In [12]:

    

plot(zoneWidth,zoneWeight)


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0x95800b8>]






    

In [13]:

    

def weightedAreaOptimise(x,inputShape):
    
    a = x[0]
    b = x[1]
    
    movedShape =af.translate(inputShape, xoff=-a, yoff=-b)
    
    weightedArea = 0.
    
    for i in range(numZones):
        weightedArea += movedShape.intersection(zone[i]).area * zoneWeight[i]
        
    return -weightedArea


    

In [14]:

    

class MyTakeStep(object):
    
    def __init__(self, stepsize=0.5):
        self.stepsize = stepsize
        
    def __call__(self, x):
        s = self.stepsize
        x[0] += np.random.uniform(-3*s, 3*s)
        x[1] += np.random.uniform(-3*s, 3*s)
        
        return x
    

def print_fun(x, f, accepted):

    global successCount_basinhopping
    global minVals_basinhopping
    global nSuccess_basinhopping

    a = x[0]
    b = x[1]
    
    f = np.around(f,4)
    
    if accepted:
        
        if np.all(np.isnan(minVals_basinhopping)):
            
            minVals_basinhopping[0] = f
            successCount_basinhopping = 1
            
            print(("first local minima of %.4f at:\n"+
                "[a,b] = [%.4f, %.4f]")
                % (f,a,b))
            
        elif (np.nanmin(minVals_basinhopping) == f):
        
            minVals_basinhopping[successCount_basinhopping] = f
            successCount_basinhopping += 1
            
            print(("agreeing local minima of %.4f at:\n"+
                "[a,b] = [%.4f, %.4f]")
                % (f,a,b))
        
        elif (np.nanmin(minVals_basinhopping) > f):
        
            minVals_basinhopping[0] = f
            successCount_basinhopping = 1
            
            print(("new local minima of %.4f at:\n"+
                "[a,b] = [%.4f, %.4f]") 
                % (f,a,b))
                
        else:
        
            print(("rejected local minima of %.4f at:\n"+
                "[a,b] = [%.4f, %.4f]") 
                % (f,a,b))
            
        
    else:
        
        print(("failed to find local minima at:\n"+
            "[a,b] = [%.4f, %.4f]") 
            % (a,b))
          

    if successCount_basinhopping >= nSuccess_basinhopping:
        return True



def optimise_middle(cutout,userNSuccess):
    
    mytakestep = MyTakeStep()
    
    global nSuccess_basinhopping
    nSuccess_basinhopping = userNSuccess
    
    global minVals_basinhopping
    minVals_basinhopping = np.empty(nSuccess_basinhopping)
    minVals_basinhopping[:] = np.nan
    
    global successCount_basinhopping
    successCount_basinhopping = 0
        
    minimizer_kwargs = {"args": (cutout,), "method": 'L-BFGS-B', 
                        "bounds": ((-5,5),(-5,5)) }
    
    x0 = np.array([0,0])

    output = basinhopping(weightedAreaOptimise,x0,
                          niter=1000,minimizer_kwargs=minimizer_kwargs,
                          take_step=mytakestep, callback=print_fun)

    return output


    

In [15]:

    

middle_points = zeros([len(cutoutRef),2])

for j in range(len(cutoutRef)):
    
    i = cutoutRef[j]
    
    print("Cutout %.f" % (i))
    output = optimise_middle(cutout[i],3)
    middle_points[j,:] = output.x
    print(" ")


    

    




Cutout 18
first local minima of -52.7387 at:
[a,b] = [-0.1359, -0.4053]
agreeing local minima of -52.7387 at:
[a,b] = [-0.1359, -0.4053]
agreeing local minima of -52.7387 at:
[a,b] = [-0.1359, -0.4053]
 
Cutout 58
first local minima of -43.4399 at:
[a,b] = [-0.0655, 0.1871]
agreeing local minima of -43.4399 at:
[a,b] = [-0.0655, 0.1871]
agreeing local minima of -43.4399 at:
[a,b] = [-0.0655, 0.1871]
 
Cutout 20
first local minima of -26.0170 at:
[a,b] = [-0.1984, -0.3417]
agreeing local minima of -26.0170 at:
[a,b] = [-0.1984, -0.3417]
agreeing local minima of -26.0170 at:
[a,b] = [-0.1984, -0.3417]
 
Cutout 33
first local minima of -39.4141 at:
[a,b] = [-0.4710, 0.3438]
agreeing local minima of -39.4141 at:
[a,b] = [-0.4710, 0.3438]
agreeing local minima of -39.4141 at:
[a,b] = [-0.4710, 0.3438]
 
Cutout 106
first local minima of -47.8430 at:
[a,b] = [-0.7973, -0.1224]
agreeing local minima of -47.8430 at:
[a,b] = [-0.7973, -0.1224]
agreeing local minima of -47.8430 at:
[a,b] = [-0.7973, -0.1224]
 
Cutout 82
first local minima of -47.0314 at:
[a,b] = [-0.4594, -0.4578]
agreeing local minima of -47.0314 at:
[a,b] = [-0.4594, -0.4578]
agreeing local minima of -47.0314 at:
[a,b] = [-0.4594, -0.4578]
 
Cutout 6
first local minima of -51.2994 at:
[a,b] = [0.2879, -0.3353]
agreeing local minima of -51.2994 at:
[a,b] = [0.2879, -0.3353]
agreeing local minima of -51.2994 at:
[a,b] = [0.2879, -0.3353]
 
Cutout 112
first local minima of -18.3410 at:
[a,b] = [0.2009, -0.0613]
agreeing local minima of -18.3410 at:
[a,b] = [0.2009, -0.0613]
agreeing local minima of -18.3410 at:
[a,b] = [0.2009, -0.0613]
 
Cutout 104
first local minima of -44.2003 at:
[a,b] = [-0.7110, 0.5318]
agreeing local minima of -44.2003 at:
[a,b] = [-0.7110, 0.5318]
agreeing local minima of -44.2003 at:
[a,b] = [-0.7110, 0.5318]
 
Cutout 34
first local minima of -51.1747 at:
[a,b] = [0.3099, -0.3208]
agreeing local minima of -51.1747 at:
[a,b] = [0.3099, -0.3208]
agreeing local minima of -51.1747 at:
[a,b] = [0.3099, -0.3208]
 
Cutout 19
first local minima of -46.2550 at:
[a,b] = [0.4125, 0.2889]
agreeing local minima of -46.2550 at:
[a,b] = [0.4125, 0.2889]
agreeing local minima of -46.2550 at:
[a,b] = [0.4125, 0.2889]
 
Cutout 73
first local minima of -52.1552 at:
[a,b] = [-0.1866, -0.3155]
agreeing local minima of -52.1552 at:
[a,b] = [-0.1866, -0.3155]
agreeing local minima of -52.1552 at:
[a,b] = [-0.1866, -0.3155]
 
Cutout 70
first local minima of -52.1909 at:
[a,b] = [-0.0335, -0.0770]
agreeing local minima of -52.1909 at:
[a,b] = [-0.0335, -0.0770]
agreeing local minima of -52.1909 at:
[a,b] = [-0.0335, -0.0770]
 
Cutout 41
first local minima of -51.0287 at:
[a,b] = [0.0219, -0.0656]
agreeing local minima of -51.0287 at:
[a,b] = [0.0219, -0.0656]
agreeing local minima of -51.0287 at:
[a,b] = [0.0219, -0.0656]
 
Cutout 16
first local minima of -48.4960 at:
[a,b] = [0.1155, -0.5297]
agreeing local minima of -48.4960 at:
[a,b] = [0.1155, -0.5297]
agreeing local minima of -48.4960 at:
[a,b] = [0.1155, -0.5296]
 
Cutout 38
first local minima of -30.7432 at:
[a,b] = [-0.0221, -0.1692]
agreeing local minima of -30.7432 at:
[a,b] = [-0.0221, -0.1692]
agreeing local minima of -30.7432 at:
[a,b] = [-0.0221, -0.1692]
 
Cutout 109
first local minima of -46.9829 at:
[a,b] = [0.4073, -0.0851]
agreeing local minima of -46.9829 at:
[a,b] = [0.4073, -0.0851]
agreeing local minima of -46.9829 at:
[a,b] = [0.4073, -0.0851]
 
Cutout 3
first local minima of -31.2122 at:
[a,b] = [0.2018, -0.1626]
agreeing local minima of -31.2122 at:
[a,b] = [0.2018, -0.1626]
agreeing local minima of -31.2122 at:
[a,b] = [0.2018, -0.1626]
 
Cutout 43
first local minima of -45.6972 at:
[a,b] = [-0.1034, 0.1104]
agreeing local minima of -45.6972 at:
[a,b] = [-0.1034, 0.1104]
agreeing local minima of -45.6972 at:
[a,b] = [-0.1034, 0.1104]
 
Cutout 14
first local minima of -27.8475 at:
[a,b] = [0.2145, -0.4803]
agreeing local minima of -27.8475 at:
[a,b] = [0.2145, -0.4803]
agreeing local minima of -27.8475 at:
[a,b] = [0.2145, -0.4803]
 
Cutout 22
first local minima of -43.5119 at:
[a,b] = [-0.7864, -0.2803]
agreeing local minima of -43.5119 at:
[a,b] = [-0.7864, -0.2803]
agreeing local minima of -43.5119 at:
[a,b] = [-0.7864, -0.2803]
 
Cutout 83
first local minima of -44.4189 at:
[a,b] = [-0.9573, -0.3625]
agreeing local minima of -44.4189 at:
[a,b] = [-0.9573, -0.3625]
agreeing local minima of -44.4189 at:
[a,b] = [-0.9573, -0.3625]
 
Cutout 32
first local minima of -52.7161 at:
[a,b] = [0.5823, 0.3608]
agreeing local minima of -52.7161 at:
[a,b] = [0.5823, 0.3608]
agreeing local minima of -52.7161 at:
[a,b] = [0.5823, 0.3608]
 
Cutout 57
first local minima of -20.5940 at:
[a,b] = [-0.2101, -0.1810]
agreeing local minima of -20.5940 at:
[a,b] = [-0.2101, -0.1810]
agreeing local minima of -20.5940 at:
[a,b] = [-0.2101, -0.1810]
 
Cutout 116
first local minima of -29.4762 at:
[a,b] = [-0.0000, -0.0000]
agreeing local minima of -29.4762 at:
[a,b] = [0.0000, 0.0000]
agreeing local minima of -29.4762 at:
[a,b] = [-0.0000, -0.0000]
 
Cutout 117
first local minima of -19.9398 at:
[a,b] = [0.0002, 0.0001]
agreeing local minima of -19.9398 at:
[a,b] = [0.0049, 0.0054]
agreeing local minima of -19.9398 at:
[a,b] = [0.0016, -0.0022]
 
Cutout 118
first local minima of -19.9398 at:
[a,b] = [0.0007, 0.0010]
agreeing local minima of -19.9398 at:
[a,b] = [-0.0002, 0.0000]
agreeing local minima of -19.9398 at:
[a,b] = [-0.0001, 0.0002]
 
Cutout 119
first local minima of -11.4815 at:
[a,b] = [0.0000, 0.0000]
agreeing local minima of -11.4815 at:
[a,b] = [-0.0000, 0.0000]
agreeing local minima of -11.4815 at:
[a,b] = [0.0000, 0.0000]
 
Cutout 120
first local minima of -11.4815 at:
[a,b] = [-0.0000, 0.0000]
agreeing local minima of -11.4815 at:
[a,b] = [-0.0000, 0.0000]
agreeing local minima of -11.4815 at:
[a,b] = [0.0000, 0.0000]
 
Cutout 121
first local minima of -11.4815 at:
[a,b] = [0.0000, 0.0000]
agreeing local minima of -11.4815 at:
[a,b] = [0.0000, -0.0000]
agreeing local minima of -11.4815 at:
[a,b] = [0.0000, 0.0000]
 

In [15]:

    




    

In [15]:

    




    

In [51]:

    

predicted_factor = zeros(len(cutoutRef))

for j in range(len(cutoutRef)):

    i = cutoutRef[j]
    
    cutoutLineString = sh.LineString(cutout[i].exterior.coords)

    steps = linspace(0,1,1000)

    xcoords = zeros(len(steps))
    ycoords = zeros(len(steps))

    for k in range(len(steps)):


        xycoords = cutoutLineString.interpolate(steps[k], normalized=True).xy

        xcoords[k] = xycoords[0][0]
        ycoords[k] = xycoords[1][0]
    
    
    x_vals = xcoords
    y_vals = ycoords

    a = middle_points[j,0]
    b = middle_points[j,1]

    x_vals = array(x_vals,'float') - a
    y_vals = array(y_vals,'float') - b

    angles = arctan2(y_vals, x_vals) * 180 / pi

    clockwise = sum(mod(diff(angles),360)) > sum(mod(-diff(angles),360))
    diffAdj = 1 - 2*clockwise

    sector_angles = mod(diffAdj*diff(angles),360)
    wrong_direction = sector_angles > 180

    while any(wrong_direction):

        wrong_direction = append(wrong_direction,False)

        x_vals = x_vals[~wrong_direction]
        y_vals = y_vals[~wrong_direction]

        angles = arctan2(y_vals, x_vals) * 180 / pi
        sector_angles = mod(diffAdj*diff(angles),360)

        wrong_direction = sector_angles > 180

    radii = hypot(x_vals,y_vals)
    sector_radii = radii[0:-1] + diff(radii)/2

    angles = arctan2(y_vals, x_vals) * 180 / pi
    sector_angles = mod(diffAdj*diff(angles),360)

    predicted_factor[j] = 1/sum(circleFit(sector_radii*2) * sector_angles / 360)


    

In [52]:

    

predicted_factor


    

    
Out[52]:





array([ 0.99380498,  1.0056196 ,  1.0357941 ,  1.0113048 ,  1.00018241,
        1.00054378,  0.99658809,  1.05380847,  1.00491797,  0.99526787,
        1.00251433,  0.99412503,  0.99379444,  0.99568008,  0.99974382,
        1.02646332,  1.00054398,  1.0258009 ,  1.00233149,  1.03208727,
        1.00707405,  1.00436753,  0.99821249,  1.04833487,  1.01098432,
        1.02946681,  1.02358438,  1.04606499,  1.03914319,  1.05465967])

In [53]:

    

difference = factor-predicted_factor

df = pd.DataFrame(transpose([cutoutRef,predicted_factor, factor,difference]))
df.columns = ['Reference','sector_predict', 'measured_factor','difference']

df.to_csv(path_or_buf = "sector_integration.csv", index=False)

df


    

    
Out[53]:






  

    

      

      
Reference
      
sector_predict
      
measured_factor
      
difference
    
  
  

    

      
0 
      
  18
      
 0.993805
      
 0.998054
      
 0.004249
    
    

      
1 
      
  58
      
 1.005620
      
 1.008999
      
 0.003380
    
    

      
2 
      
  20
      
 1.035794
      
 1.037967
      
 0.002173
    
    

      
3 
      
  33
      
 1.011305
      
 1.006789
      
-0.004515
    
    

      
4 
      
 106
      
 1.000182
      
 1.000871
      
 0.000688
    
    

      
5 
      
  82
      
 1.000544
      
 1.000867
      
 0.000323
    
    

      
6 
      
   6
      
 0.996588
      
 1.001302
      
 0.004714
    
    

      
7 
      
 112
      
 1.053808
      
 1.053904
      
 0.000095
    
    

      
8 
      
 104
      
 1.004918
      
 1.006601
      
 0.001683
    
    

      
9 
      
  34
      
 0.995268
      
 0.993330
      
-0.001937
    
    

      
10
      
  19
      
 1.002514
      
 1.001089
      
-0.001425
    
    

      
11
      
  73
      
 0.994125
      
 0.995472
      
 0.001347
    
    

      
12
      
  70
      
 0.993794
      
 0.997190
      
 0.003396
    
    

      
13
      
  41
      
 0.995680
      
 0.994829
      
-0.000851
    
    

      
14
      
  16
      
 0.999744
      
 0.999567
      
-0.000177
    
    

      
15
      
  38
      
 1.026463
      
 1.027190
      
 0.000727
    
    

      
16
      
 109
      
 1.000544
      
 1.006979
      
 0.006435
    
    

      
17
      
   3
      
 1.025801
      
 1.032704
      
 0.006904
    
    

      
18
      
  43
      
 1.002331
      
 1.007783
      
 0.005451
    
    

      
19
      
  14
      
 1.032087
      
 1.043027
      
 0.010939
    
    

      
20
      
  22
      
 1.007074
      
 1.011713
      
 0.004639
    
    

      
21
      
  83
      
 1.004368
      
 1.010552
      
 0.006184
    
    

      
22
      
  32
      
 0.998212
      
 0.993330
      
-0.004882
    
    

      
23
      
  57
      
 1.048335
      
 1.059704
      
 0.011370
    
    

      
24
      
 116
      
 1.010984
      
 1.017153
      
 0.006169
    
    

      
25
      
 117
      
 1.029467
      
 1.040668
      
 0.011201
    
    

      
26
      
 118
      
 1.023584
      
 1.030055
      
 0.006470
    
    

      
27
      
 119
      
 1.046065
      
 1.054423
      
 0.008358
    
    

      
28
      
 120
      
 1.039143
      
 1.053516
      
 0.014373
    
    

      
29
      
 121
      
 1.054660
      
 1.069714
      
 0.015054
    
  

In [54]:

    

std(difference)


    

    
Out[54]:





0.0050329850430285471

In [70]:

    

# pos = 120
# ref = cutoutRef[pos]
# ref

ref = 57
pos = find(ref == cutoutRef)[0]


    

In [71]:

    

# i = 19
# i = 106
i = ref

# cutout[i] = cutout[i].buffer(0.001,resolution=100)


scaled_fig_start(12,12)

x_vals_show = cutout[i].exterior.xy[0]
y_vals_show = cutout[i].exterior.xy[1]

x_vals_show = array(x_vals_show,'float')
y_vals_show = array(y_vals_show,'float')

a = middle_points[pos,0]
b = middle_points[pos,1]

scatter(a,b,c='red',s=100, zorder=10)

for k in range(len(x_vals_show)):
    
    plot([a,x_vals_show[k]],[b,y_vals_show[k]],'g')
    
plot(x_vals_show,y_vals_show,linewidth=3)

scaled_fig_end(12,12)


    

In [72]:

    

test = sh.LineString(cutout[i].exterior.coords)

steps = linspace(0,1,100)


xcoords = zeros(len(steps))
ycoords = zeros(len(steps))

for i in range(len(steps)):
    
    
    xycoords = test.interpolate(steps[i], normalized=True).xy
    
    xcoords[i] = xycoords[0][0]
    ycoords[i] = xycoords[1][0]


    

In [73]:

    

scaled_fig_start(12,12)

a = middle_points[pos,0]
b = middle_points[pos,1]

scatter(a,b,c='red',s=100, zorder=10)

for k in range(len(xcoords)):
    
    plot([a,xcoords[k]],[b,ycoords[k]],'g')
    
plot(xcoords,ycoords,linewidth=3)

scaled_fig_end(12,12)


    

In [74]:

    

i = 19

x_vals = xcoords
y_vals = ycoords

a = middle_points[i,0]
b = middle_points[i,1]

x_vals = array(x_vals,'float') - a
y_vals = array(y_vals,'float') - b

angles = arctan2(y_vals, x_vals) * 180 / pi
plot(angles)


    

    
Out[74]:





[<matplotlib.lines.Line2D at 0xb3f6860>]






    

In [75]:

    

clockwise = sum(mod(diff(angles),360)) > sum(mod(-diff(angles),360))
diffAdj = 1 - 2*clockwise

sector_angles = mod(diffAdj*diff(angles),360)
plot(sector_angles)


    

    
Out[75]:





[<matplotlib.lines.Line2D at 0xb44aa90>]






    

In [76]:

    

wrong_direction = sector_angles > 180

while any(wrong_direction):

    wrong_direction = append(wrong_direction,False)

    x_vals = x_vals[~wrong_direction]
    y_vals = y_vals[~wrong_direction]

    angles = arctan2(y_vals, x_vals) * 180 / pi
    sector_angles = mod(diffAdj*diff(angles),360)

    wrong_direction = sector_angles > 180
    
# wrong_direction


    

In [77]:

    

radii = hypot(x_vals,y_vals)
sector_radii = radii[0:-1] + diff(radii)/2
plot(sector_radii)


    

    
Out[77]:





[<matplotlib.lines.Line2D at 0xb487668>]






    

In [78]:

    

angles = arctan2(y_vals, x_vals) * 180 / pi
sector_angles = mod(diffAdj*diff(angles),360)
plot(sector_angles)


    

    
Out[78]:





[<matplotlib.lines.Line2D at 0x916a4e0>]






    

In [79]:

    

plot(cumsum(sector_angles),circleFit(sector_radii*2))


    

    
Out[79]:





[<matplotlib.lines.Line2D at 0xb45f9e8>]






    

In [81]:

    

1/sum(circleFit(sector_radii*2) * sector_angles / 360)


    

    
Out[81]:





1.0508186048485129

In [ ]: