/*
2-photon-gravitational-orbital-simulation.c (Nov 2025)
Program to simulate atomic orbital transitions using the gravity simulator
(rotating point-to-point n-body orbital pairing)

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

4. Geometrical origins of quantization in H atom electron transitions
https://www.doi.org/10.2139/ssrn.3703266

 * 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 *xy_out, *data_out;

#define	mpoints 67      //total_counter number points j = mpoints-1; center+1+1+1

const long double pi = 3.1415926535897932384626433832795L;
const long double alpha = 137.035999177L;




int main()
{
unsigned int k, pt, irot, jrot, points, nlevel, max_nlevel;
unsigned int i, j, maxi, maxj, prt, pk, save_interval;
long unsigned int counter, full_counter, Nsamples, prev_counter;

long double x1, y1, x0, y0, ipoints, tsim, jpoints, kr, jmax;
long double r, r0, start_r0, a0, radius, center_radius, spiral_radius;
long double rincr, lambda, rplus, length, orbital_length, lorbital_n1;
long double G, first_y, fh[5], dt_sim, fH, c, le, lp, torbital, torbital_n1 ;
long double n_theory, n_theory_old, n_spiral, n_spiral_old, n, ns, n_quantum;
long double orbital_counter, ralpha;
long double pheta, pheta_old, pheta180, pheta_orbital, angle, n_angle;
long double minx, maxx, miny, maxy, bary_x, bary_y;
long double *x, *y, *xbak, *ybak;
long double** sumx = (long double**)malloc(mpoints * sizeof(long double));
long double** sumy = (long double**)malloc(mpoints * sizeof(long double));



//////
        //xy_out = fopen("xy-data-m66-n2.txt", "w");
        data_out = fopen("simulation-data-m66-n4.txt", "w");
        //m=total mass, n=max_nlevel (from n=1 to n=max_nlevel

        nlevel = 1;
        max_nlevel = 4;
        rplus = 3.5L;
        kr = 1.0L;                                  //usually set at 1.0
        save_interval = (2 * max_nlevel * max_nlevel) * 2;
        // the radius shows a slight reduction in radius during orbit and so rplus is a correction factor
        // its value will depend on the center mass, n_spiral and n_spiral_old used to tune so radius has an integer value at prescribed angles
        // for example i = 3, at rplus = 1.3333333L; nlevel = 4.000000967   4.000000371



		x = (long double*)malloc(mpoints * sizeof(long double));
		y = (long double*)malloc(mpoints * sizeof(long double));
		for (i = 0; i < mpoints; i++) {
			sumx[i] = (long double*)malloc(mpoints * sizeof(long double));
		}
		for (i = 0; i < mpoints; i++) {
			sumy[i] = (long double*)malloc(mpoints * sizeof(long double));
		}

        maxj = 65534;
        maxi = 65535;

		jpoints = (long double)mpoints - 1.0L;
        ipoints = mpoints - 2.0L;
        points = mpoints - 1;
        jmax = (kr * ipoints) + 1.0L;
        ralpha = 2.0L * alpha;
        G = center_radius = sqrtl(ralpha);
        rincr = 1.0L/(2.0L*pi*ralpha);

        lambda = 2.0L * jmax * jmax / (ipoints * ipoints); //kr = 1.0 then lambda = 2.0L*jpoints*jpoints/(ipoints*ipoints);
        a0 = ralpha * lambda;
        r0 = (ralpha + (4.5L * rincr)) * lambda;
        printf("initial = %.0lf   %.14lf   %.9lf   %.9lf  %.3lf \n", (double)(jpoints), (double)lambda, (double)a0, (double)r0, (double)rplus);
        printf("radius... %.9lf  \n", (double)(r0 - a0));


//  ___________________________________________________________
//

    //reference data, experimental values
    fh[0] = 0.0L;
    fh[1] = 0.0L;
    fh[2] = (1415879.869694778);
    fh[3] = (3775674.709988457);
    fh[4] = (7079388.330994962);
    printf("H freq  = %.0lf   %.0lf   %.0lf \n", (double)fh[2], (double)fh[3], (double)fh[4] );

    fh[2] = (1887839.826259705);
    fh[3] = (4247634.048737014);
    fh[4] = (7551347.553061293);
    printf("H freq  = %.0lf   %.0lf   %.0lf \n\n", (double)fh[2], (double)fh[3], (double)fh[4] );

    le = 2.42631023538e-12;
    lp = 1.32140985360e-15;
    c = 299792458.0L;



//	set co-ordinates
	x[1] = r0;          //orbiting point
	y[1] = 0.0L;
	x[2] = 0.0L;        //center point of center mass
	y[2] = 0.0L;
	pt = points-2;      //circular pattern around center point
	for (k=1; k<=pt; k++) {
		angle = 2.0L*pi*(double)k/(double)pt;
		x[k+2] = cos(angle) * center_radius;
		y[k+2] = sin(angle) * center_radius;
	}


//	set parameters to start
    counter = prt = Nsamples = full_counter = prev_counter = 0;
	pheta_orbital = length = n_spiral_old = n_theory_old = pheta_old = orbital_counter = n_quantum = 0.0L;
    first_y = 99.0L;
    n_spiral = n_theory = orbital_length = 1.0L;
	minx = miny = 99.0L;
	maxx = maxy =-99.0L;





//====================================================
// SIMULATION BEGINS
//====================================================
    for (j = 0; j <= maxj; j++) {
        for (i = 1; i <= maxi; i++) {
        counter++;


//  Initialize accumulators for center mass points and orbiting point
    long double xsum_center[mpoints] = {0.0L};
    long double ysum_center[mpoints] = {0.0L};
    long double xsum_orbit = 0.0L, ysum_orbit = 0.0L;
    int orbit_pair_count = 0;

//  Main pairwise computation loop — preserves original precision
    for (irot = 1; irot < points; irot++) {
        for (jrot = irot + 1; jrot <= points; jrot++) {

            long double x1_val = x[irot];
            long double y1_val = y[irot];
            long double x2_val = x[jrot];
            long double y2_val = y[jrot];

            long double dx = x1_val - x2_val;
            long double dy = y1_val - y2_val;
            r = sqrtl(dx*dx + dy*dy) / 2.0L;
            if (r < center_radius) r = center_radius;

            x0 = (x1_val + x2_val) / 2.0L;
            y0 = (y1_val + y2_val) / 2.0L;

            angle = 1.0L / (G * r * sqrtl(r));
            pheta = atan2l(y1_val - y0, x1_val - x0);

            long double new_x1 = x0 + cosl(pheta + angle) * r;
            long double new_y1 = y0 + sinl(pheta + angle) * r;
            long double new_x2 = x0 + cosl(pheta + pi + angle) * r;
            long double new_y2 = y0 + sinl(pheta + pi + angle) * r;

        //  Accumulate for center mass points
            xsum_center[irot] += new_x1;
            ysum_center[irot] += new_y1;
            xsum_center[jrot] += new_x2;
            ysum_center[jrot] += new_y2;

        // Special: if one of the points is the orbiting point (index 1), accumulate its new position
            if (irot == 1) {
                xsum_orbit += new_x1;
                ysum_orbit += new_y1;
                orbit_pair_count++;
            }
            if (jrot == 1) {
                xsum_orbit += new_x2;
                ysum_orbit += new_y2;
                orbit_pair_count++;
            }
        }
    }


//      Update center mass points
        for (k = 2; k <= points; k++) {
            x[k] = xsum_center[k] / ipoints;
            y[k] = ysum_center[k] / ipoints;
        }

//      Update orbiting point (index 1)
            x1 = xsum_orbit / orbit_pair_count;  // should be ipoints, but let's be safe
            y1 = ysum_orbit / orbit_pair_count;


//      update path length and angle
        length += sqrtl((x1 - x[1])*(x1 - x[1]) + (y1 - y[1])*(y1 - y[1]));


//      reset center points around middle to keep the proton the same size
        pt = points - 2;
        for (k=1; k<=pt; k++) {
            angle = 2.0L*pi*(double)k/(double)pt;
            x[k+2] = cos(angle) + x[2];
            y[k+2] = sin(angle) + y[2];
        }


//      n level parameters
		if (nlevel==1) {
        //====================================================
        // ORBITAL PHASE
        //====================================================
            if (x1>maxx) maxx = x1;
			if (x1<minx) minx = x1;
			if (y1>maxy) maxy = y1;
			if (y1<miny) miny = y1;
            x[1] = x1;
            y[1] = y1;
            if (first_y <= 0.0L && y1 > 0.0L) {
                radius = sqrtl((x1 - x[2])*(x1 - x[2]) + (y1 - y[2])*(y1 - y[2]));
                printf("radius... %.9lf  %.0lf  \n", (double)(radius - a0), counter*1.0);
                orbital_counter = 2.0L* pi * ralpha * ralpha * lambda;
                orbital_length = 2.0L* pi * (a0 - a0 / jpoints);
                torbital_n1 = orbital_counter;
                lorbital_n1 = orbital_length;
                bary_x = (minx+maxx)/2.0L;
                bary_y = (miny+maxy)/2.0L;
				nlevel = 2;
                ns = nlevel * 1.0L;
                pheta_orbital = 360.0L;
                printf("nlevel = 1 ... %.0lf   %.9lf   %.9lf   %.9lf   \n\n", counter * 1.0,  (double)(torbital_n1),  (double)(torbital_n1/lambda),  (double)(lorbital_n1) );
                // reset to start position
                length = 0.0L;
                counter = 0;
                prev_counter = 0;
                r0 = (2.0L*alpha + (rplus * rincr)) * lambda;
                x[1] = r0;
                y[1] = 0.0L;
                x[2] = 0.0L;
                y[2] = 0.0L;
                pt = points-2;
                for (k=1; k<=pt; k++) {
                    angle = 2.0L*pi*(double)k/(double)pt;
                    x[k+2] = cos(angle) * center_radius;
                    y[k+2] = sin(angle) * center_radius;
                }
			}
            first_y = y1;
		}
		else {
        //====================================================
        // TRANSITION PHASE
        //====================================================
            radius = sqrtl((x1 - x[2])*(x1 - x[2]) + (y1 - y[2])*(y1 - y[2]));
            pheta = atan2l(y1 - y[2], x1 - x[2]);
            if (pheta < 0.0L) pheta180 = 360.0L + pheta*180.0L/pi;
            else pheta180 = pheta*180.0L/pi;

            x[1] = x[2] + cosl(pheta) * (radius + rincr);
            y[1] = y[2] + sinl(pheta) * (radius + rincr);

            //radius / Bohr radius
            n_spiral = radius/a0;
            n_quantum = sqrtl(radius/a0);
            n_theory = sqrtl(1.0L + (long double)(counter - 1) /(long double)orbital_counter);

            //test for spiral angle
            if (nlevel == 2) {
                if (pheta_old > pheta180) {
                    prt = 1;
                    n_angle = 120.0L;
                }
            }
            else {
                if (pheta180 >= n_angle && pheta_old < n_angle) {
                    prt = 1;
                    n = ns + 1.0L;
                    n_angle = 720.0L * (n*n - n)/(n*n) - 360.0L;
                }
            }
		}




		full_counter = (counter - 1) + orbital_counter;
		if (prt > 0){
                printf("angles = %.0lf   %.9lf   %.9lf \n", nlevel*1.0, (double)pheta180, (double)pheta_old);
                printf("nlevel = %.9lf   %.9lf   %.9lf   %.9lf    %.9lf \n", (double)n_spiral, (double)n_spiral_old, (double)n_quantum, (double)n_theory, (double)n_theory_old);
                printf("count.... %.0lf   %.9lf   %.0lf  %.9lf  ", counter*1.0,  (double) (counter/lambda),  full_counter*1.0,  (double)(full_counter/lambda) );
                if (nlevel <= 4) printf("H spectra... %.9lf  \n", (double)((full_counter/lambda) - fh[nlevel]) );
                else printf("\n");
                printf("length... %.9lf   %.9lf \n", (double)length, (double)(length/lorbital_n1) );
                printf("xy values... x1=%.9lf   y1=%.9lf     x2=%.9lf  y2=%.9lf   %.0lf \n", (double)x[1], (double)y[1], (double)x[2], (double)y[2], j*1.0);
                n = (n_spiral - 1.0L) * 4.0L * pi * c * lambda / ((le + lp) * (long double)full_counter);
                printf("analysis = %.11lf,  %.11lf,  %.11lf,  %.11lf,  %.1lf \n\n",(double)n_quantum, (double)n_spiral,  (double)(full_counter/lambda) , (double)pheta180,  (double)n);
                prt = 0;
                nlevel++;
                ns = nlevel * 1.0L;
                if (nlevel > max_nlevel) j = maxj - 1;
		}
        pheta_old = pheta180;
        n_spiral_old = n_spiral;
        n_theory_old = n_theory;



//
        if (counter > (prev_counter + save_interval) && n_quantum > 1.0L) {
                //
                    fprintf(data_out, "%.11f  %.11f  %.11f  %.11f  %.11f  %.11f  %.11f  %.11f  %.0f\n",
                        (double)n_quantum,
                        (double)(radius),
                        (double)(x1),
                        (double)(y1),
                        (double)(x[2]),
                        (double)(y[2]),
                        (double)(pheta180),
                        (double)(length),
                        (double)(full_counter)
                         );
                        //
                prev_counter = counter;
                Nsamples++;
                //fprintf(xy_out, "%.6f  %.6f\n", (double)(x1/r0), (double)(y1/r0));
        }
//


	}	//for i
	}	//for j

    dt_sim = (long double)counter / ( (long double)Nsamples * lambda );
	printf("samples ...  %.0lf  \n", full_counter*1.0);
    printf("file ...  %.0lf    %.0lf    %.11lf   %.11lf   \n", max_nlevel*1.0,  (double)jpoints,  (double)rplus , (double)lorbital_n1);
    fclose(xy_out);
    fclose(data_out);
    scanf("%d", &n);    //prevents window from closing
}