#include #include #include #include #include #include #include using namespace std; // Global arrays long long S[18230000]; long long Sy[18230000]; double ex[90000]; double ex1[90000]; double fr[90000]; int main() { srand(time(0)); long long i=0; long long g=0; double frf; double dist; double l; double d1; double st1; double d0; double st0; long long f; long long f1; long long edx; long long edx1; long long m; double p; double li; double p1; double li1; long long w =10000; double en; double en1; double alpha = 0.0; double a1=0.0; double a2=0.0; double a3=0.0; double avg=0; double cn =0; double ent =0; double intr; double b1=0.0; l = 1000000; // Universe size d0 =1822.8885; // Particle 1 size d1 =1822.8885; // Particle 2 size long long kj =1000000000; // # of random throws double d0div = d0 ; intr=.05; cout << intr <<"\n"; long mk; for (mk = 2479778; mk <= 500000000; mk++) { dist = (mk) * 2 * intr+(d0); st1 =(l/2)- (((double)mk ) * intr); st0 = (l/2)+ d1 + (((double)mk ) * intr); f = 0; f1 = 0; edx = 0; edx1 = 0; en = 0.0; en1=0; ent=0; cout << l/2<<"\n"; long long kk; for (kk = 0; kk <= 20000000; kk++) { S[kk] = 0; Sy[kk] = 0; } // Next kk //number of throws long long j; for (j = 0; j <= kj; j++) { long r= rand(); double rndm=(double)r/((double)RAND_MAX); long r1= rand(); double rndm1=(double)r1/((double)RAND_MAX); long r2= rand(); double rndm2=(double)r2/((double)RAND_MAX); long r3= rand(); double rndm3=(double)r3/((double)RAND_MAX); long rx= rand(); long ry= rand(); double rndmx=(double)rx/((double)RAND_MAX); double rndmy=(double)ry/((double)RAND_MAX); p = (d0 *(rndm)); li = ((d0+dist) *rndm2); p1 = (rndm3 * d1); li1 = (rndm1 * (dist+d0)); if ( st1+p1 + li1 > st0+ p - li) { Sy[int(w*p)]=Sy[int(w*p)]+1; en1=en1+li; f1 = f1 + 1; } else { S[int(w*p)]=S[int(w*p)]+1; en=en+li; f = f + 1; } } // Next j frf = (double)f/en; //energy of the particle en1 = (double)f1/en1; long long n; for (n = 0; n <= int(w*d1); n++) { // CALCULATE EXPECTATION VALUE AFTER INTERACTION HAS TAKEN PLACE edx = edx + (( n) * S[n]); edx1 = edx1 + (( n) * Sy[n]); } ex1[mk] = (double)edx1 / ((double)f1)- (0.5 * int(w*d1))+.5 ; ex[mk] = (double)edx / ((double)f)- (0.5 * int(w*d1))+.5 ; cout << std::setprecision(17)<< dist<<" "<