#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 = int( 1.0 + iter - log(log(sqrt(x_2+y_2)))/log(2.0));
							break;
						}
						y=2.0*x*y+Cyr;
						x=x_2-y_2+Cxr;
					}
					kolor += kolorr;
				}
				tablica[idx] = unsigned int(kolor/NP);
			}
		}
	}			
}

double Powierzchnia(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 Powierzchnia2(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);
}