SHARE
TWEET

Untitled

a guest Oct 13th, 2017 62 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2.  #include <vector>
  3.  #include <fstream>
  4.  #include <random>
  5.  #include <functional>
  6.  #include <math.h>
  7.  using namespace std;
  8.  double randnum (double aa, double bb)  //defining a function to create random numbers
  9.  {
  10.  static std::default_random_engine generator;
  11.  std::uniform_real_distribution<double> distribution (aa,bb);
  12.  return distribution(generator);
  13.  }
  14.  
  15. int main() {
  16. double dr;
  17. double a = 1.0;              
  18. number with power-law distribution
  19. int N = 10000;      
  20. double u;                    
  21. //random number with desired distribution
  22. int Particles = 100000;        
  23. int H = 500;                  
  24. int h = 200;              
  25. double C = 0.004;            
  26. int pre_y, pre_x;
  27. double A = 0.002;        
  28. double Delta = 1/7;  
  29. double b = 25.0;                
  30. number with power-law distribution
  31. std::vector<double>  phi(Particles);
  32. std::vector<double>  x(Particles);
  33. std::vector<double>  J(N);
  34. std::vector<double>  y(Particles);
  35. for(int i = 0; i < Particles; i++)  
  36.     {
  37.     x[i] = 0;
  38.     y[i] = 200;
  39.     }
  40. ofstream myfile;
  41. myfile.open ("J.txt");
  42. int alph = 0;              
  43. double B = (1-alph)/( pow( b,(1-alph))- pow( a,(1-alph)));
  44. if(alph < 1.01 && alph > 0.99)
  45. B = 1/(log(b)* pow( b, 1-alph )-log(a)* pow( a, 1-alph ));
  46.  
  47. for(int i = 0; i < N; i++)  //steps
  48. {
  49.         J[i]=0;
  50.         for (int j = 0; j < Particles; j++)  //Particles
  51.             {
  52.              u = randnum(0,1);
  53.              dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
  54.              phi[j] = randnum(0,M_PIl);
  55.              pre_x = x [j];
  56.              pre_y = y [j];
  57.              x[j] = pre_x + cos(phi[j]) * dr;
  58.              y[j] = pre_y + sin(phi[j]) * dr;
  59.             while( (sin(A * x[j]) + Delta * sin(C * x[j])/2) * h + H < y[j])
  60.                  {
  61.                  u = randnum(0,1);
  62.                  dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
  63.                  phi[j] = randnum(0,M_PIl);
  64.                  x[j] = pre_x + cos(phi[j]) * dr;
  65.                  y[j] = pre_y + sin(phi[j]) * dr;
  66.                  }
  67.                  J[i] = J[i] + cos(phi[j]);
  68.              }
  69.                  myfile << J[i]/Particles << endl;
  70.         }
  71.         myfile.close();
  72.  }
RAW Paste Data
Top