/*
Gravitational-elliptical-orbits.cpp (Nov 2025)
Program to simulate elliptical orbits using rotating point-to-point orbital pairing model
Copyright Malcolm J. Macleod (malcolm@codingthecosmos.com)


Main parameters
#define	points     //total number of points; center mass + 1 orbiting point
kr                      //a radius co-efficient
kcurve               //0 = circular, 1 = straight line, 9999... = circular
direction            //1 anticlockwise, -1 = clockwise
cycle          //number of orbits

 * 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 32+1          //for single point orbiting a center of mass


const long double pi=3.1415926535897932384626433833L;
const long double alpha=137.0359991770L;
const long double center_radius = sqrtl(2.0L*137.0359991770L);
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 = sqrt((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 = sqrt((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 = xsum + circx[ksj];
			ysum = ysum + circy[ksj];
		}
		xx1 = xsum/ipoints;
		yy1 = ysum/ipoints;

}






int main()
{
unsigned int ks, ksi, ksj, ksum, kr, cnt_orbit, cycle, ellipse;
unsigned int i, j, maxj, incr, incrmax;
long unsigned int counter, print_counter, prev_counter;
long double jmax, jpoints, xsum, ysum, bary_x, precession;
long double angle, pheta, prev_pheta, eccentricity, direction, delta_period, delta_radius, orbital_period, anomal_pheta;
long double xLe, yLe, xRe, yRe, ex, ey, kcurve, orbit_y1, orbit_ey, orbit_once;
long double x1, y1, sx[points+2][points+2], sy[points+2][points+2];
long double min_x, min_y, min_ex, semi_a, semi_b, r, x0, y0;
long double lgth1, lgth2, lgth3, lx1, lx2, lx3, ly1, ly2, ly3;
double pn;








//	initialize parameters
	//outa = fopen("e108x12x20.txt", "w");
    //outb = fopen("e2.txt", "w");
	kr = 12;
	kcurve = 108.0L;      //0 = circular, 1 = straight line, 9999... = circular
    direction = 1.0L;          //(1 anticlockwise, -1 = clockwise)
    cycle = 2;         //number of orbits to repeat
    pn = 10.0;    //data file reduced numerical size for plotting








//	general parameters
	maxj = 65534/2;
	incrmax = 65536/2;
    jpoints = ipoints + 1.0L;
	jmax = kr*ipoints + 1.0L;
    G = sqrtl(2.0L * jpoints/ipoints);

    long double calc_period = 16.0L * pi * powl(alpha, 1.5L) * powl(jmax, 3.0L) / (powl(ipoints, 2.5L) * sqrtl(jpoints));
    long double calc_radius = 2.0L * alpha * 2.0L * jmax * jmax / (ipoints * ipoints);
    long double calc_barycenter = calc_radius / jpoints;

    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("\nTheoretical Metrics:\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);
    printf("Planck units conversion (i*lp = Schwarzschild radius):\n");
    printf("Period (Planck time): %.1lf tp\n", (double)(calc_period * ipoints));
    printf("Radius (Planck length): %.6lf lp\n", (double)(calc_radius * ipoints));
    printf("Barycenter (Planck length): %.5lf lp\n\n", (double)(calc_barycenter * ipoints));



//	initialize co-ordinates
	x[1] = calc_radius;          //main point
	y[1] = 0.0L;
	x[2] = 0.0L;                 //center of center mass
	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;
	}
	for (ks=1; ks<=points; ks++) {
		//printf("%2i  %.1f   %.1f\n", ks, (double)x[ks], (double)y[ks]);
	}



	incr = incrmax + 1;
	counter = prev_counter = cnt_orbit = ellipse = 0;
	prev_pheta = min_ex = orbit_ey = min_x = min_y = 0.0L;
    lgth1 = lgth2 = lgth3 = lx1 = lx2 = lx3 = ly1 = ly2 = ly3 = 0.0L;
	x1 = xLe = xRe = ex = x[1];
	y1 = yLe = yRe = ey = y[1];
    semi_a = semi_b = calc_radius;




//calculate
	for (j=0; j<=maxj; j++) {
	for (i=0; i<=65535; i++)
	{

    counter++;
//	calculate elliptical orbiting point separately
	if (kcurve>0.0L)
	{
		// anti-clockwise turn
            arrow = direction;
                orbit_point(xLe, yLe);
                xLe = xx1;
                yLe = yy1;

		// clock-wise turn
            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);

		//reference point, doesnt influence center mass
            arrow = +1.0L;
                orbit_point(ex, ey);
                ex = xx1;
                ey = yy1;
	}

	else
	{
//		calculate circular
            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
		for (ksum=2; ksum<=points; ksum++) {
			xsum = 0.0L;
			ysum = 0.0L;
			for (ksj=1; ksj<=points; ksj++) {
				xsum = xsum + sx[ksum][ksj];
				ysum = ysum + sy[ksum][ksj];
			}
			x[ksum] = xsum/ipoints;
			y[ksum] = ysum/ipoints;
		}
		//update with new x1, y1
		x[1] = x1;
		y[1] = y1;



        //distance traveled (frame dragging analysis)
		lgth1 = lgth1 + sqrtl( (x1-lx1)*(x1-lx1) + (y1-ly1)*(y1-ly1) );
		lx1 = x1;
		ly1 = y1;
		lgth2 = lgth2 + sqrtl( (x[2]-lx2)*(x[2]-lx2) + (y[2]-ly2)*(y[2]-ly2) );
		lx2 = x[2];
		ly2 = y[2];
		lgth3 = lgth3 + sqrtl( (x[3]-lx3)*(x[3]-lx3) + (y[3]-ly3)*(y[3]-ly3) );
		lx3 = x[3];
		ly3 = y[3];




//      Update maximum and minimum
        if (x1 < min_x) {
                min_x = x1;
                min_y = y1;
        }
        r = sqrt((x1-x[2])*(x1-x[2]) + (y1-y[2])*(y1-y[2]));
        if (r < semi_b) semi_b = r;
        if (ex < min_ex) min_ex = ex;
        if (x1 > orbit_once) orbit_once = x1;


//		test for circular orbit
        if (orbit_ey < 0.0L) {
                if (ey > 0.0L) {
                        orbit_once = min_ex;
                        ellipse = 1;
                        print_counter = counter;
                        bary_x = (calc_radius + min_ex) /2.0L;
                        printf("circular orbit... %d   %.0f   %.7f    %.7f    %.7f \n\n", j, counter*1.0, (double)ex, (double)ey, (double)bary_x);
                }
        }
        orbit_ey = ey;



 //		test for anomalistic orbit (where orbit turns around)
        if (ellipse > 0) {
            if (x1 < orbit_once) {
                    cnt_orbit++;
                    orbital_period = 1.0L*(counter - prev_counter);
                    printf("orbital period (elliptical) ... %d   %d   %.0f \n", cnt_orbit, j, (double)orbital_period);
                    r = sqrt((x1-min_x)*(x1-min_x) + (y1-min_y)*(y1-min_y))/2.0L;
                    x0 = (x1 + min_x) / 2.0L;
                    y0 = (y1 + min_y) / 2.0L;
                    pheta = atan2l(y1 - y0, x1 - x0);
                    printf("anomalistic point  %.11f    %.7f   %.7f   %.7f   %.7f \n",  (double)(pheta),  (double)x1, (double)y1, (double)min_x, (double)min_y );

                    //orbital analysis
                    if (cnt_orbit <= 3){    //cnt_orbit <= ... select number of orbits prints
                            delta_radius = cbrtl(jpoints * orbital_period * orbital_period / (4.0L * pi * pi * ipoints));
                            delta_period = calc_period * (2.0L*pi + pheta)/(2.0L*pi);
                            printf("equivalent radius  %.3f  %.9f \n",  (double)delta_radius, (double)(delta_radius/calc_radius));
                            printf("equivalent period  %.3f  %.9f \n", (double)delta_period, (double)(orbital_period/delta_period));
                            printf("distance traveled %.7f   %.7f   %.7f \n", (double)lgth1, (double)lgth2, (double)lgth3);

                            eccentricity = sqrtl(1.0L - (semi_b/semi_a) * (semi_b/semi_a) );
                            printf("axis a, b %.9f     %.9f     %.9f  \n",  (double)semi_a, (double)semi_b, (double)(semi_b/semi_a));

                            anomal_pheta = pheta - prev_pheta;
                            //precession = 6.0L * pi * semi_a / (ipoints * semi_b * semi_b); //= 6.0L * pi / (semi_a * (1.0L - eccentricity * eccentricity) );
                            printf("precession angle %.9f    %.9f    %.9f \n\n", (double)anomal_pheta, (double)(anomal_pheta*ipoints), (double)eccentricity );
                    }

                    //reset
                    prev_counter = counter;
                    min_x = min_y = semi_b = calc_radius;
                    orbit_once = min_ex;
                    prev_pheta = pheta;
                    ellipse = 0;
                    lgth1 = lgth2 = lgth3 = lx1 = lx2 = lx3 = ly1 = ly2 = ly3 = 0.0L;
                    if (cnt_orbit >= cycle)  j = maxj - 1;          // terminate
                }
        }


/*
//		write to file at intervals
		if (incr > incrmax){
			incr = 0;
			fprintf(outa, "%.4lf  %.4lf   %.4lf  %.4lf   %.4lf  %.4lf\n", (double)x[1]/pn, (double)y[1]/pn, (double)x[2]/pn, (double)y[2]/pn, (double)ex/pn, (double)ey/pn);
			//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++;

//    internal orbit check
        if (counter > print_counter){
            if (counter < print_counter + 8192) {
                //fprintf(outb, "%d  %.6lf   %.6lf \n", counter, (double)x[1]/pn, (double)y[1]/pn);
                //fprintf(outb, "%.5lf   %.4lf \n", (double)x[1], (double)y[1]);
            }
        }
*/


	}	//for i
	}	//for j



	printf("end ... %.0f\n", counter*1.0);
	fclose(outa);
	fclose(outb);
	scanf("%d", &j);
}





/*
Sample gnuplot to animate results

--------------------------------------------------------------------
reset
set terminal gif animate delay 1

set title "orbit "
set output 'C:\files\Hatom.gif'
f = 'C:\files\data32.txt'

pxL = 18
pxR = 6
pyU = 89
pyD = 4
incr = 256
ed = 62743-incr

set pointsize 1
show title
set parametric
j=0.0
unset key


st = 0

set xrange [-pxL:pxR]
set yrange [-pyD:pyU]


do for[j=st:ed:incr]{

plot f u 1:2 every ::j::j w p pt 7 ps 1 lc rgb "dark-blue", f u 1:2 every ::0::j w l lc rgb "dark-blue", f u 3:4  every ::j::j w p pt 7 ps 1 lc rgb "red", f u 3:4 every ::0::j w l lc rgb "red"

}

set output

---------------------------------------------------------------------------------
#   load 'C:\files\plot-Hatom.p'



*/