//compile with: g++ TPS_rtp.cpp nrutil.cpp random_mars.cpp -o TPS_rtp -O2 -w

#include <fstream>

#include <iomanip>

#include <iostream>

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

#include <sys/time.h>

#include <math.h>

#include <string>

#include <string.h>

#include <sstream>

#include "TPS_rtp.h"

#include "nrutil.h"

#include "random_mars.h"

#include <random>

#include <iostream>

#include <ctime>

#include <cstdlib>

#include <vector>

#include <chrono>

#define PI 3.1415926535897932384626433832795028841971693993751
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Main Function
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char * argv[]) {
  if (argc != 2) {
    std::cerr << "wrong, plese use as " << argv[0] << " <number_of_screen_u_use>" << std::endl;
    return 1;
  }
  int screennum = std::stoi(argv[1]); // Convert the first argument to double
  std::cout << "it worked, screen_num: " << screennum << std::endl;
  inftyrunrate = false; // default is false for RTP. if you want infinity passive to active rate, set to true
  dt = 0.005; // time-step length, 0.005 standard
  Nsim = 20000; // max number of timesteps for a trajectory
  indexshift = Nsim / 4;
  int sample_size = 10000000;

  // system parameters
  Tinit = 1.0; //temperature in units of kB
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);

  Pe = 10; //0.0 5.0; 10.0;	// Peclet number 
  mu = 0.1;
  D_th = 1.0; // standard 1.0; 
  D = 0.1; // standard 0.10;
  //Tinit * mu;
  B = sqrt(2.0 * D * dt);
  Bth = sqrt(2.0 * D_th * dt);
  v_drive = Pe / (sqrt(3 / (4 * D * D_th))); // drive velocity of active particle

  // potential parameters
  kx = 6.0;
  ky = 20.0;
  x0 = -1.0;

  // positions, trajectories
  x_one = 0.;
  y_one = 0.;
  theta_one = 0.;
  fx_one = 0.;
  fy_one = 0.;
  U_one = 0.;
  nu_one = 0.;
  R = dmatrix(0, 6, 0, Nsim - 1); //Trajectory = (x, y, theta, fx, fy, nu)
  R_new = dmatrix(0, 6, 0, Nsim - 1); //new Trajectory
  times = dmatrix(0, 3, 0, 21 * 21);
  histosize = static_cast < int > (10000);
  Histo = ivector(0, 10000);

  // steady state distribution
  Ngrid = 104; // number of grid cells per dimension
  xmin = -1.26;
  xmax = -0.65;
  deltx = (xmax - xmin) / Ngrid;
  ymin = -0.45;
  ymax = 0.45;
  delty = (ymax - ymin) / Ngrid;
  thetamin = 0.;
  thetamax = 2. * PI;
  deltatheta = (thetamax - thetamin) / Ngrid;
  index_x = 0;
  index_y = 0;
  index_theta = 0;
  Neq = 1000000000; // number of timesteps for steady state distribution histogram
  eq_rounds = 100; // if lambda10, lambda01 are in (0.1, 1, 10), then this could be increased to 100. most data i did with 50 tho
  grid_count_run = f3tensor(0, Ngrid, 0, Ngrid, 0, Ngrid); // steady state histogram
  grid_count_tumble = f3tensor(0, Ngrid, 0, Ngrid, 0, Ngrid); // steady state histogram

  grid_dx = 0.1;
  grid_dy = 0.1;
  grid_x_left = -3.0;
  grid_x_right = 3.0;
  grid_y_left = -3.0;
  grid_y_right = 3.0;

  n_grid_x_left = static_cast < int > (std::round(grid_x_left / grid_dx));
  n_grid_x_right = static_cast < int > (std::round(grid_x_right / grid_dx));
  n_grid_y_left = static_cast < int > (std::round(grid_y_left / grid_dy));
  n_grid_y_right = static_cast < int > (std::round(grid_y_right / grid_dy));

  int dim_x = n_grid_x_right - n_grid_x_left + 1;
  int dim_y = n_grid_y_right - n_grid_y_left + 1;
  n_reactiveJ = imatrix(0, dim_x - 1, 0, dim_x - 1);
  reactiveprobdensity = dmatrix(0, dim_x - 1, 0, dim_x - 1);

  reactiveJ = f3tensor(0, dim_x, 0, dim_y, 0, 1); //f3tensor

  for (int i = 0; i < dim_x; i++) {
    for (int j = 0; j < dim_y; j++) {
      reactiveprobdensity[i][j] = 0.;
      n_reactiveJ[i][j] = 0;
      reactiveJ[i][j][0] = 0;
      reactiveJ[i][j][1] = 0;
    }
  }
  seed = screennum;
  printf("seed=%i\n", seed);
  srand(seed);
  fflush(stdout);
  random_mars = new RanMars(seed);

  B = sqrt(2.0 * D * dt);
  Bth = sqrt(2.0 * D_th * dt);
  v_drive = Pe / (sqrt(3 / (4 * D * D_th))); // drive velocity of active particle
  histosize = static_cast < int > (10000);
  printf("%f\n", dt);
  for (int ni = 0; ni <= histosize; ni++) {
    Histo[ni] = 0;
  }
  // different values of dt for RTP with inftyrunrate approaches -> RTPE with instantaneous tumbling events (fig 6)
  //brute_force_lengthdist_rtpe(100000);
  //brute_force_lengthdist(100000);

  /////////////////////////////////////////////////////////////////////////////////////////////////////
  // Sampling sampling(Pe, tr, rr, sample_size, TPS, BF, SSD);
  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
  sample_size = 1000000; // standard precision for mean TPT
  
  //iterate for different tr and rr for trajectories (fig 1) pathlength distributions (figs 2, 3) as well as reactive probability density and current
  
  tr = 1; // 0.12 1/tau
  rr = 1; // 0.12 1/tau
  sampling(Pe, tr, rr, sample_size, true, false, true); // creates steady-state distribution first, then does some TPS for parameters Pe, tr and rr. samplesize 

  return 0;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Functions
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/*********************************************************************************************************************************/
void sampling(double peclet, double tumblerate, double runrate, int sampsize, bool TPS, bool BF, bool steadystate) {
  Pe = peclet;
  tr = tumblerate;
  rr = runrate;
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  if (inftyrunrate == true) {
    tumbleprob = 0;
  }

  //if(rr==0 && tr == 0)
  print_md_trajectory(); 
  print_TPS_trajectories();
  if (steadystate == true) {
    SteadyState_distribution();
  }
  if (BF == true) {
    brute_force_lengthdist(sampsize);
  }
  if (TPS == true) {
    TPS_lengthdist(sampsize); //
  }
}

/*********************************************************************************************************************************/

void print_md_trajectory() {
  std::stringstream sstr;
  sstr << "../data/trajectories/MDtrajectory_rtp_step" << dt << "pe" << Pe << "tumblerate" << tr << "runrate" << rr << ".dat"; //
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "wb");
  int counter = 0; // print brute forced trajectory with past
  lattice_packing();
  bool r = false;
  bool t = false;
  int i = 0;
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  do {
    U_one = calculate_potential(x_one, y_one);
    //path_length++;	
    if (U_one < 2.0 && x_one < 0.) {
      r = true;
    }
    if (U_one < 2.0 && x_one > 0.) {
      t = true;
    }
    fprintf(out1, "%i\t%f %f %f %f\n", i, x_one, y_one, theta_one, nu_one);
    md_step();
  }
  while (t == false);
  fclose(out1);
}
/*********************************************************************************************************************************/
void print_TPS_trajectories() { // marker
  lattice_packing();
  R_found = false;
  do {
    find_first_reactive_R(); // writing on R[][], the very first R old

  }
  while (R_found == false);
  if (R_found == true) {
    path_length = n_max - n_min;
  } else {
    printf("no initial R found, try different seed\n");
  }
  do {
    new_R();
    path_length_new = n_max_new - n_min_new;
    // printf("x_new_nminnew = %f\n", R_new[0][n_min_new]);
  }
  while (backward == false || forward == false); // for long trajectory
  //while (backward == false || forward == false || path_length_new >100);
  std::stringstream sstr;
  sstr << "../data/trajectories/tps_trajectories_rtp_step" << dt << "pe" << Pe << "tumblerate" << tr << "runrate" << rr << ".dat"; //
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "wb");
  for (int i = 1; i < Nsim; i++) // print old and new trajectory
  {
    double xloc = R[0][i];
    double yloc = R[1][i];
    double thetaloc = R[2][i];
    double nuloc = R[5][i];
    double uloc = calculate_potential(R[0][i], R[1][i]);

    double xlocn = R_new[0][i];
    double ylocn = R_new[1][i];
    double thetalocn = R_new[2][i];
    double nulocn = R_new[5][i];
    double ulocn = calculate_potential(R_new[0][i], R_new[1][i]);
    fprintf(out1, "%i\t%f %f %f %f %f \t%f %f %f %f %f\n", i, xloc, yloc, thetaloc, nuloc, uloc, xlocn, ylocn, thetalocn, nulocn, ulocn);
  }
  fclose(out1);
}

/*********************************************************************************************************************************/
void brute_force_lengthdist(int steps) {
  for (int ni = 0; ni <= histosize; ni++) {
    Histo[ni] = 0;

  }
  std::stringstream sstr;
  sstr << "../data/length_distributions/length_distribution_histo_md_rtp_step" << dt << "pe" << Pe << "tumblerate" << tr << "runrate" << rr << "samplesize" << steps << ".dat"; //
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "wb");

  printf("\n\nbrute force sampling\npe=%f\t tr=%f\t rr=%f\n", Pe, tr, rr);
  printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n");
  printf("progress:\n");
  struct timeval tp;
  gettimeofday( & tp, NULL);
  long int ms = tp.tv_sec * 1000 + tp.tv_usec / 1000;

  int counter = 0;
  lattice_packing();
  bool r = false;
  bool t = false;
  int i = 0;
  double avgruns = 0;
  double avgtumbles = 0;
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  if (inftyrunrate == true) {
    tumbleprob = 0;
  }
  do {
    U_one = calculate_potential(x_one, y_one);
    path_length++;
    if (U_one < 2.0 && x_one < 0.) {
      r = true;
      if (t == true) {
        //fprintf(out1, "%i\n", path_length);
        if (path_length <= histosize) {
          Histo[path_length]++;
        } else {
          Histo[histosize]++;
        }
        //				printf("trajectory found\n");
        counter++;
        t = false;
        i = counter;
        if ((100 * i) % (steps) == 0) {
          printf("%f percent\n", (double)(i * 100) / (steps));
        }
      }
      path_length = 0;
    }
    if (U_one < 2.0 && x_one > 0.) {
      t = true;
      if (r == true) {
        if (path_length <= histosize) {
          Histo[path_length]++;
        } else {
          Histo[histosize]++;
        }
        //				printf("trajectory found\n");
        counter++;
        r = false;
        i = counter;
        if ((100 * i) % (steps) == 0) {
          printf("%f percent\n", (double)(i * 100) / (steps));
        }
      }
      path_length = 0;
    }
    //printf("%i\t%f %f %f %f\n", i, x_one, y_one, theta_one, nu_one);
    md_step();
  }
  while (counter < steps);
  for (int lu = 0; lu <= histosize; lu++) //marker
  {
    fprintf(out1, "%i\n", Histo[lu]);
  }
  fclose(out1);
  struct timeval tpend;
  gettimeofday( & tpend, NULL);
  long int msend = tpend.tv_sec * 1000 + tpend.tv_usec / 1000;
  long int difftime = msend - ms;
  printf("Brute Force Sampling Time: %d ms \n", difftime);
  md_time = difftime;
  printf("done");

}

/*********************************************************************************************************************************/
void brute_force_lengthdist_rtpe(int steps) {

  for (int ni = 0; ni <= histosize; ni++) {
    Histo[ni] = 0;

  }

  std::stringstream sstr;
  sstr << "../data/length_distributions/length_distribution_histo_md_rtpE_step" << dt << "pe" << Pe << "tumblerate" << tr << "samplesize" << steps << ".dat"; //
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "wb");

  printf("\n\nbrute force sampling RTPE \npe=%f\t tr=%f\n", Pe, tr);
  printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n");
  printf("progress:\n");
  struct timeval tp;
  gettimeofday( & tp, NULL);
  long int ms = tp.tv_sec * 1000 + tp.tv_usec / 1000;

  int counter = 0;

  lattice_packing();

  bool r = false;
  bool t = false;
  int i = 0;
  double avgruns = 0;
  double avgtumbles = 0;
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  if (inftyrunrate == true) {
    tumbleprob = 0;

  }

  do {
    U_one = calculate_potential(x_one, y_one);
    path_length++;
    if (U_one < 2.0 && x_one < 0.) {
      r = true;
      if (t == true) {
        //fprintf(out1, "%i\n", path_length);
        if (path_length <= histosize) {

          Histo[path_length]++;
        } else {
          Histo[histosize]++;
        }
        //				printf("trajectory found\n");
        counter++;
        t = false;
        i = counter;
        if ((100 * i) % (steps) == 0) {
          printf("%f percent\n", (double)(i * 100) / (steps));
        }
      }
      path_length = 0;
    }
    if (U_one < 2.0 && x_one > 0.) {
      t = true;
      if (r == true) {
        if (path_length <= histosize) {
          Histo[path_length]++;
        } else {
          Histo[histosize]++;
        }
        //				printf("trajectory found\n");
        counter++;
        r = false;
        i = counter;
        if ((100 * i) % (steps) == 0) {
          printf("%f percent\n", (double)(i * 100) / (steps));
        }
      }
      path_length = 0;
    }
    //printf("%i\t%f %f %f %f\n", i, x_one, y_one, theta_one, nu_one);
    md_step_RTPE();
  }
  while (counter < steps);
  for (int lu = 0; lu <= histosize; lu++) //marker
  {
    fprintf(out1, "%i\n", Histo[lu]);
  }
  fclose(out1);
  struct timeval tpend;
  gettimeofday( & tpend, NULL);
  long int msend = tpend.tv_sec * 1000 + tpend.tv_usec / 1000;
  long int difftime = msend - ms;
  printf("Brute Force Sampling Time RTPE: %d ms \n", difftime);
  md_time = difftime;
  printf("done");

}
/*********************************************************************************************************************************/
void TPS_lengthdist(int steps) {
  int counting = 0;
  double acc_ratio = 0;
  double avgrun = 0;
  double avgtumble = 0;
  int prevruns = 0;
  int prevtumbles = 0;

  for (int ni = 0; ni <= histosize; ni++) {
    Histo[ni] = 0;
  }

  read_steady_state_dist();

  lattice_packing();
  acceptance_ratio = 0;
  // track TPS time
  struct timeval tp;
  gettimeofday( & tp, NULL);
  long int ms = tp.tv_sec * 1000 + tp.tv_usec / 1000;

  // output of path lengths to file
  std::stringstream sstr;

  sstr << "../data/length_distributions/length_distribution_histo_TPS_rtp_step" << dt << "pe" << Pe << "tumblerate" << tr << "runrate" << rr << "samplesize" << steps << ".dat";
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "wb");
  printf("\n\ntransition path sampling\npe=%f\t tr=%f\t rr=%f\n", Pe, tr, rr);
  printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n");
  printf("progress:\n");

  mean_length = 0;
  for (int i = 0; i < steps; i++) {
    if (i % (steps / 100) == 0) {
      //  printf("%f percent\n", (double)(i * 100) / (steps));
      printf("%i percent \n", i / (steps / 100));
    }

    new_R(); // 
    //counting++;
    //printf("old R starts at %i and ends at %i length %i and the random chosen step is %i \n", n_min, n_max, n_max-n_min, ind); 
    if (backward == false || forward == false) {
      avgrun += prevruns;
      avgtumble += prevtumbles;
      path_length = n_max - n_min;
      if (counting == 0) {
        //        printf("no trajectory generated, write to file old path length = %i\n",path_length);
        counting = 0;
        //  fprintf(out1, "%i\n", path_length);
        if (path_length <= histosize) {

          Histo[path_length]++;
        } else {
          Histo[histosize]++;
        }
        mean_length += path_length;
      }

    }
    if (backward == true && forward == true) {
      accept_or_reject();

      if (accept == true) {
        acc_ratio++;
        avgrun += NumOfRuns;
        avgtumble += NumOfTumbles;
        prevruns = NumOfRuns;
        prevtumbles = NumOfTumbles;
        NumOfTumbles = 0;
        NumOfRuns = 0;

        if (counting == 0) {
          // 			printf("new trajectory accepted, write to file new path length = %i\n",path_length_new);
          counting = 0;
          //  fprintf(out1, "%i\n", path_length_new);
          if (path_length <= histosize) {
            Histo[path_length]++;
          } else {
            Histo[histosize]++;
          }
          mean_length += path_length_new;
        }
        //printf("new R %i to file \n\n", path_length_new);

      }
      if (accept == false) {
        avgrun += prevruns;
        avgtumble += prevtumbles;
        if (counting == 0) {
          //		  			printf("new trajectory rejected, write to file old path length = %i\n",path_length);
          counting = 0;
          //  fprintf(out1, "%i\n", path_length);
          if (path_length <= histosize) {
            Histo[path_length]++;
          } else {
            Histo[histosize]++;
          }
          mean_length += path_length;

        }
        //printf("old R %i to file \n\n", path_length);
      }
    }
    if (i % (steps / 10) == 0) {
      calcReactiveFlux();
    }
  }
  compute_final_averages2d();
  print_reactiveprobabilitydensity2d();
  print_reactivecurrent2d();

  for (int lu = 0; lu <= histosize; lu++) //marker
  {
    fprintf(out1, "%i\n", Histo[lu]);
  }
  mean_length = mean_length / steps;
  acceptance_ratio = acc_ratio / steps;
  fclose(out1);
  struct timeval tpend;
  gettimeofday( & tpend, NULL);
  long int msend = tpend.tv_sec * 1000 + tpend.tv_usec / 1000;
  long int difftime = msend - ms;
  printf("TPS Sampling Time: %d ms \n", difftime);
  tps_time = difftime;
  printf("done");
}
/*********************************************************************************************************************************/
void accept_or_reject() {
  v_drive = Pe / (sqrt(3 / (4 * D * D_th)));
  path_length = n_max - n_min;
  path_length_new = n_max_new - n_min_new;
  double xi = R[0][n_min]; // position 0 on old R, 
  double yi = R[1][n_min];
  double thetai = R[2][n_min];
  double nui = R[5][n_min];
  index_x = (xi - xmin) / deltx; // dropping the remainder with (int)
  index_y = (yi - ymin) / delty;
  index_theta = (thetai - thetamin) / deltatheta;
  double rho0_old = 0.;

  if (nui == 1.) {
    rho0_old = grid_count_run[index_x][index_y][index_theta];
  }
  if (nui == 0.) {

    rho0_old += grid_count_tumble[index_x][index_y][index_theta];

  }
  xi = R_new[0][n_min_new]; // position 0 on new R
  yi = R_new[1][n_min_new];
  thetai = R_new[2][n_min_new];
  nui = R_new[5][n_min_new];
  index_x = (xi - xmin) / deltx; // dropping the remainder with (int)
  index_y = (yi - ymin) / delty;
  index_theta = (thetai - thetamin) / deltatheta;
  double rho0_new = 0.;
  if (nui == 1) {
    rho0_new = grid_count_run[index_x][index_y][index_theta]; //
  }
  if (nui == 0) {
    rho0_new = grid_count_tumble[index_x][index_y][index_theta]; // 
  }

  // boltzmann distribution of initial state
  double pi_old = exp(kx * (((R[0][n_min] * R[0][n_min]) - (x0 * x0)) * (((R[0][n_min] * R[0][n_min])) - (x0 * x0))) + (ky / 2. * (R[1][n_min] * R[1][n_min])));
  double pi_new = exp(kx * (((R_new[0][n_min_new] * R_new[0][n_min_new]) - (x0 * x0)) * (((R_new[0][n_min_new] * R_new[0][n_min_new])) - (x0 * x0))) + (ky / 2. * (R_new[1][n_min_new] * R_new[1][n_min_new])));

  //probability of dynamics
  double p1 = 0.; // initialization
  for (int i = n_min; i < ind; i++) {
    if (R[5][i] == 1.) {
      p1 += // backward
        ((-v_drive) / (2 * D) *
          (
            ((cos(R[2][i])) * ( // like luigi: cos(R[2][i + 1])
              R[0][i] - R[0][i + 1] -
              mu * dt * R[3][i + 1])) +
            ((sin(R[2][i])) * ( // like luigi: sin(R[2][i + 1])
              R[1][i] - R[1][i + 1] -
              mu * dt * R[4][i + 1]))
          )
        );
      p1 -= (
        //forward
        (-v_drive) / (2 * D) *
        (
          (cos(R[2][i]) * (
            R[0][i] - R[0][i + 1] +
            mu * dt * R[3][i])) +
          (sin(R[2][i]) * (
            R[1][i] - R[1][i + 1] +
            mu * dt * R[4][i]))
        )
      );
    }
  }
  for (int j = n_min_new; j < ind; j++) // new R
  {
    if (R_new[5][j] == 1.) {
      p1 +=
        (
          (-v_drive) / (2 * D) //forward
          *
          (
            cos(R_new[2][j]) * (
              R_new[0][j] - R_new[0][j + 1] +
              mu * dt * R_new[3][j]) +
            sin(R_new[2][j]) * (
              R_new[1][j] - R_new[1][j + 1] +
              mu * dt * R_new[4][j])
          )

        );
      p1 -=
        (
          //backward
          (-v_drive) / (2 * D) *
          (
            ((cos(R_new[2][j])) * ( // like luigi : cos(R_new[2][j + 1])
              R_new[0][j] - R_new[0][j + 1] -
              mu * dt * R_new[3][j + 1])) +
            ((sin(R_new[2][j])) * ( // like luigi : sin(R_new[2][j + 1])
              R_new[1][j] - R_new[1][j + 1] -
              mu * dt * R_new[4][j + 1]))
          )
        );
    }
  }
  p1 = exp(p1);

  double probability = 0;
  double condition = 0;

  probability = (double)(path_length) / ((double) path_length_new) * rho0_new / rho0_old * pi_old / pi_new * p1;

  probability *= transition_prob(R_new[5][n_min_new], R[5][n_min]);
  probability /= transition_prob(R[5][n_min], R_new[5][n_min_new]);

  condition = rand() / ((double) RAND_MAX);
  if (probability >= condition) {
    accept = true;
    if (Nsim - n_max > indexshift && n_min > indexshift) // make sure to not reach the edges of the trajectory array
    {
      for (int i = 0; i < Nsim; i++) {
        R[0][i] = R_new[0][i];
        R[1][i] = R_new[1][i];
        R[2][i] = R_new[2][i];
        R[3][i] = R_new[3][i];
        R[4][i] = R_new[4][i];
        R[5][i] = R_new[5][i];
        if (n_min_new <= i && i <= n_max_new) {
          //printf("%f\n",R_new[5][i]);  	
        }
      }
      //printf("\n");
      n_max = n_max_new;
      n_min = n_min_new;

    }

    if (Nsim - n_max <= indexshift) // shift to lower array index if trajectory comes too near to the max index
    {
      for (int i = 0; i < Nsim; i++) {
        R[0][i] = 0;
        R[1][i] = 0;
        R[2][i] = 0;
        R[3][i] = 0;
        R[4][i] = 0;
        R[5][i] = 0;
      }
      for (int i = n_min_new; i < n_max_new; i++) {
        R[0][i - indexshift] = R_new[0][i];
        R[1][i - indexshift] = R_new[1][i];
        R[2][i - indexshift] = R_new[2][i];
        R[3][i - indexshift] = R_new[3][i];
        R[4][i - indexshift] = R_new[4][i];
        R[5][i - indexshift] = R_new[5][i];
      }
      n_max = n_max_new - indexshift;
      n_min = n_min_new - indexshift;

    }

    if (n_min <= indexshift) // shift to higher array index if trajectory comes too near to the min index
    {
      for (int i = 0; i < Nsim; i++) {
        R[0][i] = 0;
        R[1][i] = 0;
        R[2][i] = 0;
        R[3][i] = 0;
        R[4][i] = 0;
        R[5][i] = 0;
      }
      for (int i = n_min_new; i < n_max_new; i++) {
        R[0][i + indexshift] = R_new[0][i];
        R[1][i + indexshift] = R_new[1][i];
        R[2][i + indexshift] = R_new[2][i];
        R[3][i + indexshift] = R_new[3][i];
        R[4][i + indexshift] = R_new[4][i];
        R[5][i + indexshift] = R_new[5][i];
      }
      n_max = n_max_new + indexshift;
      n_min = n_min_new + indexshift;

    }

  } else {
    accept = false;
    //printf("path rejected\n length = %i; probability %f; condition %f\n",path_length_new, probability, condition);
    //printf("(path_length)/(path_length_new)*rho0_new/rho0_old*pi_old/pi_new * exp (p1 + p2 + p3 + p4) = \n = %i/%i*%f/%f*%f/%f * epx( %f + %f + %f + f)\n\n",path_length,path_length_new,rho0_new,rho0_old,pi_old,pi_new, p1 , p2 , p3 , p4);
  }
}

/*********************************************************************************************************************************/
void read_steady_state_dist() // 
{
  float myvariable = 0;
  std::stringstream sstr;
  sstr << "../data/rtp_steadystate/SteadyStateDistribution_rtp_step" << dt << "pe" << Pe << "rounds" << eq_rounds << "tumblerate" << tr << "runrate" << rr << ".dat";
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "r");
  for (int yi = 0; yi < Ngrid; yi = yi + 1) {
    for (int xi = 0; xi < Ngrid; xi = xi + 1) {
      for (int thetai = 0; thetai < Ngrid; thetai = thetai + 1) {
        fscanf(out1, "%f", & myvariable);
        grid_count_run[xi][yi][thetai] = myvariable;
        if (myvariable < 5.0) {}
        fscanf(out1, "%f", & myvariable);
        grid_count_tumble[xi][yi][thetai] = myvariable;
      }
    }
  }
  fclose(out1);
}
/*********************************************************************************************************************************/
void SteadyState_distribution() {
  // min and max values chosen such that xmin is a bit smaller and xmax a bit larger, such that the grid is a bit larger then the reactant/target.
  printf("tr=%f rr=%f pe=%f\n", tr, rr, Pe);
  for (int yi = 0; yi < Ngrid; yi++) {
    for (int xi = 0; xi < Ngrid; xi++) {
      for (int thetai = 0; thetai < Ngrid; thetai++) {
        grid_count_run[xi][yi][thetai] = 0;
        grid_count_tumble[xi][yi][thetai] = 0;
      }
    }
  }

  // Constructing the filename
  std::stringstream sstr;
  sstr << "../data/rtp_steadystate/SteadyStateDistribution_rtp_step" << dt << "pe" << Pe << "rounds" << eq_rounds << "tumblerate" << tr << "runrate" << rr << ".dat";
  const std::string filename = sstr.str();

  // Check if file exists
  FILE * existingFile = fopen(filename.c_str(), "r");
  if (existingFile != NULL) {
    fclose(existingFile);
    std::cout << "File " << filename << " already exists. Skipping...\n";
    return; // Skip file creation
  }

  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  for (int j = 0; j < eq_rounds; j++) {
    printf("cycle %i of %i, particle at position (x,y) = (%f, %f) \n", j, eq_rounds, x_one, y_one);
    for (int i = 0; i <= Neq; i++) {
      if (i % (Neq / 100) == 0) {
        printf("progress of ssd: step %i of %i in round %i of %i particle at position (x,y,theta) = (%f, %f, %f) \n", i, Neq, j + 1, eq_rounds, x_one, y_one, theta_one);
      }
      md_step();
      if (xmin <= x_one && xmax >= x_one && ymax >= y_one && ymin <= y_one) {
        index_x = (x_one - xmin) / deltx;
        index_y = (y_one - ymin) / delty;
        index_theta = (theta_one - thetamin) / deltatheta;
        if (nu_one == 1.) {
          grid_count_run[index_x][index_y][index_theta] = grid_count_run[index_x][index_y][index_theta] + 1;
        }
        if (nu_one == 0.) {
          //			index_theta = 1;
          grid_count_tumble[index_x][index_y][index_theta] = grid_count_tumble[index_x][index_y][index_theta] + 1;
        }
        //printf("x=%f xindex=%i y=%f yindex=%i gridcount=%f\n", x_one, index_x, y_one, index_y,grid_count_run[index_x][index_y][index_theta]);
      }
      if (-xmax <= x_one && -xmin >= x_one && ymax >= y_one && ymin <= y_one) {
        x_one = -x_one;
        index_x = (x_one - xmin) / deltx; // dropping the remainder with (int)
        index_y = (y_one - ymin) / delty;
        theta_one = PI - theta_one;
        if (theta_one < 0.) {
          theta_one = theta_one + 2 * PI;
        }
        index_theta = (theta_one - thetamin) / deltatheta;
        if (nu_one == 1.) {
          grid_count_run[index_x][index_y][index_theta] = grid_count_run[index_x][index_y][index_theta] + 1;
        }
        if (nu_one == 0.) {
          //			index_theta = 1;
          grid_count_tumble[index_x][index_y][index_theta] = grid_count_tumble[index_x][index_y][index_theta] + 1;
        }
      }
    }
  }

  // Open file for writing
  const char * cstr = filename.c_str(); // Use filename instead of tmp
  FILE * out1 = fopen(cstr, "wb");
  //fprintf(out1, "x range \ty range\ttheta range\t count\n");
  for (int yi = 0; yi < Ngrid; yi++) {
    for (int xi = 0; xi < Ngrid; xi++) {
      for (int thetai = 0; thetai < Ngrid; thetai++) {
        fprintf(out1, "%f\n", grid_count_run[xi][yi][thetai]);
        grid_count_run[xi][yi][thetai] = 0;
        fprintf(out1, "%f\n", grid_count_tumble[xi][yi][thetai]);
        grid_count_tumble[xi][yi][thetai] = 0;
      }
    }
  }
  fclose(out1);
}
/*********************************************************************************************************************************/
//2d case of run and tumble  motion in a potential for brute force simulations
//
void md_step() //
{

  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  if (inftyrunrate == true) {
    tumbleprob = 0;
  }
  runcond = (rand() / ((double) RAND_MAX));
  tumblecond = (rand() / ((double) RAND_MAX));

  v_drive = Pe / (sqrt(3 / (4 * D * D_th)));
  // integrate with Euler-Maruyama scheme RTP
  x_one = x_one + mu * fx_one * dt + nu_one * v_drive * cos(theta_one) * dt + B * random_mars -> gaussian();
  y_one = y_one + mu * fy_one * dt + nu_one * v_drive * sin(theta_one) * dt + B * random_mars -> gaussian();

  if (nu_one == 1.) {

    if (runprob < runcond) { // 1 -> 0
      nu_one = 0.;
      theta_one = (theta_one + Bth * random_mars -> gaussian());
      theta_one = theta_one - 2. * PI * floor(theta_one / (2. * PI));
    } else { // 1 -> 1
      nu_one = 1.;
      theta_one = (theta_one + Bth * random_mars -> gaussian());
      theta_one = theta_one - 2. * PI * floor(theta_one / (2. * PI));
    }
  } else if (nu_one == 0.) {
    if (tumbleprob < tumblecond) {
      nu_one = 1.; // 0 -> 1
      theta_one = (rand() / ((double) RAND_MAX)) * 2 * PI;
    } else {
      nu_one = 0.; // 0 -> 0
      theta_one = (theta_one + Bth * random_mars -> gaussian());
      theta_one = theta_one - 2. * PI * floor(theta_one / (2. * PI));
    }
  }

  //reculculate forces, now that particles have moved
  calculate_forces(x_one, y_one);
}
/*********************************************************************************************************************************/
//2d case of run and tumble  motion in a potential for brute force simulations
//
void md_step_RTPE() // run and tumble with tumble Event
{
  runprob = exp(-1. * tr * dt);
  tumblecond = (rand() / ((double) RAND_MAX));
  v_drive = Pe / (sqrt(3 / (4 * D * D_th)));
  // integrate with Euler-Maruyama scheme RTP
  x_one = x_one + mu * fx_one * dt + nu_one * v_drive * cos(theta_one) * dt + B * random_mars -> gaussian();

  y_one = y_one + mu * fy_one * dt + nu_one * v_drive * sin(theta_one) * dt + B * random_mars -> gaussian();
  theta_one = (theta_one + Bth * random_mars -> gaussian());
  theta_one = theta_one - 2. * PI * floor(theta_one / (2. * PI));

  if (runprob < tumblecond) {
    nu_one = 1.;
    theta_one = (rand() / ((double) RAND_MAX)) * 2 * PI;
  }
  //reculculate forces, now that particles have moved
  calculate_forces(x_one, y_one);
}

/*********************************************************************************************************************************/
void md_step_forward(int i) //for TPS or safing trajectory
{
  runcond = (rand() / ((double) RAND_MAX));
  tumblecond = (rand() / ((double) RAND_MAX));
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  v_drive = Pe / (sqrt(3 / (4 * D * D_th)));
  if (R_new[5][i] == 1.) {
    R_new[0][i + 1] = R_new[0][i] + mu * R_new[3][i] * dt + v_drive * cos(R_new[2][i]) * dt + B * random_mars -> gaussian();
    R_new[1][i + 1] = R_new[1][i] + mu * R_new[4][i] * dt + v_drive * sin(R_new[2][i]) * dt + B * random_mars -> gaussian();
    R_new[2][i + 1] = (R_new[2][i] + Bth * random_mars -> gaussian());
    R_new[2][i + 1] = R_new[2][i + 1] - 2. * PI * floor(R_new[2][i + 1] / (2. * PI));
    if (runprob < runcond) { // 1 -> 0
      R_new[5][i + 1] = 0.;
    } else { // 1 -> 1
      R_new[5][i + 1] = 1.;
    }
  } else if (R_new[5][i] == 0.) {
    R_new[0][i + 1] = R_new[0][i] + mu * R_new[3][i] * dt + B * random_mars -> gaussian();
    R_new[1][i + 1] = R_new[1][i] + mu * R_new[4][i] * dt + B * random_mars -> gaussian();
    R_new[2][i + 1] = (R_new[2][i] + Bth * random_mars -> gaussian());
    R_new[2][i + 1] = R_new[2][i + 1] - 2. * PI * floor(R_new[2][i + 1] / (2. * PI));

    if (tumbleprob < tumblecond) {
      R_new[5][i + 1] = 1.; // 0 -> 1
      R_new[2][i + 1] = (rand() / ((double) RAND_MAX)) * 2 * PI;
    } else {
      R_new[5][i + 1] = 0.; // 0 -> 0
    }
  }
  //reculculate forces, now that particles have moved
  R_new[3][i + 1] = (-4) * kx * R_new[0][i + 1] * (R_new[0][i + 1] * R_new[0][i + 1] - x0 * x0);
  R_new[4][i + 1] = (-1) * ky * R_new[1][i + 1];
}

/*********************************************************************************************************************************/
//2d case of run and tumble  motion in a potential backward in time
//
void md_step_back(int i) // for TPS
{
  v_drive = Pe / (sqrt(3 / (4 * D * D_th)));
  runcond = (rand() / ((double) RAND_MAX));
  tumblecond = (rand() / ((double) RAND_MAX));
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);
  R_new[2][i] = R_new[2][i + 1] + Bth * random_mars -> gaussian();
  R_new[2][i] = R_new[2][i] - 2. * PI * floor(R_new[2][i] / (2. * PI));
  if (R_new[5][i + 1] == 1.) {
    if (runprob < runcond) {
      R_new[5][i] = 0.; // 1 -> 0
      R_new[2][i] = (rand() / ((double) RAND_MAX)) * 2 * PI; //TODO
      R_new[0][i] = R_new[0][i + 1] + mu * R_new[3][i + 1] * dt + B * random_mars -> gaussian(); // like luigi: cos( R_new[2][i+1])
      R_new[1][i] = R_new[1][i + 1] + mu * R_new[4][i + 1] * dt + B * random_mars -> gaussian(); // like luigi: sin( R_new[2][i+1])
    } else { // 1 -> 1
      R_new[5][i] = 1.;
      R_new[0][i] = R_new[0][i + 1] + mu * R_new[3][i + 1] * dt - R_new[5][i] * v_drive * cos(R_new[2][i]) * dt + B * random_mars -> gaussian(); // like luigi: cos( R_new[2][i+1])
      R_new[1][i] = R_new[1][i + 1] + mu * R_new[4][i + 1] * dt - R_new[5][i] * v_drive * sin(R_new[2][i]) * dt + B * random_mars -> gaussian(); // like luigi: sin( R_new[2][i+1])
    }
  } else if (R_new[5][i + 1] == 0.) {
    if (tumbleprob < tumblecond) {
      R_new[5][i] = 1.; // 0 -> 1
      R_new[0][i] = R_new[0][i + 1] + mu * R_new[3][i + 1] * dt - R_new[5][i] * v_drive * cos(R_new[2][i]) * dt + B * random_mars -> gaussian(); // like luigi: cos( R_new[2][i+1])
      R_new[1][i] = R_new[1][i + 1] + mu * R_new[4][i + 1] * dt - R_new[5][i] * v_drive * sin(R_new[2][i]) * dt + B * random_mars -> gaussian(); // like luigi: sin( R_new[2][i+1])
    } else {
      R_new[5][i] = 0.; // 0 -> 0
      R_new[0][i] = R_new[0][i + 1] + mu * R_new[3][i + 1] * dt + B * random_mars -> gaussian(); // like luigi: cos( R_new[2][i+1])
      R_new[1][i] = R_new[1][i + 1] + mu * R_new[4][i + 1] * dt + B * random_mars -> gaussian(); // like luigi: sin( R_new[2][i+1])
    }
  }
  //reculculate forces, now that particles have moved
  R_new[3][i] = (-4) * kx * R_new[0][i] * (R_new[0][i] * R_new[0][i] - x0 * x0);
  R_new[4][i] = (-1) * ky * R_new[1][i];
}
/*********************************************************************************************************************************/
// force due to potential
void calculate_forces(double x, double y) {
  fx_one = (-4) * kx * x * (x * x - x0 * x0);
  fy_one = (-1) * ky * y;
}

/*********************************************************************************************************************************/

void new_R() {
  runprob = exp(-1. * tr * dt);
  tumbleprob = exp(-1. * rr * dt);

  double U = 0.;
  for (int i = 0; i < Nsim; i++) {
    R_new[0][i] = 0; // x new
    R_new[1][i] = 0; // y new
    R_new[2][i] = 0; // theta new
    R_new[3][i] = 0; // fx new
    R_new[4][i] = 0; // fy new
    R_new[5][i] = 0; // nu new
  }
  forward = false;
  backward = false;
  // chose a point on R at random
  ind = (rand() / ((double) RAND_MAX)) * (n_max - n_min) + n_min;

  //  printf("old R starts at %i and ends at %i length is %i and the random chosen step is %i \n", n_min, n_max,n_max- n_min, ind);
  R_new[0][ind] = R[0][ind];
  R_new[1][ind] = R[1][ind];
  R_new[2][ind] = R[2][ind];
  R_new[3][ind] = R[3][ind];
  R_new[4][ind] = R[4][ind];
  R_new[5][ind] = R[5][ind];

  //forward_shooting
  for (int i = ind; i < Nsim; i++) // forward shooting
  {
    U_one = calculate_potential(R_new[0][i], R_new[1][i]);
    // target hit by forward shoot?
    if (U_one < 2.0 && R_new[0][i] > 0.) {
      //      printf("forward shooting hit target area xpos=%f , nmax = %i\n", R_new[0][i], i);
      forward = true;
      n_max_new = i; // 
      break;
    }
    if (U_one < 2.0 && R_new[0][i] < 0.) {
      //      printf("forward shooting hit reactant area xpos=%f nmax = %i\n", R_new[0][i],i);
      forward = false;
      break;
    }
    md_step_forward(i);
  }
  //backward_shooting
  for (int i = ind; i > 0; i--) {
    U_one = calculate_potential(R_new[0][i], R_new[1][i]); // reactant hit by backward shoot?
    if (U_one < 2.0 && R_new[0][i] < 0.) {
      //      printf("backward shooting hit reactant area xpos=%f, n_min=%i U = %f\n", R_new[0][i], i,U_one);
      backward = true;
      n_min_new = i;
      break;
    }
    if (U_one < 2.0 && R_new[0][i] > 0.) {
      //     printf("backward shooting hit target area xpos=%f, n_min=%i\n", R_new[0][i], i);
      backward = false;
      break;
    }
    md_step_back(i - 1);
  }
}

/*********************************************************************************************************************************/
void find_first_reactive_R() {
  lattice_packing();
  for (int i = 0; i < Nsim; i++) {
    R[0][i] = 0.;
    R[1][i] = 0.;
    R[2][i] = 0.;
    R[3][i] = 0.;
    R[4][i] = 0.;
    R[5][i] = 0.;
  }
  n_min = 0;
  R_found = false;
  U_one = calculate_potential(x_one, y_one);
  calculate_forces(x_one, y_one);

  for (int i = indexshift; i < Nsim - indexshift; i++) {
    md_step();
    U_one = calculate_potential(x_one, y_one);
    if (x_one < 0.0 && U_one < 2.0) {
      i = indexshift;
      n_min = i;
    }
    R[0][i] = x_one;
    R[1][i] = y_one;
    R[2][i] = theta_one;
    R[3][i] = fx_one;
    R[4][i] = fy_one;
    R[5][i] = nu_one;

    if (x_one > 0.0 && U_one < 2.0 && R_found == false) // condition for reaching the target: the potential is smaller then a given value, and x nees to be positive (else we would just find the start)
    {
      n_max = i;
      R_found = true;
      break;
    }
  }

  for (int i = n_max + 1; i < Nsim; i++) {
    R[0][i] = 0.;
    R[1][i] = 0.;
    R[2][i] = 0.;
    R[3][i] = 0.;
    R[4][i] = 0.;
    R[5][i] = 0.;

  }
  for (int i = 0; i < n_min; i++) {
    R[0][i] = 0.;
    R[1][i] = 0.;
    R[2][i] = 0.;
    R[3][i] = 0.;
    R[4][i] = 0.;
    R[5][i] = 0.;

  }
}

/*********************************************************************************************************************************/
// initialize in reactant area and only in 2d
void lattice_packing() {
  nu_one = floor(rand() / ((double) RAND_MAX) * 2);

  //printf("nu_one = %f\n", nu_one);
  x_one = x0; //+ 2*(0.5 - rand()/((double) RAND_MAX))*0.2;
  y_one = 0.; // + 2*(0.5 double transition- rand()/((double) RAND_MAX))*0.28955;
  theta_one = PI; //rand()/((double) RAND_MAX)*2.0*PI;
  calculate_forces(x_one, y_one);
  for (int j = 0; j < 80000; j++) {

    if (inftyrunrate == true) {
      md_step_RTPE();

    } else {
      md_step();

    }

  }

  do {
    if (inftyrunrate == true) {

      md_step_RTPE();
    } else {
      md_step();
    }
    calculate_potential(x_one, y_one);
  }
  while (U_one > 2.0 || x_one > 0.);
  // printf("Initialized... (x,y,theta, U) = (%f, %f, %f, %f, %f)\n", x_one, y_one, theta_one, U_one, nu_one);
}
/*********************************************************************************************************************************/
//calculate the potential energy
double calculate_potential(double x, double y) {
  U_one = kx * (((x * x) - (x0 * x0)) * (((x * x)) - (x0 * x0))) + (ky / 2 * (y * y));
  return U_one;
}
/*********************************************************************************************************************************/
double transition_prob(double nu_i, double nu_ii) {
  // tr = 10, rr = 0.1, dt = 0.01,
  double probab = 0.;
  if (nu_i == 1. && nu_ii == 1.) // run to run
  {
    //  probab = 1-exp(-tr * dt);
    probab = exp(-tr * dt);
  }
  if (nu_i == 1. && nu_ii == 0.) // run to tumble
  {
    //  probab = (exp(-tr * dt));
    probab = 1 - (exp(-tr * dt));
  }
  if (nu_i == 0. && nu_ii == 1.) // tumble to run
  {
    //  probab = exp(-rr * dt);
    probab = 1 - exp(-rr * dt);
  }
  if (nu_i == 0. && nu_ii == 0.) // tumble to tumble
  {
    //  probab = 1. - exp(-rr * dt);
    probab = exp(-rr * dt);
  }

  return (probab);
}

/*********************************************************************************************************************************/
void print_trajectories(int step, char folder[256]) // somehow causes segfault if done as function
{
  std::stringstream sstr;
  sstr << "../data/" << folder << "/trajectories_rtp" << step << ".dat";
  const std::string tmp = sstr.str();
  const char * cstr = tmp.c_str();
  FILE * out1 = fopen(cstr, "wb");
  for (int i = 1; i < Nsim; i++) {
    double xloc = R[0][i];
    double yloc = R[1][i];
    double thetaloc = R[2][i];
    double xlocn = R_new[0][i];
    double ylocn = R_new[1][i];
    double thetalocn = R_new[2][i];
    fprintf(out1, "%i\t%f %f %f %f\t%f %f %f %f\n", i, xloc, yloc, thetaloc, xlocn, ylocn, thetalocn);
  }
}

/*********************************************************************************************************************************/
void calcReactiveFlux() {
  //TODO deubug
  double vx, vy;
  int i, j;

  for (int k = n_min; k < n_max; k++) {
    vx = (R[0][k + 1] - R[0][k]) / dt;
    vy = (R[1][k + 1] - R[1][k]) / dt;

    i = static_cast < int > (std::round(R[0][k] / grid_dx));
    j = static_cast < int > (std::round(R[1][k] / grid_dy));

    if (i >= n_grid_x_left && i <= n_grid_x_right && j >= n_grid_y_left && j <= n_grid_y_right) {
      n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left]++;
      reactiveprobdensity[i - n_grid_x_left][j - n_grid_y_left] += 1.0;
      reactiveJ[i - n_grid_x_left][j - n_grid_y_left][0] += vx;
      reactiveJ[i - n_grid_x_left][j - n_grid_y_left][1] += vy;
    }
  }
}
/*********************************************************************************************************************************/
void compute_final_averages2d() {
  int somma = 0;

  for (int i = n_grid_x_left; i <= n_grid_x_right; ++i) {
    for (int j = n_grid_y_left; j <= n_grid_y_right; ++j) {
      somma += n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left];

    }
  }

  for (int i = n_grid_x_left; i <= n_grid_x_right; ++i) {
    for (int j = n_grid_y_left; j <= n_grid_y_right; ++j) {
      if (n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left] != 0) {
        reactiveJ[i - n_grid_x_left][j - n_grid_y_left][0] /= n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left];
        reactiveJ[i - n_grid_x_left][j - n_grid_y_left][1] /= n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left];
      }
      reactiveprobdensity[i - n_grid_x_left][j - n_grid_y_left] /= somma;

    }
  }
}
/*********************************************************************************************************************************/

void print_reactiveprobabilitydensity2d() {
  std::ostringstream filename;
  filename << "../data/ReactiveProbabilityDensity_RTP_Pe" << Pe <<
    "_RR" << rr <<
    "_TR" << tr <<
    "_timestep" << dt <<
    "_seed" << seed <<
    "_newversion_.dat";

  // Öffne die Datei mit dem dynamisch erstellten Namen
  std::ofstream outFile(filename.str());
  if (!outFile) {
    std::cerr << "Error opening file for writing!" << std::endl;
    return;
  }

  outFile << "# " << grid_x_left << " " << grid_x_right << " " << grid_dx << std::endl;
  outFile << "# " << grid_y_left << " " << grid_y_right << " " << grid_dy << std::endl;

  for (int i = n_grid_x_left; i <= n_grid_x_right; ++i) {
    for (int j = n_grid_y_left; j <= n_grid_y_right; ++j) {
      outFile << std::fixed << std::setprecision(3) << (i * grid_dx) << " " <<
        (j * grid_dy) << " " <<
        std::fixed << std::setprecision(5) << reactiveprobdensity[i - n_grid_x_left][j - n_grid_y_left] << std::endl;
    }
  }

  outFile.close();
}
/*********************************************************************************************************************************/
void print_reactivecurrent2d() {
  std::ostringstream filename;
  filename << "../data/ReactiveCurrent_RTP_Pe" << Pe <<
    "_RR" << rr <<
    "_TR" << tr <<
    "_timestep" << dt <<
    "_seed" << seed <<
    "_newversion.dat";
  std::ofstream outFile(filename.str());
  if (!outFile) {
    std::cerr << "Error opening file for writing!" << std::endl;
    return;
  }

  outFile << "# " << grid_x_left << " " << grid_x_right << " " << grid_dx << std::endl;
  outFile << "# " << grid_y_left << " " << grid_y_right << " " << grid_dy << std::endl;

  for (int i = n_grid_x_left; i <= n_grid_x_right; ++i) {
    for (int j = n_grid_y_left; j <= n_grid_y_right; ++j) {
      if (n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left] > 0) {
        outFile << std::fixed << std::setprecision(3) << (i * grid_dx) << " " <<
          (j * grid_dy) << " " <<
          std::fixed << std::setprecision(4) <<
          reactiveJ[i - n_grid_x_left][j - n_grid_y_left][0] << " " <<
          reactiveJ[i - n_grid_x_left][j - n_grid_y_left][1] << " " <<
          std::setw(9) << n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left] << " " <<
          std::fixed << std::setprecision(5) <<
          reactiveprobdensity[i - n_grid_x_left][j - n_grid_y_left] << std::endl;
      } else {
        outFile << std::fixed << std::setprecision(3) << (i * grid_dx) << " " <<
          (j * grid_dy) << " " <<
          std::fixed << std::setprecision(4) << "0.0000 0.0000 " <<
          std::setw(9) << n_reactiveJ[i - n_grid_x_left][j - n_grid_y_left] << " 0.00000" << std::endl;
      }
    }
  }

  outFile.close();
}