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 &lt;iostream&gt;
#include &lt;vector&gt;
#include &lt;fstream&gt;
#include &lt;random&gt;
#include &lt;cmath&gt;
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&lt;double&gt; 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&lt;double&gt;  x(Particles); 
    std::vector&lt;double&gt;  pre_x(Particles);  //x_copmpnent of i'th particle position in previous step
    std::vector&lt;double&gt;  pre_y(Particles);
    std::vector&lt;double&gt;  J(N); 
    std::vector&lt;double&gt;  y(Particles); 
    for(int i = 0; i &lt; Particles; i++)  //Initializing initial angle
        {
        pre_x[i] = 0;
        pre_y[i] = 200;
        } 

    ofstream myfile;
    myfile.open ("J.txt");
    if(alph &lt; 1.01 &amp;&amp; alph &gt; 0.99)
        B = 1  //renormalazation constant for distribution
        for(int i = 0; i &lt; N; i++)  //time steps
            {
            J[i] = 0;
            for (int j = 0; j &lt; 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 &lt; y[j] || y[j] &lt; 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&lt;&lt;J[i]/Particles&lt;&lt;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.
Here are some things that may help you improve your code.

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 &lt;cmath&gt; rather than &lt;math.h&gt;

Use the new style &lt;cmath&gt; rather than the C-style &lt;math.h&gt; 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 &lt;functional&gt; 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 &amp;operator++();
    double phi() const { return phi_; }
public:
    static double get_dr(); 
    void next(const Particle &amp;other);
    bool bad() const;
    double x_;
    double y_;
    double phi_;
};


Now main can look like this:

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

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


Isn't that nice?

Simplify the mathematics

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

Particle &amp;Particle::operator++() {
    Particle testp;
    for (testp.next(*this); testp.bad(); testp.next(*this)) 
    { /* 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 &amp;other) {
    double dr = get_dr();
    phi_ = randnum(0, M_PIl);
    x_ = static_cast&lt;int&gt;(other.x_) + std::cos(phi_) * dr;
    y_ = static_cast&lt;int&gt;(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 &lt; 1.01 &amp;&amp; alph &gt; 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:

bool Particle::bad() const {
    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 &lt; 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:

bool Particle::bad() const {
    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 &lt; y_ || 
        (std::sin(A * x_) + Delta * std::sin(C * x_) / 2) * h + H &lt; 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.