/*
Gravitational-orbital-simulation-2body.c (May 2025)
Simulates gravitational orbits using the rotating point-to-point orbital pairing model

3. Gravitational orbits from n-body rotating particle-particle orbital pairs
https://www.doi.org/10.2139/ssrn.3444571


 * Copyright 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.
*/


#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <float.h>
#include <string.h>

FILE *outa, *outb;

#define MAX_POINTS 16+1             //total number of points in simulation

// More digits for maximum precision
const long double pi = 3.1415926535897932384626433832795L;
const long double alpha = 137.0359991770L;
const long double center_radius = 16.5551200042L;


int main()
{
    int k, points, irot, jrot, incr, save_interval;
    unsigned int i, j, maxi, maxj;

    // Using aligned memory for better cache performance
    long double __attribute__((aligned(64))) x[MAX_POINTS+1];
    long double __attribute__((aligned(64))) y[MAX_POINTS+1];
    long double __attribute__((aligned(64))) sumx[MAX_POINTS+1][MAX_POINTS+1];
    long double __attribute__((aligned(64))) sumy[MAX_POINTS+1][MAX_POINTS+1];

    long double calc_period, calc_radius, calc_velocity, calc_barycenter, calc_length;
    long double jpoints, ipoints, jmax, wavelength, radius, G, imax, lambda;
    long double xsum, ysum, r, first_y,  angle;
    long double lgth_orbit, lgth_center, last_x1, last_y1;
    long double x0, y0, ang, pheta, x1, y1, x2, y2, min_distance, start_distance;
    long double minx, maxx, miny, maxy,  kr, counter, kx;
    long double dalpha = 2.0L * alpha, ralpha = sqrtl(2.0L * alpha);




    // Open files
    outa = fopen ("data-1.txt", "w");
    outb = fopen ("data-2.txt", "w");



    // Main parameters ..........................................................................................
    points = MAX_POINTS;                             // number of points in simulation
    kr = 12.0L;                           // wavelength radius coefficient (for the orbiting point), 32 is about the limit for long double precision



    // Initialize
    jpoints = (long double)points;
    ipoints = jpoints - 1.0L;
    jmax = (kr * ipoints) + 1.0L;
    wavelength = 2.0L * jmax * jmax / (ipoints * ipoints);
    G = sqrtl(2.0L * jpoints/ipoints);                                                           // for atomic orbitals G = sqrtl(2.0L*alpha);
    maxi = 65535;


    // Calculate theoretical values
    calc_period = 16.0L*pi*alpha*sqrt(alpha)*jmax*jmax*jmax/(ipoints*ipoints*sqrt(ipoints)*sqrt(jpoints));
    calc_radius = radius = 2.0L * alpha * wavelength;
    calc_velocity = sqrtl(ipoints / (radius * jpoints));
    calc_barycenter = radius/jpoints;
    calc_length = 2.0L * pi * (radius - calc_barycenter);



    // Adjust save interval to disk as required
    maxj = (int)(1.05L*(calc_period / 65535.0L));
    save_interval = (int)(calc_period / 8192.0L);



    // Print initial values
    printf("points %.0f \n", (double)jpoints);
    printf("kr %.0f \n", (double)kr);
    printf("jmax %.0f \n", (double)jmax);
    printf("wavelength %.6f \n", (double)wavelength);
    printf("G      %.6f \n\n", (double)G);
    // Print calculated values
    printf("calculated values \n");
    printf("calculated period %.7f \n", (double)calc_period);
    printf("calculated radius %.6f \n", (double)calc_radius);
    printf("calculated orbit length %.7f \n", (double)calc_length);
    printf("calculated orbit velocity %.9f \n", (double)calc_velocity);
    printf("calculated barycenter %.7f \n", (double)calc_barycenter);
    printf("jmax/jpoints ratio %.7f \n", (double)(jmax/jpoints));
    printf("save interval %d \n", save_interval);
    printf("maxj %u \n", maxj);



    // Set orbiting point
    x[1] = radius;
    y[1] = 0.0L;


    // Set center mass points equidistant from the center .................................................................
    for (k = 1; k < points; k++) {
        angle = 2.0L*pi*k/ipoints;
        x[k+1] = cosl(angle) * center_radius;
        y[k+1] = sinl(angle) * center_radius;
                        //printf("%.5f  %.5f  \n",  (double)x[k+1], (double)y[k+1]);
    }





    // Initialize variables for simulation
    incr = 0;
    first_y = 99.0L;

	lgth_orbit = 0.0L;
    last_x1 = x[1];
    last_y1 = y[1];

	minx = 99.0L;
    maxx = -99.0L;
    miny = 99.0L;
    maxy = -99.0L;




    // Start simulation
    fprintf(outa, "%.3f  %.3f\n",  (double)x[1], (double)y[1]);
    fprintf(outb, "%.3f  %.3f\n",  (double)x[2], (double)y[2]);
    printf("xy start values  %.5f     %.5f     %.5f     %.5f\n\n",  (double)x[1], (double)y[1], (double)x[2], (double)y[2]);


    for (j = 0; j <= maxj; j++) {
        for (i = 1; i <= maxi; i++) {

            // Clear sum arrays
            for (irot = 1; irot <= points; irot++) {
                for (jrot = 1; jrot <= points; jrot++) {
                    sumx[irot][jrot] = 0.0L;
                    sumy[irot][jrot] = 0.0L;
                }
            }

            // Main computation loop - this is the most intensive part
            for (irot = 1; irot < points; irot++) {
                for (jrot = irot + 1; jrot <= points; jrot++) {
                    x1 = x[irot];
                    y1 = y[irot];
                    x2 = x[jrot];
                    y2 = y[jrot];

                    // Optimize this calculation
                    long double dx = x1 - x2;
                    long double dy = y1 - y2;
                    long double dist_squared = dx*dx + dy*dy;
                    r = sqrtl(dist_squared) / 2.0L;
                    if (r < center_radius) r = center_radius;

                    // Calculate midpoint once
                    x0 = (x1 + x2) / 2.0L;
                    y0 = (y1 + y2) / 2.0L;

                    // inverse sqrt calculation
                    ang = 1.0L / (G * r * sqrtl(r));
                    pheta = atan2l(y1 - y0, x1 - x0);

                    // Precompute trig values
                    long double cos_pheta_ang = cosl(pheta + ang);
                    long double sin_pheta_ang = sinl(pheta + ang);
                    long double cos_pheta_pi_ang = cosl(pheta + pi + ang);
                    long double sin_pheta_pi_ang = sinl(pheta + pi + ang);

                    sumx[irot][jrot] = x0 + cos_pheta_ang * r;
                    sumy[irot][jrot] = y0 + sin_pheta_ang * r;
                    sumx[jrot][irot] = x0 + cos_pheta_pi_ang * r;
                    sumy[jrot][irot] = y0 + sin_pheta_pi_ang * r;
                }
            }

            // Update positions
            for (irot = 1; irot <= points; irot++) {
                xsum = 0.0L;
                ysum = 0.0L;

                for (jrot = 1; jrot <= points; jrot++) {
                    xsum += sumx[irot][jrot];
                    ysum += sumy[irot][jrot];
                }

                x[irot] = xsum / ipoints;
                y[irot] = ysum / ipoints;
            }



            // Update length
            lgth_orbit += sqrtl( (x[1]-last_x1)*(x[1]-last_x1) + (y[1]-last_y1)*(y[1]-last_y1) );
            last_x1 = x[1];
            last_y1 = y[1];

            // Update min/max for barycenter calculation
            if (x[1] > maxx) maxx = x[1];
            if (x[1] < minx) minx = x[1];
            if (y[1] > maxy) maxy = y[1];
            if (y[1] < miny) miny = y[1];



            // End of 1 orbit check
            if (first_y <= 0.0L && y[1] > 0.0L) {
                printf("simulation results \n");
                counter = (j*65535+i)*1.0L;
                printf("simulation period %.0f      %.9f\n", (double)counter, (double)(calc_period/counter) );
                printf("simulation orbit length %.9f      %.11f\n", (double)lgth_orbit,  (double)(calc_length/lgth_orbit) );
                printf("simulation orbit radius %.9f\n", (double)lgth_orbit/(2.0L*pi));
                printf("simulation barycenter x = %.6f    y = %.6f \n", (double)((minx+maxx)/2.0L), (double)((miny+maxy)/2.0L));
                printf("xy end values %.5f     %.5f     %.5f     %.5f\n\n",  (double)x[1], (double)y[1], (double)x[2], (double)y[2]);
                printf("%.0f  \n", j*1.0);

                // this will end the simulation, to continue, comment out
                i = 65535;
                j = 65535;
            }
            first_y = y[1];



            // Save data for graphs at specified intervals
            if (incr >= save_interval) {
                fprintf(outa, "%.3f       %.3f\n", (double)x[1], (double)y[1]);
                fprintf(outb, "%.3f       %.3f\n", (double)x[2], (double)y[2]);
                incr = 0;
            }
            incr++;

        }
    }



    // Close all files
    fclose(outa);
    fclose(outb);
    printf("end \n");
    scanf("%d", &j);        //keeps the window open

    return 0;
}