Description
The code does particle simulation. In this code, some particles take steps with a certain distribution length. if particles go beyond a curve, the code put them back.
Code
#include <iostream>
#include <vector>
#include <fstream>
#include <random>
#include <cmath>
using namespace std;
double randnum (double aa, double bb) //defining a function to create random numbers
{
static std::default_random_engine generator;
std::uniform_real_distribution<double> distribution (aa,bb);
return distribution(generator);
}
int main() {
double dr;
double a = 1.0; //minimum value for generated random number with power-law distribution
int N = 100000; //number of time steps
double u; //uniform random number used to create random number with desired distribution
int Particles = 1000; //number of particles
int H = 500; //potential parameter
int h = 200; //potential parameter
double C = 0.004; //potential parameter
double A = 0.002; //Potential Parameter
double Delta = 1/7; //Assymetric Parameter
double b = 25.0; //maximum value for generated random number with power-law distribution
double phi; //first angle
std::vector<double> x(Particles);
std::vector<double> pre_x(Particles); //x_copmpnent of i'th particle position in previous step
std::vector<double> pre_y(Particles);
std::vector<double> J(N);
std::vector<double> y(Particles);
for(int i = 0; i < Particles; i++) //Initializing initial angle
{
pre_x[i] = 0;
pre_y[i] = 200;
}
ofstream myfile;
myfile.open ("J.txt");
if(alph < 1.01 && alph > 0.99)
B = 1 //renormalazation constant for distribution
for(int i = 0; i < N; i++) //time steps
{
J[i] = 0;
for (int j = 0; j < Particles; j++) //Particles
{
u = randnum(0,1);
dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
phi = randnum(0,M_PIl);
x[j] = pre_x[j] + cos(phi) * dr;
y[j] = pre_y[j] + sin(phi) * dr;
while( (sin(A * x[j]) + Delta * sin(C * x[j])/2) * h + H < y[j] || y[j] < 0)
{
u = randnum(0,1);
dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
phi = randnum(0,M_PIl);
x[j] = pre_x[j] + cos(phi) * dr;
y[j] = pre_y[j] + sin(phi) * dr;
}
pre_x[j] = x[j];
pre_y[j] = y[j];
J[i] = J[i] + cos(phi);
}
myfile<<J[i]/Particles<<endl; //Outputs array to txtFile
}
myfile.close();
}
Concerns:
code seems slow
general code quality
It is slow or in other words, I think is possible to make it faster. I need both "N" and "Particles" to be large values. With this condition, how can I make the code to run faster? Is there any change which could improve the run-time? Any other feedback is welcome as well.