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

lab4.cpp

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

#define U 6
#define RMAX 100 // 300
#define K .1     // .5
#define Dx 50
#define Dy 70
#define Dz 80
#define t 1
#define e .01    // .009

// some pascal -> c compatibility defines
#define random ((double)rand() / RAND_MAX)
#define sqr(x) pow(x, 2)

double norm(double mx, double dx) {
	int n = 100, i;
	double s = 0;
	for (i = 0; i < n; ++i)
		s += random;
	return sqrt(12./n) * (s-n/2.) * dx + mx;
}

int main(void) {
	double m[6][3];
	double Ne, Sx, Sx2, Dx2, S, G, x, y, z, R, W, mx = 0, my = 0, mz = 0;
	int i, j, N = 1000;
	
	m[0][0] = 40;    m[0][1] = 40;    m[0][2] = 40;
	m[1][0] = 70;    m[1][1] = 70;    m[1][2] = 70;
	m[2][0] = 40;    m[2][1] = 80;    m[2][2] = 20;
	m[3][0] = 20;    m[3][1] = 40;    m[3][2] = 90;
	m[4][0] = 10;    m[4][1] = 50;    m[4][2] = 50;
	m[5][0] = 80;    m[5][1] = 40;    m[5][2] = 20;
	
	for (i = 0; i < U; ++i) {
		mx += m[i][0];
		my += m[i][1];
		mz += m[i][2];
	}
	
	mx /= U;
	my /= U;
	mz /= U;
	
	Sx = 0;
	Sx2 = 0;
	Ne = 10;
	
	for (N = 0; N < Ne; ++N) {
		x = norm(mx, Dx);
		y = norm(my, Dy);
		z = norm(mz, Dz);
		S = 0;
		for (j = 0; j < U; ++j) {
			R = sqrt(sqr(x-m[j][0])+sqr(y-m[j][1])+sqr(z-m[j][2]));
			if (R <= RMAX) {
				G = exp((-1./K)*pow(R/RMAX, 4));
				if (G > random) ++S;
			}
		}
		Sx += S;
		Sx2 += S*S;
		if (N >= 5)
		{
			Dx2 = Sx2/N - sqr(Sx/N);
			Ne = Dx2*sqr(t)/sqr(e);
		}
	}
	W = Sx/N;
	std::cout << "W = " << W << std::endl;
	std::cout << "N = " << N << std::endl;
}