// system parameters double mu; // mobility double Pe; // peclet number double v_drive; // effetive velocity double tr; // tumble rate double rr; // run rate int runsteps; // number of steps in current run int tumblesteps; // number of steps in current tumble int seed; int histosize; int runsteps_back; // number of steps in current run int tumblesteps_back; // number of steps in current tumble int steps_ind; // TODO global int steps_min; // int steps_ind_new; // int steps_min_new; // bool in_target; bool tumble; // system in tumble phase or not? bool run; // system in run phase or not? double runcond; // does is probability larger then this condition for run double runprob; // probability of run double tumblecond; double tumbleprob; double Tinit; // temparature double kx; // potential parameter double ky; // potential parameter double x0; // potential parameter double D; // diffusion double D_th; // rotational diffusion double B; // sqrt(2 D dt) double Bth; // sqrt(2 D_th dt) // simulation parameters int Ngrid; // number of grid cells int Nsim; // maximum steps in a trajectory double dt; // time-step length MD // trajectory parameters double x_one; // x, y positions of the particles double y_one; double theta_one; double vartheta_one; double fx_one; // x, y components of force on the particles double fy_one; double nu_one; // phase, 1 = run, 2 = tumble double U_one; double ** R; // old trajectory, (x,y,theta, fx,fy,nu) int ind; // randomly chosen index / shooting point bool forward; // forward shooting hit the target bool backward; // backward shooting hit the reactant bool R_found; // checks if shooting finds new trajectory double ** R_new; // new trajectory, (x,y,theta, fx,fy,nu) int ** n_reactiveJ; double ** times; double ** reactiveprobdensity; float *** reactiveJ; int *Histo; // lengthdist int n_min; // index of step when the particle leaves the reactant of old trajectory int n_min_new; // index of step when the particle leaves the reactant of new trajectory int n_max; // index of step when the particle reachtes the target of old trajectory int n_max_new; // index of step when the particle reachtes the target of new trajectory int indexshift; int path_length; //length of old path int path_length_new; //length of new path double mean_length; //average path length double correlation_lifespan; //1/e zeit der korrelationsfunktion double acceptance_ratio; bool accept; // steady state distribution int Neq; // Number of equilibrations steps, number of simulation steps int eq_rounds; // Number of times we do Neq for steady state float *** grid_count_run; // grid of positions that count if particle is in cell float *** grid_count_tumble; // grid of positions that count if particle is in cell double xmin; // grid around the reactant basin double xmax; double deltx; double ymin; double ymax; double delty; double grid_x_left , grid_x_right, grid_y_left , grid_y_right , grid_dy , grid_dx ; int n_grid_x_left ,n_grid_x_right ,n_grid_y_left ,n_grid_y_right; double thetamin; double thetamax; double deltatheta; int index_x; // index of grid count int index_y; int index_theta; long int md_time; long int tps_time; int NumOfTumbles; int NumOfRuns; int numRows; bool inftyrunrate; int sampsize_comparetimes; // Declare global vectors to store times // fucntions class RanMars *random_mars; void md_step_RTPE(); float* read_pathlengths(double peclet, double tumblerate, double runrate, int sampsize); double* PathLengthCorreltaion(double peclet, double tumblerate, double runrate, int sampsize); void write_trajectory_to_file(FILE * outFile); void sampling_TPS_trajectories(double peclet, double tumblerate, double runrate, int sampsize); void length_correltaion(double dt, double pe, double tr, double rr, double samplesize); void test_output_for_read_gridcount(); void read_steady_state_dist(); void SteadyState_distribution(); void print_md_trajectory(); void print_TPS_trajectories(); void compare_time(); void sampling(double peclet, double tumblerate, double runrate, int sampsize, bool TPS, bool BF, bool steadystate); double transition_prob_n_or_more(double nu, int Nmin); void accept_or_reject(); void new_R(); void lattice_packing(); void save_R(); // TODO rn i copy the rng seed from find reactive trajectory. maybe more elegant void find_reactive_R(); void find_first_reactive_R(); double calculate_potential(double x, double y); double transition_prob(double nu_i, double nu_ii); double transition_prob_back(double nu_i, double nu_ii); void md_step(); void md_step_back(int i); //for tps and saving trajectory void md_step_forward(int i); //for tps and saving trajectory void runstep(); void tumblestep(); void calculate_forces(double x, double y); int phaselength(double phase, int min); void print_trajectories(int step, char folder[256]); double calculateMean(const std::string& file_path); void brute_force_lengthdist(int steps); void brute_force_lengthdist_rtpe(int steps); void TPS_lengthdist(int steps); void go_to_R(); void write_times(); void sample_brute_force_trajectories(); void sample_TPS_trajectories(); void calcReactiveFlux(); void compute_final_averages2d(); void print_reactiveprobabilitydensity2d(); void print_reactivecurrent2d();