#include "stdafx.h"
#include "Mandelbrot01.h"
#include "Mandelbrot01Dlg.h"
#include "Oblicz.h"
#include <cmath>
#include <omp.h>

double drand(double range_min, double range_max)
{
	double u = double(rand())/(RAND_MAX+1)*(range_max-range_min);
	return u;
}

void oblicz(unsigned int* tablica, int rozmiarX, int rozmiarY, int LiczbaWatkow, double Zx, double Zy, 
			double scaleX, double scaleY, int np)
{
	int idx , kolor;
	double x,y;
	double Cx,Cy;
	
	double skalaX=1.0/rozmiarX/scaleX;
	double skalaY=1.0/rozmiarY/scaleY;

	srand(time(NULL));

	#pragma omp parallel num_threads(LiczbaWatkow) default(none) shared(tablica, rozmiarX, rozmiarY, scaleX, scaleY)\
	shared(Zy,Zx,skalaX,skalaY)	private(x,y, Cx, Cy, kolor, idx)
	{
		int Xpol = int(rozmiarX/2.0+.5);
		int Ypol = int(rozmiarY/2.0+.5);
		#pragma omp for
		for (int iy=-Ypol;iy<Ypol;++iy)
		{
			Cy=Zy + iy*skalaY;
			for (int ix=-Xpol;ix<Xpol;++ix)
			{
				idx = ix+Xpol +(iy+Ypol)*rozmiarX;
				Cx = Zx + ix*skalaX;
				kolor = 0;
				for (int ip=0;ip<NP;ip++)
				{
					double dx=drand(0.0,1.0e-7);
					double dy=drand(0.0,1.0e-7);
					double Cxr=Cx+dx;
					double Cyr=Cy+dy;
					int kolorr = 0;
					x=0.0;y=0;0;
					for (int iter = 0;iter<maxiter;++iter)
					{
						double x_2=x*x; 
						double y_2=y*y;
						if (x_2+y_2>4.0) 
						{
							kolorr = 1 + iter;
							break;
						}
						y=2.0*x*y+Cyr;
						x=x_2-y_2+Cxr;
					}
					kolor += kolorr;
				}
				tablica[idx] = unsigned int(kolor/NP);
			}
		}
	}			
}

double Powierchnia(unsigned int* tablica, int rozmiarX, int rozmiarY, int LiczbaWatkow, double scaleX, double scaleY)
{
	int wsuma;	
	double dx=1.0/scaleX;
	double dy=1.0/scaleY;
	int area = 0;

	#pragma omp parallel num_threads(LiczbaWatkow) default(shared) private(wsuma)
	{
		wsuma=0;
		#pragma omp for 
		for (int idx = 0;idx<rozmiarX*rozmiarY;idx++)
			if( tablica[idx] == 0) wsuma++;
		#pragma omp critical( blok_krytyczny)
		{
			area += wsuma;
		}
	}
	return area*dx*dy/(rozmiarX*rozmiarY);
}

double Powierchnia2(unsigned int* tablica, int rozmiarX, int rozmiarY, int LiczbaWatkow, double scaleX, double scaleY)
{
	int wsuma=0;	
	double dx=1.0/scaleX;
	double dy=1.0/scaleY;

	#pragma omp parallel num_threads(LiczbaWatkow) default(none) shared(rozmiarX, rozmiarY, tablica, wsuma) 
	{
		#pragma omp for reduction(+:wsuma)
		for (int idx = 0;idx<rozmiarX*rozmiarY;idx++)
			if( tablica[idx] == 0) wsuma++;
	}
	return wsuma*dx*dy/(rozmiarX*rozmiarY);
}

void PowierchniaS(double* wynik, unsigned int* tablica, int rozmiarX, int rozmiarY, double scaleX, double scaleY)
{
	double dx=1.0/scaleX;
	double dy=1.0/scaleY;
	double faktor = dx*dy/(rozmiarX*rozmiarY);

	for (int idx=0;idx<rozmiarX*rozmiarY;idx++)
		wynik[tablica[idx]]++;
	for (int idx=0;idx<maxiter+1;idx++)
		wynik[idx]*=faktor;
	Sleep(100);
}
void Brzeg(unsigned int* brzeg, unsigned int* tablica, int rozmiarX, int rozmiarY)
{
	for (int ix=1;ix<rozmiarX-1;ix++)
	{
			for (int iy=1;iy<rozmiarY-1;iy++)
		{
			int idx = iy*rozmiarX+ix;
			int otoczenie=0;
			if(tablica[idx]==0)
			{
				otoczenie += tablica[idx-rozmiarX-1];
				otoczenie += tablica[idx-rozmiarX];
				otoczenie += tablica[idx-rozmiarX+1];
				otoczenie += tablica[idx-1];
				otoczenie += tablica[idx+1];
				otoczenie += tablica[idx+rozmiarX-1];
				otoczenie += tablica[idx+rozmiarX];
				otoczenie += tablica[idx+rozmiarX+1];
			}
			if(otoczenie>0) 
				brzeg[idx]=100;
		}
	}
	Sleep(100);
}

void Sekcje(unsigned int* brzeg, unsigned int* tablica, int rozmiarX, int rozmiarY, int LiczbaWatkow, double Zx, double Zy, 
			double scaleX, double scaleY, CString& str0,  CString& str1, CString& str2)
{
	double* pow;
	int lw;
	pow = new double[maxiter+1];
	for (int i = 0;i<maxiter+1;i++)pow[i]=0.0;
	omp_set_dynamic(0);
	#pragma omp parallel num_threads(LiczbaWatkow) default(shared)
	{
		int tid = omp_get_thread_num();
		#pragma omp sections
		{
			lw = omp_get_num_threads();
			#pragma omp section
			{
				double Start1 = omp_get_wtime();
				PowierchniaS(pow, tablica, rozmiarX, rozmiarY, scaleX, scaleY);
				double Koniec1 = omp_get_wtime();
				str0.Format(_T("Pow.: Wątek  %d z %d, czas %1.4f "),
					tid,lw,Koniec1-Start1);
				str2.Format(_T("Powierzchnia = %1.3f;"),pow[0]);
			}
			#pragma omp section
			{
				double Start1 = omp_get_wtime();
				Brzeg(brzeg, tablica, rozmiarX, rozmiarY);
				double Koniec1 = omp_get_wtime();
				str1.Format(_T("Brzeg: Wątek  %d z %d, czas %1.4f"),
					tid,lw,Koniec1-Start1);
			}
		}
	}
	delete pow;
}