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