A tale of two rabbits

import pygame
import math

# Initialize Pygame
pygame.init()

# Set up the display
WIDTH = 800
HEIGHT = 600
screen = pygame.display.set_mode((WIDTH, HEIGHT))
pygame.display.set_caption("A tale of two rabbits")

# Colors
WHITE = (255, 255, 255)
BLACK = (0, 0, 0)
RED = (255, 0, 0)
BLUE = (0, 0, 255)

def interpolate_color(color1, color2, factor):
    r = int(color1[0] + (color2[0] - color1[0]) * factor)
    g = int(color1[1] + (color2[1] - color1[1]) * factor)
    b = int(color1[2] + (color2[2] - color1[2]) * factor)
    return (r, g, b)

# Animation parameters
num_sides = 24
radius = 150
growth_rate = 0.5
rotation_angle = 0

# Function to calculate polygon points
def get_polygon_points(center_x, center_y, radius, num_sides, rotation):
    points = []
    for i in range(num_sides):
        angle = rotation + (2 * math.pi * i / num_sides)
        x = center_x + radius * math.cos(angle)
        y = center_y + radius * math.sin(angle)
        points.append((int(x), int(y)))
    return points

def relax_points(points, step_size):
    new_points = points.copy()
    # center_x = sum([point[0] for point in points]) / len(points)
    # center_y = sum([point[1] for point in points]) / len(points)
    # target_center_x = WIDTH // 2
    # target_center_y = HEIGHT // 2
    for i in range(1, len(points) - 1):
        x, y = points[i]
        x_next, y_next = points[i + 1]
        x_prev, y_prev = points[i - 1]
        # Place the point in the middle of x_prev and x_next, but leave the curvature intact
        dot_product = (x - x_prev) * (x_next - x_prev) + (y - y_prev) * (y_next - y_prev)
        alpha = 0.5 - dot_product / ((x_next - x_prev) ** 2 + (y_next - y_prev) ** 2)
        x_new = x + alpha * (x_next - x_prev) * step_size
        y_new = y + alpha * (y_next - y_prev) * step_size
        # x_new = x_new + (target_center_x - center_x) * step_size
        # y_new = y_new + (target_center_y - center_y) * step_size
        new_points[i] = (x_new, y_new)

    return new_points

# Maximize polygon area subject to perimeter + len(x0,x1) = 10
def gradient_descent(points, permimeter_penalty, step_size):
    new_points = points.copy()
    perimeter = 0
    for i in range(len(points)):
        x, y = points[i]
        x_next, y_next = points[(i + 1) % len(points)]
        length = math.sqrt((x_next - x) ** 2 + (y_next - y) ** 2)
        if i != len(points) - 1:
            perimeter += length*2
        else:
            perimeter += length
    for i in range(len(points)):
        x, y = points[i]
        x_next, y_next = points[(i + 1) % len(points)]
        x_prev, y_prev = points[(i - 1) % len(points)]
        d_next = math.sqrt((x_next - x) ** 2 + (y_next - y) ** 2)
        d_prev = math.sqrt((x_prev - x) ** 2 + (y_prev - y) ** 2)
        prev_multiplier = 1 if i != 0 else 0.5
        next_multiplier = 1 if i != len(points) - 1 else 0.5
        grad_perimeter_x = (x - x_prev) / d_prev * prev_multiplier + (x - x_next) / d_next * next_multiplier
        grad_perimeter_y = (y - y_prev) / d_prev * prev_multiplier + (y - y_next) / d_next * next_multiplier
        grad_perimeter_sign = 1 if perimeter < 2000 else -1
        area_grad_x = (y_next - y_prev)
        area_grad_y = (x_prev - x_next)
        x_new = x + step_size * (
                area_grad_x + permimeter_penalty * grad_perimeter_x * grad_perimeter_sign)
        y_new = y + step_size * (
                area_grad_y + permimeter_penalty * grad_perimeter_y * grad_perimeter_sign)
        new_points[i] = (x_new, y_new)
    return new_points

def compute_symmetry_points(points):
    """Symmetry with respect to (points[0], points[len(points)-1])"""
    new_points = points.copy()
    x0, y0 = points[0]
    x1, y1 = points[len(points) - 1]

    # Vector of the line
    dx = x1 - x0
    dy = y1 - y0
    # Length squared of the line
    l2 = dx*dx + dy*dy

    for i in range(len(points)):
        x, y = points[i]
        # Vector from line start to point
        px = x - x0
        py = y - y0
        # Project point onto line
        dot = (px*dx + py*dy) / l2
        # Projection point
        proj_x = x0 + dot * dx
        proj_y = y0 + dot * dy
        # Reflect point across projection
        new_x = 2 * proj_x - x
        new_y = 2 * proj_y - y
        new_points[i] = (new_x, new_y)

    return new_points


# Main game loop
running = True
clock = pygame.time.Clock()

points = get_polygon_points(WIDTH // 4, HEIGHT // 2, radius, num_sides, rotation_angle)
step_size = 0.005

while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

    # Clear screen
    screen.fill(BLACK)

    # Calculate polygon points
    center_x = WIDTH // 2
    center_y = HEIGHT // 2
    points = gradient_descent(points, permimeter_penalty=500, step_size=step_size)
    points = relax_points(points, 0.1)
    step_size *= 0.999
    symmetry_points = compute_symmetry_points(points)

    # Draw polygon
    pygame.draw.polygon(screen, WHITE, points, 2)
    pygame.draw.polygon(screen, WHITE, symmetry_points, 2)

    # Draw corner points with gradient from red to blue
    for i, point in enumerate(points):
        color = interpolate_color(RED, BLUE, i / (len(points) - 1))
        pygame.draw.circle(screen, color, point, 5)
        pygame.draw.circle(screen, color, symmetry_points[i], 5)

    # Update display
    pygame.display.flip()

    # Control frame rate
    clock.tick(60)

pygame.quit()