#include <iostream>
#include <discpp.h>

static Dislin dlin;

#define H_0 70*1000/3.08567758E22
#define GYR 365.25*24.0*60.0*60.0*1E9
#define totalSize sizeBack + sizeForward + 1

const double GYRH_0 = GYR*H_0;
char legendText[255] = { '/0' };
int legendIndex = 1;

void friedmannPlot(double Omega_m, double Omega_v, float nGy, float pGy, char* curveName, char* col = "WHITE",float w_v = -1.0f, double dT = 0.00001){
	int sizeBack = GYRH_0 / dT*fabs(nGy), sizeForward = GYRH_0 / dT*fabs(pGy);
	double *a = new double[totalSize], *t = new double[totalSize];
	memset(a, 0, sizeof(double)*(totalSize));
	memset(t, 0, sizeof(double)*(totalSize));
	a[sizeBack] = 1;
	t[sizeBack] = 0;
	double a_0_min = INFINITY, t_0_min = NAN;
	double a_0_max = -INFINITY, t_0_max = NAN;
	double a_z = NAN, t_z = NAN; // a = 1/2 corresponds to z = 1 because 1+z = a_obs / a_emitted

	std::cout << curveName << "(Omega_m = " << Omega_m << ", Omega_v = " << Omega_v << "):" << std::endl;

	auto helperLoop = [&](int i, int loopDir){
		double a_0 = 1,	t_0 = 0, temp;
		int acceldir = loopDir;
		for (; (i >= 0) && (i < totalSize); i += loopDir){
			temp = (Omega_m / a_0 + Omega_v*pow(a_0,2-3*(1+w_v)) + (1 - Omega_m - Omega_v));
			if (acceldir > 0 && temp < 0){
				acceldir = -1;
			}
			else if (acceldir < 0 && temp < 0){
				acceldir = 1;
			}

			if (a_0 > 0)
				a_0 += dT*acceldir*sqrt(abs(temp));
			t_0 += loopDir*dT;

			if (a_0 <= 0.5 && isnan(t_z)){
				t_z = t_0;
				a_z = a_0;
			}

			if (a_0 < 0){
				a_0 = 0;
				std::cout << "\ta_0 = 0 at: " << t_0 / GYRH_0 << "Gyr" << std::endl;
			}

			if (a_0 < a_0_min){
				a_0_min = a_0;
				t_0_min = t_0;
			}
			if (a_0 > a_0_max){
				a_0_max = a_0;
				t_0_max = t_0;
			}

			a[i] = a_0;
			t[i] = t_0 / GYRH_0;
		}
	};

	helperLoop(sizeBack - 1, -1);
	helperLoop(sizeBack + 1, 1);

	std::cout << "\ta_0_max = " << a_0_max << " at: " << t_0_max / GYRH_0 << "Gyr. \t" << "a_0_min = " << a_0_min << " at: " << t_0_min / GYRH_0 << "Gyr." << std::endl;
	std::cout << "\ta_Z = " << a_z << " at: " << t_z / GYRH_0 << "Gyr." << std::endl;
	
	std::cout << std::endl;

	dlin.color(col);
	dlin.leglin(legendText, curveName, legendIndex++);
	dlin.curve(t, a, totalSize);

	delete[] a;
	delete[] t;
};

int main() {
	dlin.metafl("PNG");
	dlin.filmod("COUNT");
	dlin.scrmod("REVERS");
#ifdef _DEBUG
	dlin.metafl("XWIN");
#endif
	auto dislinInit = [](){
		dlin.disini();
		//dlin.nochek();
		dlin.pagera();
		dlin.name("Time t (Gyr)", "X");
		dlin.name("Scale Factor a(t)", "Y");
		legendIndex = 1;
		dlin.complx();
		
	};
	
	int nGYr, pGyr;

	//a)
	dislinInit();
	nGYr = -15;
	pGyr = 10;
	dlin.graf(nGYr, pGyr, nGYr, 5, 0, 2.78, 0, 0.5);

	dlin.legini(legendText, 5, 50);
	dlin.legtit("Legend");

	friedmannPlot(0.3, 0.7, nGYr, pGyr, "Concordance", "WHITE"); // Concordance
	friedmannPlot(1.0, 0.0, nGYr, pGyr, "Critical Matter", "RED"); // Critical Matter
	friedmannPlot(0.3, 0.0, nGYr, pGyr, "Low Density", "BLUE"); // Low Density
	friedmannPlot(0.0, 0.0, nGYr, pGyr, "Empty", "GREEN"); // Empty
	friedmannPlot(0.1, 3.0, nGYr, pGyr, "Bounce", "CYAN"); // Bounce

	dlin.color("WHITE");
	dlin.legpos(1000, 650);
	dlin.legend(legendText,legendIndex-1);
	

	dlin.disfin();

	//b)
	dislinInit();
	nGYr = -80;
	pGyr = 10;
	dlin.graf(nGYr, pGyr, nGYr, 10, 0, 2.2707, 0, 0.5);

	dlin.legini(legendText, 4, 50);
	dlin.legtit("Legend");

	friedmannPlot(1.0, 0.0, nGYr, pGyr, "Critical Matter", "WHITE"); // Critical Matter
	friedmannPlot(0.0, 0.0, nGYr, pGyr, "Empty", "RED"); // Empty
	friedmannPlot(0.0, 1.0, nGYr, pGyr, "Vacuum Only", "BLUE"); // Vacuum Only
	friedmannPlot(0.3, 1.713, nGYr, pGyr, "Loitering", "GREEN"); // Loitering

	dlin.color("WHITE");
	dlin.legpos(1000, 650);
	dlin.legend(legendText, legendIndex-1);


	dlin.disfin();

	//c)	
	dislinInit();
	nGYr = -10;
	pGyr = 55;
	dlin.graf(nGYr, pGyr, nGYr, 10, 0, 2, 0, 0.5);

	dlin.legini(legendText, 2, 50);
	dlin.legtit("Legend");

	friedmannPlot(3.0, 0.0, nGYr, pGyr, "Crunch", "WHITE"); // Crunch
	friedmannPlot(3.0, 0.1, nGYr, pGyr, "Crunch with Vacuum", "RED"); // Crunch with Vacuum

	dlin.color("WHITE");
	dlin.legpos(1000, 1350);
	dlin.legend(legendText, legendIndex-1);

	dlin.disfin();

	//Bonus
	dislinInit();
	nGYr = -15;
	pGyr = 10;
	dlin.graf(nGYr, pGyr, nGYr, 10, 0, 2, 0, 0.5);

	dlin.legini(legendText, 2, 50);
	dlin.legtit("Legend");

	friedmannPlot(0.3, 0.7, nGYr, pGyr, "Concordance with w_v = -1", "WHITE");
	friedmannPlot(0.3, 0.7, nGYr, pGyr, "Concordance with w_v = -0.9", "RED",-0.9f);

	dlin.color("WHITE");
	dlin.legpos(800, 550);
	dlin.legend(legendText, legendIndex - 1);

	dlin.disfin();

	system("PAUSE");
	return(0);
}


