// Kod opracowany na potrzeby Laboratorium 1:
// Generacja syntetycznych fluktuacji turbulentnych na wlocie do obszaru obliczeniowego dla metod LES i metod hybrydowych RANS-LES
// Autor: Bartosz Górecki

// Zalacz potrzebne biblioteki
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <iostream>
#include <time.h>
#include "winbgi2.h"
#include "Wykres.h"

// Zdefiniuj potrzebne stale
#define amp 1.452762113
#define pi 3.141592653589793

// Funkcje skladowe programu

void WydrukujParametry(int nmodes, double wew1fct, double visc, double qm, double sli, double up, double epsm, int nstep)
{
	// Wydrukuj ustawienia generatora fluktuacji
	printf("\t**************************************************\n");
	printf("\t GENERACJA SYNTETYCZNYCH FLUKTUACJI TURBULENTNYCH\n");
	printf("\t**************************************************\n");
	printf("\tAutor: Bartosz Gorecki\n\n");

	printf("Liczba modow:\t\t\t\t\t%d\n", nmodes);
	printf("Stosunek skal k1/ke:\t\t\t\t%lf\n", wew1fct);
	printf("Lepkosc kinematyczna:\t\t\t\t%lg\n", visc);
	printf("Energia kinetyczna fluktuacji:\t\t\t%lf\n", qm);
	printf("Skala dlugosci:\t\t\t\t\t%lf\n", sli);
	printf("Odchylenie standardowe skladowej predkosci:\t%lf\n", up);
	printf("Epsilon:\t\t\t\t\t%lf\n", epsm);
	printf("Liczba realizacji:\t\t\t\t%d\n", nstep);
}

void GenerujSiatke(double *&wsp, double d, int n)
{
	// Funkcja generuje rownomierna siatke wezlow w 1D

	wsp = new double[n];
	wsp[0] = 0;
	for(int i = 1; i<n; ++i)
		wsp[i] = wsp[i-1] + d;
}

double MinimalnaOdleglosc(double *wsp, int n)
{
	// Funkcja znajduje minimalna odleglosc miedzy wezlami siatki 1D przechowywanej w wektorze wsp

	double dmin = 1e6;
	for(int i = 0; i<n-1; ++i)
	{
		if(dmin > (wsp[i+1] - wsp[i]))
			dmin = wsp[i+1] - wsp[i];
	}

	return dmin;
}

void LosujKaty(double *fi, double *psi, double* alfa, double* teta, int nmodes)
{
	// Funkcja generuje losowe katy o rownomiernym rozkladzie liczb losowych
	
	double ang;

	// Wygeneruj fi, psi oraz alfa z przedzialu od 0 do 2*pi
	for(int i = 0; i<nmodes; ++i)
	{
		fi[i] = ((double)rand())/((double)RAND_MAX)*2*pi;
		psi[i] = ((double)rand())/((double)RAND_MAX)*2*pi;
		alfa[i] = ((double)rand())/((double)RAND_MAX)*2*pi;
		ang = ((double)rand())/((double)RAND_MAX);
		teta[i] = acos(1.-ang/0.5);
	}
}

int main(int argc, char* argv[])
{
	// Deklaracje zmiennych potrzebnych w programie
	int step;
	double dminy, dminz, dxmin;
	double wnrn, wnre, wnreta, wnr1, dkl;
	double xc, yc, zc, utrp, vtrp, wtrp;
	double kxi, kyi, kzi, sx, sy, sz, kx, ky, kz, rk;
	double arg, tfunk, e, utn;

	// Zdefiniuj parametry i ustawienia
	int nmodes = 150;
	double wew1fct = 2;
	double visc = 15.2e-6;
	double qm = 4;
	double delta = 0.5; // zakladamy grubosc warstwy przysciennej rowna 0.5
	double sli = 0.1*delta; // uzywamy skali dlugosci rownej 1/10 grubosci warstwy przysciennej
	double up = sqrt(2*qm/3);
	double epsm = (pow(qm,1.5))/sli;

	int nstep = 200; // wygeneruj nstep realizacji

	WydrukujParametry(nmodes, wew1fct, visc, qm, sli, up, epsm, nstep);
	
	// Zdefiniuj siatke
	int nx = 2; // liczba wezlow siatki w kierunku x
	int ny = 2; // liczba wezlow siatki w kierunku y
	int nz = 63; // liczba wezlow siatki w kierunku z
	double dx = 1;
	double dy = 1;
	double dz = 0.004;

	// Generuj siatke w kierunku odpowiednio x, y i z
	double *x, *y, *z;

	GenerujSiatke(x, dx, nx);
	GenerujSiatke(y, dy, ny);
	GenerujSiatke(z, dz, nz);

	// Znajdz minimalna odleglosc miedzy wezlami siatki

	dminy = MinimalnaOdleglosc(y, ny);
	dminz = MinimalnaOdleglosc(z, nz);
	dxmin = std::min(dminy, dminz);

	// Zainicjalizuj generator liczb losowych
	srand(time(NULL));

	// Zaalokuj i wyzeruj tablice
	double *wnr, *alfa, *fi, *teta, *psi, *wnrf, *dkn, *kxio, *kyio, *kzio;
	double *sxio, *syio, *szio, *u, *v, *w;

	wnr = new double[nmodes];
	alfa = new double[nmodes];
	fi = new double[nmodes];
	teta = new double[nmodes];
	psi = new double[nmodes];
	wnrf = new double[nmodes+1];
	dkn = new double[nmodes];
	kxio = new double[nmodes];
	kyio = new double[nmodes];
	kzio = new double[nmodes];
	sxio = new double[nmodes];
	syio = new double[nmodes];
	szio = new double[nmodes];
	u = new double[nz];
	v = new double[nz];
	w = new double[nz];

	// Zadeklaruj dwuwymiarowe tablice na przechowanie pol zaburzen dla kazdego kroku czasowego
	double **ut, **vt, **wt;

	ut = new double*[nz-1];
	vt = new double*[nz-1];
	wt = new double*[nz-1];

	for(int i = 0; i<nz-1; ++i)
	{
		ut[i] = new double[nstep];
		vt[i] = new double[nstep];
		wt[i] = new double[nstep];
	}
	// Koniec deklaracji dwuwymiarowych tablic

	srand(time(NULL));
	
	for(int it = 0; it<nstep; ++it)
	{
		// Generuj losowe katy
		LosujKaty(fi, psi, alfa, teta, nmodes);

		// Najwyzsza liczba falowa
		wnrn = 2*pi/dxmin;

		// Liczba falowa zwiazana z modem o najwiekszej energii
		wnre = 9*pi*amp/(55*sli);

		// Liczba falowa uzyta w czlonie lepkim w spektrum von Karmana
		wnreta = pow((epsm/pow(visc, 3.0)), 0.25);

		// Najmniejsza liczba falowa
		wnr1 = wnre/wew1fct;

		// Krok miedzy liczbami falowymi
		dkl = (wnrn-wnr1)/nmodes;

		// Wyznacz rownomierny ciag liczb falowych
		for(int m = 0; m<=nmodes; ++m)
			wnrf[m] = wnr1+dkl*(double)m;

		// Wyznacz liczby falowe w srodkach podzialow
		for(int m = 0; m<nmodes; ++m)
		{
			wnr[m] = 0.5*(wnrf[m] + wnrf[m+1]);
			dkn[m] = wnrf[m+1] -wnrf[m];
		}

		// Wektory falowe z losowych katow
		for(int m = 0; m<nmodes; ++m)
		{
			kxio[m] = sin(teta[m])*cos(fi[m]);
			kyio[m] = sin(teta[m])*sin(fi[m]);
			kzio[m] = cos(teta[m]);

			// sigma (s = sigma) z losowych katow. Sigma jest wersorem, ktory
			// wyznacza kierunek syntetycznego wektora predkosci (u, v, w)
			sxio[m]=cos(fi[m])*cos(teta[m])*cos(alfa[m])-sin(fi[m])*sin(alfa[m]);
			syio[m]=sin(fi[m])*cos(teta[m])*cos(alfa[m])+cos(fi[m])*sin(alfa[m]);
			szio[m]=-sin(teta[m])*cos(alfa[m]);
		}

		// Petla po siatce

		for(int k = 0; k<nz-1; ++k)
		{
			for(int j = 0; j<ny-1; ++j)
			{
				// Wspolrzedne srodka komorki
				xc = 0.5;
				yc = 0.5*(y[j] + y[j+1]);
				zc = 0.5*(z[k] + z[k+1]);

				// Zainicjalizuj fluktuacje turbulentne zerem
				utrp = 0;
				vtrp = 0;
				wtrp = 0;

				// Petla po wszystkich liczbach falowych
				for(int m = 0; m<nmodes; ++m)
				{
					kxi = kxio[m];
					kyi = kyio[m];
					kzi = kzio[m];

					sx = sxio[m];
					sy = syio[m];
					sz = szio[m];

					kx = kxi*wnr[m];
					ky = kyi*wnr[m];
					kz = kzi*wnr[m];
					rk = sqrt(kx*kx + ky*ky + kz*kz);

					// Jesli liczba falowa rk jest mniejsza niz najwieksza dopuszczalna,
					// wygeneruj fluktuacje
					if(rk < wnrn)
					{
						arg = kx*xc + ky*yc + kz*zc + psi[m];

						tfunk = cos(arg);

						// Spektrum von Karmana
						e = amp/wnre*pow(wnr[m]/wnre,4.0)/(pow(1+pow(wnr[m]/wnre,2.0),17./6.))*exp(-2*pow(wnr[m]/wnreta,2.0));
						//printf("e = %e\n", e);
						utn = sqrt(e*up*up*dkn[m]);

						// Syntetyczne pole predkosci
						utrp += 2*utn*tfunk*sx;
						vtrp += 2*utn*tfunk*sy;
						wtrp += 2*utn*tfunk*sz;
					}
				}

				// Przechowaj fluktuacje w tablicach
				
				// Zdefiniuj, ktora z realizacji chcesz przechowac w tablicach u, v, w
				if(it == 5)
				{
					u[k] = utrp;
					v[k] = vtrp;
					w[k] = wtrp;
				}

				ut[k][it] = utrp;
				vt[k][it] = vtrp;
				wt[k][it] = wtrp;
			}
		}
	}
	
	graphics(800, 500);

	// Zapisz realizacje z kolejnych krokow czasowych do pliku (tylko skladowa u_prim) //////////////////////////////////////
	FILE *f = fopen("Bez_korelacji_w_czasie.xls", "w");

	for(int i = 0; i<nz-1; ++i)
	{
		fprintf(f, "%lf\t", z[i]);
		for(int j = 0; j<nstep; ++j)
		{
			fprintf(f, "%lf\t", ut[i][j]);
		}
		fprintf(f, "\n");
	}

	fclose(f);
	// Koniec zapisu /////////////////////////////////////////////////////////////////////////////////////////////////////////
	
	// CWICZENIA

	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// Ad. 1. Obejrzyj spektrum von Karmana
	
	double *spektrum, *logomega;
	spektrum = new double[nmodes];
	logomega = new double[nmodes];

	for(int i = 0; i<nmodes; ++i)
	{
		// Uzupelnij: spektrum[i] = ...
	}
	// Uzupelnij: Rysuj wykres ze skala liniowa Wykres(...);

	for(int i = 0; i<nmodes; i++)
	{
		// Uzupelnij: logomega[i] = ...
		// Uzupelnij: spektrum[i] = ...
	}
	// Uzupelnij: Rysuj wykres ze skala logarytminczna Wykres(...)

	delete [] spektrum;
	delete [] logomega;
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// Ad. 2. Przyjrzyj się jednej realizacji dla różnych parametrów

	// Obejrzyjmy realizacje zapisana w tablicach u, v, w (zmodyfikuj linie 266, aby wybrac realizacje)
	// Uzupelnij: Rysuj wykres Wykres(...)

	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// Ad. 3. Stworz korelacje w czasie
	double a, b;
	double dt = 0.01;
	double T = 1;

	for(int m = 1; m<nstep; ++m)
	{
		// Uzupelnij: oblicz wspolczynnik a, a = ...
		// Uzupelnij: oblicz wspolczynnik b, b = ...

		for(int k = 0; k<nz-1; ++k)
		{
			// Uzupelnij: wygeneruj skorelowane w czasie pole fluktuacji dla skladowej u ut[k][m] = ...
			// Uzupelnij: wygeneruj skorelowane w czasie pole fluktuacji dla skladowej v vt[k][m] = ...
			// Uzupelnij: wygeneruj skorelowane w czasie pole fluktuacji dla skladowej w wt[k][m] = ...
		}
	}
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	// KONIEC CWICZEN

	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	// Zapisz fluktuacje z kolejnych krokow czasowych do pliku (tylko skladowa u_prim) //////////////////////////////////////
	FILE *g = fopen("Z_korelacja_w_czasie.xls", "w");

	for(int i = 0; i<nz-1; ++i)
	{
		fprintf(g, "%lf\t", z[i]);
		for(int j = 0; j<nstep; ++j)
		{
			fprintf(g, "%lf\t", ut[i][j]);
		}
		fprintf(g, "\n");
	}

	fclose(g);
	// Koniec zapisu /////////////////////////////////////////////////////////////////////////////////////////////////////////

	// Usun tablice
	delete [] x;
	delete [] y;
	delete [] z;
	delete [] fi;
	delete [] teta;
	delete [] psi;
	delete [] wnrf;
	delete [] dkn;
	delete [] kxio;
	delete [] kyio;
	delete [] kzio;
	delete [] sxio;
	delete [] syio;
	delete [] szio;
	delete [] u;
	delete [] v;
	delete [] w;

	// Usun tablice dwuwymiarowe
	for(int i = 0; i<nz-1; ++i)
	{
		delete [] ut[i];
		delete [] vt[i];
		delete [] wt[i];
	}

	delete [] ut;
	delete [] vt;
	delete [] wt;

	return 0;
}