"""
Gravitational-orbital-modeling
Newton-vs-Orbital_Orbital.py

Orbital derived 3 body orbit for comparison (see Newton-vs-Orbital_Newton.py)

Copyright 2025 Malcolm J. Macleod (malcolm@codingthecosmos.com)

* This software is free to use, modify, and distribute under creative commons 
 * https://creativecommons.org/licenses/by-nc-sa/4.0/ (CC CC BY-NC-SA 4.0). 
 * Any derivative works or redistributions must include appropriate credit to 
 * Malcolm Macleod and this permission notice shall
 * be included in all copies or substantial portions of the Software.
 
"""

import math

MAX_POINTS = 3
outa = open("data-3b-long.txt", "w")
outb = open("data-3b-plot.txt", "w")

pi = math.pi
alpha = 137.0359991770
center_radius = 1.0*math.sqrt(2.0*alpha)




def process_points(points_x, points_y, G):
    """
    Process points through rotation and averaging
    :param points_x: List of x coordinates
    :param points_y: List of y coordinates
    :return: Updated (x, y) lists
    """
    n = len(points_x)
    if n < 2:
        return points_x, points_y
    
    # Storage for rotated positions
    new_x = [[0.0] * (n-1) for _ in range(n)]
    new_y = [[0.0] * (n-1) for _ in range(n)]
    counters = [0] * n
    
    # Process all unique point pairs
    for i in range(n):
        for j in range(i + 1, n):
            # Calculate midpoint
            x0 = (points_x[i] + points_x[j]) / 2.0
            y0 = (points_y[i] + points_y[j]) / 2.0
            
            # Calculate half-distance with minimum radius check
            dx = points_x[i] - points_x[j]
            dy = points_y[i] - points_y[j]
            r = math.sqrt(dx*dx + dy*dy) / 2.0
            if r < center_radius:
                r = center_radius
            
            # Calculate rotation angle
            ang = 1.0 / (G * r * math.sqrt(r))
            pheta = math.atan2(points_y[i] - y0, points_x[i] - x0)
            
            # Calculate new positions
            new_x_i = x0 + math.cos(pheta + ang) * r
            new_y_i = y0 + math.sin(pheta + ang) * r
            new_x_j = x0 + math.cos(pheta + ang + math.pi) * r
            new_y_j = y0 + math.sin(pheta + ang + math.pi) * r
            
            # Store new positions
            new_x[i][counters[i]] = new_x_i
            new_y[i][counters[i]] = new_y_i
            counters[i] += 1
            
            new_x[j][counters[j]] = new_x_j
            new_y[j][counters[j]] = new_y_j
            counters[j] += 1
    
    # Calculate average positions
    updated_x = []
    updated_y = []
    for i in range(n):
        avg_x = sum(new_x[i]) / (n - 1)
        avg_y = sum(new_y[i]) / (n - 1)
        updated_x.append(avg_x)
        updated_y.append(avg_y)
    
    return updated_x, updated_y





def main():
    points = MAX_POINTS
    jpoints = float(points)
    G = math.sqrt(2.0 * jpoints / (jpoints - 1.0))
    st = 99

    #save_interval = 2**6
    maxj = 1118213 + 256     #save_interval
    

   
    # Initialize positions
    x = [0.0] * (points + 1)
    y = [0.0] * (points + 1)
    r0 = 2.0*alpha
    x[1] = r0*6.355*2.0#6.3675     360
    y[1] = 0.0
    x[2] = math.cos(pi*2.0/3.0)*r0
    y[2] = math.sin(pi*2.0/3.0)*r0
    x[3] = math.cos(pi*4.0/3.0)*r0
    y[3] = math.sin(pi*4.0/3.0)*r0
    
    # Print initial parameters
    print(f"points: {points}")
    print(f"maxj: {maxj:.0f}")
    print(f"x1: {x[1]:.5f}")
    
    
    # Initialize simulation variables
    #outa.write(f"{x[1]:.3f} {y[1]:.3f}\n")
    #outb.write(f"{x[2]:.3f} {y[2]:.3f}\n")

    # Simulation loop
    for j in range(1, maxj + 1):
        # Extract current points (excluding index 0)
        current_x = x[1:points+1]
        current_y = y[1:points+1]
        
        # Process points through rotation and averaging
        updated_x, updated_y = process_points(current_x, current_y, G)
        
        # Update positions
        for i in range(1, points + 1):
            x[i] = updated_x[i-1]
            y[i] = updated_y[i-1]
            
        # check for orbit completion
        if (st < 0):
            if (y[1] > 0):
                print(f"period: {j:.0f}")
        st = y[1]
        
        # Save data at specified intervals
        outa.write(f"{j:.0f} {x[1]:.3f} {y[1]:.3f} {x[2]:.3f} {y[2]:.3f} {x[3]:.3f} {y[3]:.3f}\n")  #orbiting point
        #if j % save_interval == 0:
        #    outa.write(f"{j:.0f} {x[1]:.3f} {y[1]:.3f} {x[2]:.3f} {y[2]:.3f} {x[3]:.3f} {y[3]:.3f}\n")  #orbiting point
        #    outb.write(f"{x[1]:.3f} {y[1]:.3f}\n") 
    
    outa.close()
    outb.close()
    print("end")


if __name__ == "__main__":
    main()