//g++ calc_msdlog_dis.cpp random_mars.cpp  -o calc_msdlog_dis -O2
//./calc_msdlog_dis 93757359 100 200 

// code snippets from :
//https://kops.uni-konstanz.de/entities/publication/9954ee47-48c8-4142-864e-f0cd6f93f9f0
//https://kops.uni-konstanz.de/entities/publication/22436d1a-9ff1-4779-aad3-94ea9fe6e625

#include <stdio.h> 
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string>
#include <sstream>
#include <complex>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include "random_mars.h"
#include <random> // Include this for the random number generation
#define PI 3.1415926535897

class RanMars *random_mars;
int seed;
// Global variables

double x_before, fx_before, fx_after;

int n_max, t_samp, resolution, min_rate;
double dt;
int64_t *counter;
double **msd;
double **isf;
double **isf_im;
double ****xo;
int ****imageo; //for periodic boundary conditions
int *length, *rate;
int **sample;

// für den langen code dann  int n_part = 1000; // Number of particles
int n_part = 1000; 
// Random number generator for Gaussian distribution
std::default_random_engine generator;
std::normal_distribution<double> distribution(0.0, 1.0);


void Init(int n_part, int n_step) {
    
    n_max = 0;
    t_samp = 0;
    length = NULL;
    rate = NULL;
    sample = NULL;
    counter = NULL;
    msd = NULL;
    xo = NULL;
    imageo = NULL;
    
    
    resolution = 8; // Number of points per decade
    min_rate = 10; // Min number of samples
    double a = pow(10.0, 1.0 / resolution);
    double e = log(min_rate) / log(n_step);
    int n_min = static_cast<int>(-log(a - 1) / log(a)) + 1;
    n_max = static_cast<int>(log(n_step) / log(a));
    n_max-=n_min;
    length = new int[n_max];
    rate = new int[n_max];
    sample = new int*[n_max];

    for (int n = 0; n < n_max; n++) {
        length[n] = static_cast<int>(pow(a, n + 1 + n_min) - pow(a, n_min));
        rate[n] = static_cast<int>(ceil(pow(length[n], e)));
        sample[n] = new int[rate[n]];
        for (int nn = 0; nn < rate[n]; nn++) {
            sample[n][nn] = length[n] * nn / rate[n];
        }
    }


    // Allocate memory
    counter = new int64_t[n_max];
    msd = new double*[n_max];
    isf = new double*[n_max];
    isf_im= new double*[n_max];
    xo = new double***[n_max];
    imageo = new int***[n_max];

    for (int n = 0; n < n_max; n++) {
        msd[n] = new double[4];
        isf[n] = new double[6];
        isf_im[n] = new double[6];
        for (int d = 0; d < 4; d++) {
            msd[n][d] = 0.0;
        }
        for (int d = 0; d < 6; d++) {
            isf[n][d] = 0.0;
            isf_im[n][d] = 0.0;
        }
        counter[n] = 0;
        xo[n] = new double**[rate[n]];
        imageo[n] = new int**[rate[n]];
        for (int nn = 0; nn < rate[n]; nn++) {
            xo[n][nn] = new double*[n_part];
            imageo[n][nn] = new int*[n_part];
            for (int i = 0; i < n_part; i++) {
                xo[n][nn][i] = new double[3];
                imageo[n][nn][i] = new int[3];
            }
        }
    }
    //Monte carlo
    x_before=0.0;
	fx_before=0.0;
	fx_after=0.0;
}



void Sample(int t_samp, double **x, int **image, const double *boxlength) {
    double q;
    double Qnu;
    double Qmu;
    const std::complex<double> j(0.0, 1.0);

    for (int n = 0; n < n_max; n++) {
        for (int nn = 0; nn < rate[n]; nn++) {
            if ((t_samp % length[n]) == sample[n][nn]) {
                if (t_samp > sample[n][nn]) {
                    // Calculate MSD if the time origin already exists
                    for (int i = 0; i < n_part; i++) {

                        double dx = x[i][0] - xo[n][nn][i][0];                 
                        msd[n][0] += dx ;
                        msd[n][1] += dx * dx;
                        msd[n][2] += dx * dx * dx;
                        msd[n][3] += dx * dx * dx * dx;

                        const double factors[] = {1*2*PI, 1*2*PI, 1.5*2*PI, 1.5*2*PI, 3.5*2*PI, 3.5*2*PI, 2.0*2*PI, 1.0*2*PI, 0.5*2*PI, 1.5*2*PI, 0.5*2*PI, -0.5*2*PI};
                        for (int j = 0; j < 6; j=j+1) {
                            isf[n][j]      += std::cos((-factors[2 * j] * x[i][0] + factors[2 * j+1] * xo[n][nn][i][0]));
                            isf_im[n][j]      += std::sin((-factors[2 * j] * x[i][0] + factors[2 * j+1] * xo[n][nn][i][0]));
                        }
                      
                    counter[n]++;
                    }
                    
                }
                // Set new time origin
                for (int i = 0; i < n_part; i++) {
                    for (int d = 0; d < 3; d++) {
                        xo[n][nn][i][d] = x[i][d];
                        imageo[n][nn][i][d] = image[i][d];
                    }
                }
            }
        }
    }
}

void Finish(const char *filename) {
    FILE *fp = fopen(filename, "w");
    if (fp == NULL) {
        perror("Error opening file");
        return;
    }

    for (int n = 0; n < n_max; n++) {
        for (int d = 0; d < 4; d++) {
            msd[n][d] /= counter[n];
        }
        for (int d = 0; d < 6; d++) {
            isf[n][d] /= counter[n];
            isf_im[n][d] /= counter[n];
        }
        double time = length[n] * dt;

    fprintf(fp, "%10.8f\t", time);
    for (int i = 0; i < 4; i++) {
        fprintf(fp, "%10.8f\t", msd[n][i]);
    }
    for (int i = 0; i < 6; i++) {
        fprintf(fp, "%10.8f\t", isf[n][i]);
    }
        for (int i = 0; i < 6; i++) {
        fprintf(fp, "%10.8f\t", isf_im[n][i]);
    }
    fprintf(fp, "%ld\t%d\n", counter[n], n);

      
    }

    fclose(fp);

    // Free allocated memory
    for (int n = 0; n < n_max; n++) {
        for (int nn = 0; nn < rate[n]; nn++) {
            for (int i = 0; i < n_part; i++) {
                delete[] imageo[n][nn][i];
                delete[] xo[n][nn][i];
            }
            delete[] imageo[n][nn];
            delete[] xo[n][nn];
        }
        delete[] xo[n];
        delete[] imageo[n];
        delete[] msd[n];
        delete[] isf[n];
        delete[] isf_im[n];
        delete[] sample[n];
    }
    delete[] xo;
    delete[] imageo;
    delete[] msd;
    delete[] isf;
    delete[] isf_im;
    delete[] counter;
    delete[] length;
    delete[] rate;
    delete[] sample;
}




int main(int argc, char *argv[]) {
    int seed = atoi(argv[1]);
    double u = std::atof(argv[2]);
    double f = std::atof(argv[3]);
	random_mars = new RanMars(seed);
    int n_step = 1000000; // Number of time steps //used  10000000
    dt = 0.0001; // Time step size used 0.001
	double B = sqrt (2.0*dt );
    // Initialize the algorithm
    Init(n_part, n_step);

    // Define x, image, and boxlength arrays (must be provided)
    double **fx = new double*[n_part];
    double **x = new double*[n_part];
    int **image = new int*[n_part];
    double boxlength[3] = {100.0, 100.0, 100.0}; // Example box length

    for (int i = 0; i < n_part; i++) {
        x[i] = new double[3];
        fx[i] = new double[3];
        image[i] = new int[3];
        // Initialize x and image to some initial condition
        for (int d = 0; d < 3; d++) {
            x[i][d] = 0.0;
            fx[i][d]=0.0;
            image[i][d] = 0;
        }
    }

    // Example sampling loop
    for (int t_samp = 0; t_samp < n_step; t_samp++) {
        // Simulate Brownian motion
        for (int i = 0; i < n_part; i++) {

                //double randi = distribution(generator); // Generate Gaussian random number
                
                double randi=random_mars->gaussian ();
                fx[i][0] = u*2.0*PI*sin(2*PI*x[i][0])+f*2.0*PI; 
                x[i][0] += fx[i][0]*dt+B* randi; // Update position with random displacement
                x[i][1] += B* randi; // Update position with random displacement
                //x[i][0] = u/(2*PI)*sin(2.0*PI*x[i][0]); 
                fx[i][0] = u*2.0*PI*sin(2*PI*x[i][0])+f*2.0*PI; 


                //Monte-Carlo version   
                /*
                double randi_x=random_mars->gaussian();
                x_before=x[i][0];
                fx_before=fx[i][0];
                x[i][0]+=B*randi_x;	
                x[i][1]+=B*randi_x;	
                //fx[i][0] = u*sin(2*PI*x[i][0]);
                fx[i][0] = u*sin(x[i][0]);
                
                fx_after=fx[i][0];

                double deltaE= fx_after-fx_before;
                if (deltaE> 0.0){
                    double randi_exp= random_mars->uniform();

                    if  (exp ( - deltaE) < randi_exp) {
                        x[i][0]=x_before;
                        fx[i][0]=fx_before;
            
                    }
                }*/
           
        }
         
                
        // Sample MSD at the current time step
        Sample(t_samp, x, image, boxlength);
    }

    // Finalize and output the results



    std::ostringstream fname;
    fname << "data_observables_variable_name/ISF/"
        << "msd_u" << std::fixed << std::setprecision(2) << u
        << "_f" << std::fixed << std::setprecision(2) << f << ".txt";

    Finish(fname.str().c_str());
    // Free x and image memory
    for (int i = 0; i < n_part; i++) {
        delete[] x[i];
        delete[] image[i];
        delete[] fx[i];
    }
    delete[] x;
    delete[] image;
    delete[] fx;
    return 0;
}