#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <algorithm>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
         // constnatns (attractive)
#define jump 0.05//this is how many time steps it takes to change alpha,for the radicle jump this is twice this size
#define time_jump 100
#define N 100 //number of particles
#define time_end  10000
#define ensambles 1

//int time_end=100000;
using namespace std;
typedef boost::mt19937 base_generator_type;
typedef boost::variate_generator< base_generator_type&, boost::uniform_int<> > gen_type;
typedef boost::variate_generator< base_generator_type&, boost::uniform_real<> > gen_type1;

             //defning the random numbers generators
          unsigned int s=time(NULL);
          base_generator_type generator(s);
          boost::uniform_real<> uni_distq(0,1);
          gen_type1 uniq( generator,  uni_distq );
          boost::uniform_int<> uni_dist(0,N-1);
          gen_type uni( generator,  uni_dist );

void dissipation(double energy[][5][time_end],int n1,int n2,int i,int t,double alpha){
                        double temp;
                        temp=energy[n1][i][t]+energy[n2][i][t];
                        double z=uniq();
                        energy[n1][i][t]=z*alpha*temp;
                        energy[n2][i][t]=(1-z)*alpha*temp;
                        for(int k=0;k<N;k++){if(k!=n1 &&k!=n2)
                                                      energy[k][i][t]+=(1-alpha)*temp/(N-2);
                                                      }
     }
void print_alpha(double energy[][5][time_end],int i,int t){
             string f_type = "alpha.txt";
             char f_name[8];
             sprintf(f_name,"%02d",i);
             string file = f_name + f_type;
             ofstream myfile ( file.c_str() );
             for(int k=0;k<N;k++)
             myfile<<(double) energy[k][i][t]<<" ";
             myfile.close();
     }
void alpha_change(double *alpha,t){
                 if((max(t,100) % min(t,100))==0){
                                  *alpha[0]+=jump;
                                  *alpha[2]=uniq();
                                  }
                 if((max(t,150) % min(t,150))==0){
                                  cout<<t<<"  the current time  "<<endl;
                                  *alpha[1]+=0.06*jump;
                                  *alpha[3]=0.3*alpha[2];
                                  for(int g=0;g<5;g++)
                                          print_alpha(energy,g,t);
                                          }
int main(int argc, char *argv[]){

    for(int n=0; n<ensambles ; n++) { //ensambles
            //prepering the paramters and the files output;
          double energy[N][5][time_end];  /* 0&1-dist slow/fast,2&3-rand slow/fast,4-boltzman*/
          double alpha[4];
           //giving them the inital energy
           for(int j=0;j<5;j++){
                   alpha[j]=0;
                   for(int i=0;i<N;i++)
                            energy[i][j][0]=1;
                           }
    for(int t=0;t<time_end;t++)
    {      //time
            //getting the particles
            int n1=uni();
            int n2=uni();
            if(n1==n2){
                        if(n1==N)
                        n2=N-uni();
                        else
                        n2=abs(n1-N);
                        }
            //the disspation procceses
            for(int i=0;i<5;i++){
                    dissipation(energy,n1,n2,i,t,alpha[i]);
                    }
            //changing the alpha
            alpha_change(alpha,t);

            }  //time

                                     }  //ensambles
               system("pause");
               return EXIT_SUCCESS;

               } ;