/*
elliptical-orbits.c
Simulates gravitational elliptical 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 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.
*/

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

FILE *outa, *outb;

#define points 56+1          // single point 1 orbiting a center of mass

const long double pi = 3.1415926535897932384626433833L;
const long double alpha = 137.0359991770L;
const long double center_radius = 16.55512000421622;		//sqrtl(2.0L * alpha);
const long double ipoints = points - 1.0L;
long double xx1, xx2, yy1, yy2;
long double arrow, G, x[points+1], y[points+1];




void rotate_full(long double xr1, long double yr1, long double xr2, long double yr2)
{
    long double r, x0, y0, pheta;

    r = sqrtl((xr1-xr2)*(xr1-xr2) + (yr1-yr2)*(yr1-yr2)) / 2.0L;
    if (r < center_radius) r = center_radius;
    x0 = (xr1 + xr2) / 2.0L;
    y0 = (yr1 + yr2) / 2.0L;
    pheta = atan2l(yr1 - y0, xr1 - x0) + (arrow / (G * r * sqrtl(r)));

    xx1 = x0 + cos(pheta) * r;
    yy1 = y0 + sin(pheta) * r;
    xx2 = x0 + cos(pheta + pi) * r;
    yy2 = y0 + sin(pheta + pi) * r;
}

void rotate_half(long double xr1, long double yr1, long double xr2, long double yr2)
{
    long double r, x0, y0, pheta;

    r = sqrtl((xr1-xr2)*(xr1-xr2) + (yr1-yr2)*(yr1-yr2)) / 2.0L;
    x0 = (xr1 + xr2) / 2.0L;
    y0 = (yr1 + yr2) / 2.0L;
    pheta = atan2l(yr1 - y0, xr1 - x0) + (arrow / (G * r * sqrtl(r)));

    xx1 = x0 + cos(pheta) * r;
    yy1 = y0 + sin(pheta) * r;
}

void orbit_point(long double xorbit, long double yorbit)
{
    long double circx[points+1], circy[points+1], xsum, ysum;
    int ksj, npoints = points;

    for (ksj = 2; ksj <= npoints; ksj++) {
        rotate_half(xorbit, yorbit, x[ksj], y[ksj]);
        circx[ksj] = xx1;
        circy[ksj] = yy1;
    }
    xsum = 0.0L;
    ysum = 0.0L;
    for (ksj = 2; ksj <= npoints; ksj++) {
        xsum += circx[ksj];
        ysum += circy[ksj];
    }
    xx1 = xsum / ipoints;
    yy1 = ysum / ipoints;
}




int main()
{
    unsigned int ks, ksi, ksj, ksum, kr, cnt_orbit, cycle, prt;
    unsigned int i, j, maxj, prev_counter;
    long unsigned int counter, turn_counter[99], start_count, incr, incrmax;
    long double jmax, jpoints, xsum, ysum;
    long double angle, eccentricity[99], direction, delta_period[99];
    long double xLe, yLe, xRe, yRe, ex, ey, kcurve, orbit_ey;
    long double x1, y1, sx[points+2][points+2], sy[points+2][points+2];
    long double semi_a, semi_b, r, delta_t, delta_half, half_turn[99];
    long double turn_pheta, GR, precession[99], delta_precession[99], turn_x, turn_y, periapsis, apoapsis;
    long double dp, dt, e, cyc, semi_ave;
    long double prev_x, prev_y, speed, min_speed, vx, vy, rx, ry;
	long double calc_period, calc_radius, calc_barycenter;
    double pn;

    //outa = fopen("plotfile_M8-kr12-kc2.txt", "w");
    //outb = fopen("datafile_M8-kr12-kc2.txt", "w");


    // Main parameters
    kr = 12;
    kcurve = 2.0L;          // 0 = circular, 1 = straight line, 9999... = circular
    direction = 1.0L;       // 1 anticlockwise, -1 = clockwise
    cycle = 4;              // number of half orbits to repeat
    pn = 1000.0;        //for data file



    // General parameters
    jpoints = ipoints + 1.0L;
    jmax = kr * ipoints + 1.0L;
    G = sqrtl(2.0L * jpoints / ipoints);

    //    calc_period = 16.0L*pi*alpha*sqrt(alpha)*jmax*jmax*jmax/(ipoints*ipoints*sqrt(ipoints)*sqrt(jpoints));
    calc_period = 16.0L * pi * powl(alpha, 1.5L) * powl(jmax, 3.0L) / (powl(ipoints, 2.5L) * sqrtl(jpoints));
    calc_radius = 2.0L * alpha * 2.0L * jmax * jmax / (ipoints * ipoints);
    calc_barycenter = calc_radius / jpoints;
    incrmax = (long unsigned int)(calc_period / 8192.0L);    //select required number of lines in data file

    printf("Initial Parameters:\n");
    printf("points: %d\n", points);
    printf("kr: %d\n", kr);
    printf("kcurve: %.0lf\n", (double)kcurve);
    printf("G: %.5lf\n", (double)G);
    printf("direction: %.0lf\n", (double)direction);
    printf("cycles: %.0lf\n", (double)cycle);
    printf("\nTheoretical Metrics (circular orbit):\n");
    printf("Calculated period: %.1lf\n", (double)calc_period);
    printf("Calculated radius: %.6lf\n", (double)calc_radius);
    printf("Calculated barycenter: %.5lf\n\n", (double)calc_barycenter);

/*
    fprintf(outb, "points: %d\n", points);
    fprintf(outb, "kr: %d\n", kr);
    fprintf(outb, "kcurve: %.0lf\n", (double)kcurve);
    fprintf(outb, "G: %.5lf\n", (double)G);
    fprintf(outb, "direction: %.0lf\n", (double)direction);
    fprintf(outb, "cycles: %.0lf\n", (double)cycle);
    fprintf(outb, "Calculated period: %.1lf\n", (double)calc_period);
    fprintf(outb, "Calculated radius: %.6lf\n", (double)calc_radius);
    fprintf(outb, "Calculated barycenter: %.5lf\n\n", (double)calc_barycenter);
*/

    // Initialize coordinates
    x[1] = calc_radius;
    y[1] = 0.0L;
    x[2] = 0.0L;
    y[2] = 0.0L;

    for (ks = 3; ks <= points; ks++) {
        angle = (2.0L * pi / (points - 2.0L)) * (ks - 3.0L);
        x[ks] = cos(angle) * center_radius;
        y[ks] = sin(angle) * center_radius;
    }
    x1 = xLe = xRe = ex = prev_x = x[1];
    y1 = yLe = yRe = ey = prev_y = y[1];
    maxj = 65535;
    orbit_ey = 0.0L;

 //initialize precession variables
    counter = cnt_orbit = prt = start_count = prev_counter = incr = 0;
    periapsis = min_speed = calc_radius * 2.0L;
    apoapsis = semi_ave = 0.0L;
    for (ks = 0; ks < 99; ks++) {
        turn_counter[ks] = 0;
        precession[ks] = eccentricity[ks] = delta_precession[ks] = half_turn[ks] = 0.0L;
    }




// ==================================
// Main simulation loop
// ==================================
    for (j = 0; j <= maxj; j++) {
        for (i = 0; i <= 65535; i++) {
            counter++;

            // Calculate elliptical orbiting point
            if (kcurve > 0.0L) {
                arrow = direction;
                orbit_point(xLe, yLe);
                xLe = xx1;
                yLe = yy1;

                arrow = -1.0L * direction;
                orbit_point(xRe, yRe);
                xRe = xx1;
                yRe = yy1;

                x1 = (kcurve * xLe + xRe) / (kcurve + 1.0L);
                y1 = (kcurve * yLe + yRe) / (kcurve + 1.0L);

                arrow = +1.0L;
                orbit_point(ex, ey);
                ex = xx1;
                ey = yy1;
            } else {
                arrow = +1.0L;
                orbit_point(x1, y1);
                x1 = xx1;
                y1 = yy1;
            }

            // Calculate center mass for a normal orbit
            arrow = +1.0L;
            sx[points][points] = 0.0L;
            sy[points][points] = 0.0L;
            for (ksi = 1; ksi < points; ksi++) {
                sx[ksi][ksi] = 0.0L;
                sy[ksi][ksi] = 0.0L;
                for (ksj = ksi + 1; ksj <= points; ksj++) {
                    rotate_full(x[ksi], y[ksi], x[ksj], y[ksj]);
                    sx[ksi][ksj] = xx1;
                    sy[ksi][ksj] = yy1;
                    sx[ksj][ksi] = xx2;
                    sy[ksj][ksi] = yy2;
                }
            }

            // Sum contributions
            for (ksum = 2; ksum <= points; ksum++) {
                xsum = 0.0L;
                ysum = 0.0L;
                for (ksj = 1; ksj <= points; ksj++) {
                    xsum += sx[ksum][ksj];
                    ysum += sy[ksum][ksj];
                }
                x[ksum] = xsum / ipoints;
                y[ksum] = ysum / ipoints;
            }
            x[1] = x1;
            y[1] = y1;




// ==================================
// Analysis section:
// ==================================
/*
In this test program the orbit focus is in the center and so it is not a standard elliptical orbit but a fly-by, and so apoapsis and periapsis occur twice per 360 turn.
Basic approximations; in this test model the number of cycles will be <12 and so the apoapsis will only occur in the 1st and 3rd quadrants
*/

            //detects minimum velocity
            vx = x1 - prev_x;
            vy = y1 - prev_y;
            speed = vx*vx + vy*vy;
            if (speed < min_speed){
                    min_speed = speed;
                        turn_x = x1;
                        turn_y = y1;
                        turn_pheta = atan2l(y1 - y[2], x1 - x[2]);
                        turn_counter[cnt_orbit] = counter;
                        start_count++;
            }
            prev_x = x1;
            prev_y = y1;


            //detects apoapsis and periapsis
            rx = x1 - x[2];
            ry = y1 - y[2];
            r = sqrtl(rx * rx + ry * ry);
            if (r > apoapsis)
            {
                apoapsis = r;
                angle = atan2l(y1 - y[2], x1 - x[2]);
            }
            if (r < periapsis) periapsis = r;


            //detects 180 orbit, turn occurred in 1st quadrant
            if (orbit_ey > 0.0L  && ey < 0.0L) {
                    prt = 1;
                    precession[cnt_orbit] = turn_pheta;
            }
            //detects 360 orbit, turn occurred in 3rd quadrant
            if (orbit_ey < 0.0L && ey > 0.0L) {
                    prt = 1;
                    precession[cnt_orbit] = pi + turn_pheta;
            }
            orbit_ey = ey;


            //output once every half circle
            if (prt ==1 ) {
                    half_turn[cnt_orbit] = counter - prev_counter;
                    prev_counter = counter;
                    if (cnt_orbit > 0) {
                        semi_a = apoapsis;
                        semi_b = periapsis;
                        semi_ave = semi_ave + semi_a;
                        eccentricity[cnt_orbit] = sqrtl(1.0L - (semi_b * semi_b) / (semi_a * semi_a) );

                        delta_half = (long double)half_turn[cnt_orbit - 1] ;
                        delta_t = (long double)(turn_counter[cnt_orbit] - turn_counter[cnt_orbit - 1]) ;
                        delta_period[cnt_orbit] = (pi * delta_t / delta_half) - pi;
                        delta_precession[cnt_orbit] = precession[cnt_orbit] - precession[cnt_orbit - 1];

                        printf("anomalistic point  %d   %d   %d    %.7f      %.7f    %.3f    %.3f  \n", cnt_orbit, turn_counter[cnt_orbit], counter, (double)turn_pheta,  (double)angle, (double)turn_x, (double)turn_y);
                        //fprintf(outb, "anomalistic point  %d   %d   %d    %.7f    %.5f    %.5f  \n", cnt_orbit, turn_counter[cnt_orbit], counter, (double)turn_pheta, (double)turn_x, (double)turn_y);

                        printf("period  %d    %.1f    %.1f    %.7f \n", start_count,  (double)delta_t, (double)delta_half,  (double)delta_t/(double)delta_half);
                        //fprintf(outb, "period  %d    %.1f    %.1f    %.7f \n", start_count,  (double)delta_t, (double)delta_half,  (double)delta_t/(double)delta_half);

                        printf("eccentricity %.9f     %.9f     %.9f  \n",  (double)semi_a, (double)semi_b,  (double)eccentricity[cnt_orbit]);
                        //fprintf(outb, "eccentricity %.9f     %.9f     %.9f  \n",  (double)semi_a, (double)semi_b,  (double)eccentricity[cnt_orbit]);

                        printf("precession angle %.9f     %.9f    %.9f     \n\n",  (double)precession[cnt_orbit],  (double)delta_precession[cnt_orbit], (double)delta_period[cnt_orbit]);
                        //fprintf(outb, "precession angle %.9f     %.9f    %.9f  \n\n",  (double)precession[cnt_orbit],  (double)delta_precession[cnt_orbit], (double)delta_period[cnt_orbit]);
                 }
                else {
                        printf("precession angle  %d  %d   %.7f   %.7f   %.0f \n\n",  cnt_orbit,  start_count,  (double)precession[0],  (double)angle,  (double)half_turn[0]);
                        //fprintf(outb, "precession angle  %d  %d  %.7f   %.0f   %.3f   %.3f\n\n",  cnt_orbit,  start_count,  (double)precession[0],  (double)half_turn[0],  (double)turn_x, (double)turn_y );
                 }

                prt = start_count = 0;
                periapsis = min_speed = 9999999.9L;
                apoapsis  = 0.0L;
                cnt_orbit++;
                if (cnt_orbit > cycle) j = maxj - 1;  // terminate
            }


// ==================================
// write to files
/*==================================
		if (incr > incrmax){
			incr = 0;
			fprintf(outa, "%d  %.4lf  %.4lf   %.4lf  %.4lf   %.4lf  %.4lf\n", counter, (double)x[1]/pn, (double)y[1]/pn, (double)x[2]/pn, (double)y[2]/pn, (double)ex/pn, (double)ey/pn);
		}
        incr++;
*/


        }
    }




// ==================================
// Summary output
// ==================================
    dp = dt = e = cyc = 0.0L;
    for (ks = 2; ks <= cycle; ks++) {
        printf("%d   %d   %.9f   %.9f    %.9f   %.9f  \n",  ks, turn_counter[ks] - turn_counter[ks - 1],  (double)precession[ks], (double)delta_precession[ks], (double)delta_period[ks], (double)eccentricity[ks]);
        dp = dp + delta_precession[ks];
        dt = dt + delta_period[ks];
        e = e + eccentricity[ks];
    }
    cyc = (long double)(cycle - 1);
    dt = dt / cyc;
    dp = dp / cyc;
    e = e / cyc;
    semi_ave = semi_ave / (cyc + 1.0L);
    printf("results: %.9f    %.9f   %.9f   %.7f  \n", (double)dp, (double)dt, (double)e, (double)semi_ave);

    GR = dt / (3.0L * pi / (semi_ave * (1.0L - e * e ))); //half orbit
    printf("GR comparison velocity ... %.5f\n", (double)GR );
    GR = dp / (3.0L * pi / (semi_ave * (1.0L - e * e ))); //half orbit
    printf("GR comparison apoapsis ... %.5f\n", (double)GR );



    fclose(outa);
    fclose(outb);
    scanf("%d", &j);  // Pause before exit
}