#include <iostream>
 #include <vector>
 #include <fstream>
 #include <random>
 #include <functional>
 #include <math.h>
 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;
number with power-law distribution
int N = 10000;
double u;
//random number with desired distribution
int Particles = 100000;
int H = 500;
int h = 200;
double C = 0.004;
int pre_y, pre_x;
double A = 0.002;
double Delta = 1/7;
double b = 25.0;
number with power-law distribution
std::vector<double>  phi(Particles);
std::vector<double>  x(Particles);
std::vector<double>  J(N);
std::vector<double>  y(Particles);
for(int i = 0; i < Particles; i++)
    {
    x[i] = 0;
    y[i] = 200;
    }
ofstream myfile;
myfile.open ("J.txt");
int alph = 0;
double B = (1-alph)/( pow( b,(1-alph))- pow( a,(1-alph)));
if(alph < 1.01 && alph > 0.99)
B = 1/(log(b)* pow( b, 1-alph )-log(a)* pow( a, 1-alph ));

for(int i = 0; i < N; i++)  //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[j] = randnum(0,M_PIl);
             pre_x = x [j];
             pre_y = y [j];
             x[j] = pre_x + cos(phi[j]) * dr;
             y[j] = pre_y + sin(phi[j]) * dr;
            while( (sin(A * x[j]) + Delta * sin(C * x[j])/2) * h + H < y[j])
                 {
                 u = randnum(0,1);
                 dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
                 phi[j] = randnum(0,M_PIl);
                 x[j] = pre_x + cos(phi[j]) * dr;
                 y[j] = pre_y + sin(phi[j]) * dr;
                 }
                 J[i] = J[i] + cos(phi[j]);
             }
                 myfile << J[i]/Particles << endl;
        }
        myfile.close();
 }