Simulating motion of particles

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.

• It is simulating motion of particles @Incomputable Oct 13 '17 at 7:19
• I've edited the title and fixed a stray comment, but if the indentation doesn't reflect exactly what's in your IDE then I'd suggest you re-paste the code block straight from the IDE - then select the whole block and press Ctrl+K to automatically add the leading 4 spaces that will make it render as a code block. If the code is indented that way in your editor, ...well that makes a low-hanging fruit for reviewers to pick =) Oct 13 '17 at 13:44
• Feel free to further edit your introduction to present your code to reviewers - tell us about it, what it does / what problem it's solving. See how to get the best value out of Code Review for more information. Oct 13 '17 at 13:46
• @coddercode, I fixed all of the problems with comments. The code does compile now, although I wasn't sure enough to run it on my machine. Please fill in the dots that I left after edit, to describe the code. You can use edit button at the bottom of the post for that. After the filling those in, the post should be ready to go! It might look like it is cumbersome at first, but that is usually how it would be done in a more formal world. Oct 13 '17 at 13:51
• I have added the description. What should I do to find someone who optimizes my code? @Incomputable Oct 13 '17 at 13:56

Don't abuse using namespace std

Putting using namespace std at the top of every program is a bad habit that you'd do well to avoid. Once you remove that statement, you will have to add in the namespace where appropriate. For example, fstream would be changed to std::fstream, endl to std:endl, etc.

Break up the code into smaller functions

In the current code, almost everything is done in main(). Rather than having everything in one long function, it would be easier to read and maintain if each discrete step were its own function.

Eliminate "magic numbers"

This code is littered with "magic numbers," that is, unnamed constants such as 0.004, 200, 25.0, etc. Generally it's better to avoid that and give such constants meaningful names. That way, if anything ever needs to be changed, you won't have to go hunting through the code for all instances of "200" and then trying to determine if this particular 200 is one that needs to be changed or just some other constant that happens to have the same value.

Consider improving names

I think randnum is not a bad function name, but most of the variable names are single letter names that are not at all descriptive of what they're intended to represent. More descriptive variable names would help readers of the code understand what it's doing.

Write results to the std::cout

Rather than have a single hard-coded file name, write the data to std::cout. If the user then wants to save a copy, doing so is then a simple matter of command line output redirection.

Use <cmath> rather than <math.h>

Use the new style <cmath> rather than the C-style <math.h> for two reasons. First, it is more idiomatic modern C++, but also because it uses namespaces. Consistent with that, you shoud also use std::cos rather than just cos. It's the same for log, pow and sin.

Don't use std::endl when '\n' will do

Using std::endl emits a \n and flushes the stream. Unless you really need the stream flushed, you can improve the performance of the code by simply emitting '\n' instead of using the potentially more computationally costly std::endl.

Use only necessary #includes

The #include <functional> line is not necessary and can be safely removed.

Use constexpr where practical

In main, almost all of the variables are actually used as constants, so it would make sense to at least declare them as const and preferably constexpr if you're using a C++11 compliant compiler (and if you're not, you really ought to update your compiler).

Use objects

The purpose of the code is to simulate objects called particles. It's written in C++. Why doesn't the code use C++ objects named Particle? Here's a Particle class we can use:

class Particle {
public:
Particle() : x_{0}, y_{200}, phi_{0} {}
Particle &operator++();
double phi() const { return phi_; }
public:
static double get_dr();
void next(const Particle &other);
double x_;
double y_;
double phi_;
};

Now main can look like this:

int main() {
constexpr int N{10000};
std::vector<double> J(N, 0);
constexpr int numParticles{1000};
std::vector<Particle> particles(numParticles);

for (auto &step : J) {
for (auto &p : particles) {
++p;
step += std::cos(p.phi());
}
std::cout << step / numParticles << '\n';
}
}

Isn't that nice?

Simplify the mathematics

Let's look at how the class functions above are implemented. First is the operator++:

Particle &Particle::operator++() {
Particle testp;
{ /* do nothing */ }
return *this = testp;
}

This has the effect of the inner loop in the original code. The next function is this:

void Particle::next(const Particle &other) {
double dr = get_dr();
phi_ = randnum(0, M_PIl);
x_ = static_cast<int>(other.x_) + std::cos(phi_) * dr;
y_ = static_cast<int>(other.y_) + std::sin(phi_) * dr;
}

The get_dr function is this:

double Particle::get_dr() {
constexpr double a = 1.0;
constexpr double b = 25.0;
constexpr int alph = 0;
constexpr double B = (alph < 1.01 && alph > 0.99) ?
1 / (std::log(b) * std::pow(b, 1 - alph) - std::log(a) * std::pow(a, 1 - alph))
: (1 - alph) / (std::pow(b, (1 - alph)) - std::pow(a, (1 - alph)));
return std::pow(std::pow(a, 1 - alph) + randnum(0, 1) * (1 - alph) / B, 1 / (1 - alph));
}

Finally the simple version of bad() is this:

constexpr int H = 500;
constexpr int h = 200;
constexpr double C = 0.004;
constexpr double A = 0.002;
constexpr double Delta = 1 / 7;
return
(std::sin(A * x_) + Delta * std::sin(C * x_) / 2) * h + H < y_;
}

However, the maximum value of sin is 1, so the upper bound of the left hand side of that inequality can be expressed as a constant. Since comparing against a constant is usually much faster than invoking std::sin we can use this knowledge as follows:

constexpr int H = 500;
constexpr int h = 200;
constexpr double C = 0.004;
constexpr double A = 0.002;
constexpr double Delta = 1 / 7;
constexpr double maxval = (1 + Delta * 1 / 2) * h + H;
return maxval < y_ ||
(std::sin(A * x_) + Delta * std::sin(C * x_) / 2) * h + H < y_;
}

The short circuit expression evaluation means we only invoke the computationally expensive sin function is the comparison against maxval fails. Testing on my machine shows that this makes the program about 20% faster. You could squeeze out still more time by noting that x need only be calculated if the first test fails. Generally, you want to arrange things so that simpler tests are done first.

It may be advantageous to look at using a different coordinate scheme (polar vs. rectangular for instance) if the resulting equations simplify. I suspect they might.

• When I remove using namespace std, the code doen't run and there are errors such as: ‘ofstream’ was not declared in this scope Oct 14 '17 at 3:49
• COuld you please update your answer? There was a problem in my previous code and I have changed it know Oct 14 '17 at 3:52
• Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers. Oct 14 '17 at 14:43
• ok I don't do that anymore. Have updated your answer according to edited question? Oct 14 '17 at 16:55
• No, but all of the principles outline in my answer apply to the code old or new. Oct 14 '17 at 17:25

Your code is very hard to read. This is mostly because most of your variables have cryptic names.

In a physics simulation, it's ok to name natural constants such as c or G by their name, but that's not what you do here. How is any reader going to tell h and H apart? Why are some variables uppercase and most others lowercase? What does each of the variables mean?

Performance and std::endl don't go together.

Comparing an int variable with 0.99 is nonsense.