#A Python program to generate the coordinates of first 100 polygons

import numpy as np
import math as m

#Get centre of a n-gon from two vertices
def get_centre2(p1, p2, n):
    x1 = p1[0]
    y1 = p1[1]
    x2 = p2[0]
    y2 = p2[1]
    
    ax = (x2 - x1)/2
    ay = (y2 - y1)/2
    
    bx1 =  - ay / m.tan(m.pi/n)
    by1 = ax / m.tan(m.pi/n)
    
    bx2 = ay / m.tan(m.pi/n)
    by2 = - ax / m.tan(m.pi/n)
    
    if ax*by1 - ay*bx1 < 0:
        bx = bx1
        by = by1
    else:
        bx = bx2
        by = by2
    
    cx = p1[0] + ax + bx
    cy = p1[1] + ay + by
    
    return [cx, cy]

#get all vertex coordinates
def get_points(p1, p2, n):
    centre = get_centre2(p1,p2,n)

    
    points = []
    points.append(p1)
    points.append(p2)
    

    
    for i in np.arange(2,n,1):
        new_point_x = np.cos(i * 2*m.pi/n) * (p1[0] - centre[0]) + np.sin(i * 2*m.pi/n) * (p1[1] - centre[1]) + centre[0]
        new_point_y = -np.sin(i * 2*m.pi/n) * (p1[0] - centre[0]) + np.cos(i * 2*m.pi/n) * (p1[1] - centre[1]) + centre[1]
        
        points.append([new_point_x, new_point_y])

    return [centre, points]

#create candidate edges for next n-gon in series
def create_candidates(array):
    new_array = []
    for i in array:
        new_array.append(i)
    new_array.append(array[0])
    candidates = []
    for i in np.arange(0, len(array) - 1, 1):
        candidate_pair = []
        candidate_pair.append(new_array[i+2])
        candidate_pair.append(new_array[i+1])
        candidates.append(candidate_pair)
    
    return candidates

#initialise polygons, create first two (triangle, and square)
polygons = []
polygons.append([[0, 0], [1, 0], [0.5, 3**0.5/2]])
polygons.append([[0, 0], [1, 0], [1, -1], [0, -1]])

#initialise indices; a(n)
index = []
index.append(1)


all_polygons = [Polygon(j) for j in polygons]
union = unary_union(all_polygons)


#run loop
for i in np.arange(5,100,1):
    candidate_vertices = create_candidates(polygons[i-4])[int(np.floor((i-4)*4/10)):]
    vertex_index = 2+int(np.floor((i-4)*4/10))#makes it a bit faster

    for candidate in candidate_vertices:        
        candidate_polygon_centre0, candidate_polygon_points0 = get_points(candidate[0],candidate[1],i)
        candidate_polygon0 = Polygon(candidate_polygon_points0)
        
        if np.abs(unary_union([union, candidate_polygon0]).area - (union.area + candidate_polygon0.area)) < 10**(-5):
            new_polygon_points = candidate_polygon_points0
            polygons.append(new_polygon_points)
            index.append(vertex_index)
            union = unary_union([union,Polygon(new_polygon_points)])
            break
        else:
            vertex_index += 1
            continue