Лабораторная №2

lab2.cpp

#include <cmath>
#include <cstdlib>
#include <iostream>

/**
 * Test data:
 * 70 120 .5 .7 50 70 20 10 30
 * W1 = 0.292
 * W2 = 0.145
 */

/**
 * if 1, distribution is normal,
 * else it is uniform
 */

#define SEQ1_NORMAL 0
#define SEQ2_NORMAL 1

#define EPS 0.01

double rselect (double mx, double dx, int isnormal) {
	// mx is mean value, dx is variance
	if (isnormal) {
		double S = 0;
		int i, n = 100;
		for (i = 0; i < n; ++i)
			S += (double)rand() / RAND_MAX;
		return sqrt (12.0/n)*(S-n/2.0)*dx+mx;
	}
	else return mx+dx*(2.0*rand() / RAND_MAX - 1);
}

int main (void) {
	double Rmax, K, dx, dy, dz, Romax, Ko, f1, f2, ta = 1;
	double S, f, G, x, y, z, R, W;
	double P, Emax, Go, Gf, r, e, M, mz, alpha;

	int i, N = 100000;
	std::cout << "Enter Rmax, Romax, K, Ko, dx, dy, dz, f1, f2:" << std::endl;
	std::cin >> Rmax >> Romax >> K >> Ko >> dx >> dy >> dz >> f1 >> f2;
	S = 0;
	i = 0;
	
	/*** PART 1 ***/
	// calculating an integral
	while (i < N) {
		#if SEQ1_NORMAL
			x = rselect(0, dx, 1);
			y = rselect(0, dx, 1);
			z = rselect(0, dx, 1);
		#else
			x = rselect(0, dx*3, 0);
			y = rselect(0, dy*3, 0);
			z = rselect(0, dz*3, 0);
		#endif
		f = pow(2*M_PI, -1.5) / (dx*dy*dz);
		f *= exp(-0.5*((x/dx)*(x/dx)+(y/dy)*(y/dy)+(z/dz)*(z/dz)));
		R = sqrt(x*x+y*y+z*z);
		
		// coordinate rule
		if (R <= Rmax)
			G = exp((-1.0/K)*(R/Rmax)*(R/Rmax)*(R/Rmax)*(R/Rmax));
		else G = 0;
		
		#if SEQ1_NORMAL
			S += G;
		#else
			S += f*G;
		#endif
		i++;
	}
	#if SEQ1_NORMAL
		W=S/N;
	#else
		W=216.0*dx*dy*dz*S/N;
	#endif
	if (W > 1) W=1;
	std::cout << "W1 = " << W << std::endl;
	
	/*** PART 2 ***/
	Emax = M_PI * (f2-f1) / 360;
	S = 0; N = 0; M = 1000;
	while (N < M) {
		++N;
		x = rselect(0, dx, SEQ2_NORMAL);
		y = rselect(0, dy, SEQ2_NORMAL);
		r = sqrt(x*x+y*y);
		f = M_PI*0.5*(f1+f2)/180;
		mz = r*cos(f)/sin(f);
		z = rselect(mz, dz, SEQ2_NORMAL);
		R = sqrt(x*x+y*y+z*z);
		alpha = atan(r/z);
		e = fabs(alpha-f);
		
		// coordinate rule for f part
		if (R <= Rmax)
			Gf = exp((-1.0/K)*(R/Rmax)*(R/Rmax)*(R/Rmax)*(R/Rmax));
		else Gf = 0;
		
		// coordinate rule for o part
		if (R <= Romax && e < Emax && z >= 0)
			Go = exp((-1.0/Ko)*(R/Romax)*(R/Romax)*(R/Romax)*(R/Romax))
			* exp(-(e/Emax)*(e/Emax)*(e/Emax)*(e/Emax));
		else Go = 0;
		
		if (Gf > Go) S += Gf;
		else S += Go;
		
		P = S/N;
		if (N > 100) M = (ta/EPS)*(ta/EPS)*P*(1-P);
	}
	W = S/N;
	std::cout << "W2 = " << W << std::endl;
	return 0;
}