#pragma once
#include <vector>
#include <iostream>

double norm(double mu, double std, double t);
std::vector <double> generateNormD(double mu, double std, double h);
std::vector <double> generateNormCDF(double mu, double std, double h);
std::vector <double> generateLysisDistr(double mu, double std, double h);
double findMax(std::vector<double> vec);
double findMin(std::vector<double> vec);
double summ(std::vector<double> vec);

extern double h;		//step size
extern double interval;	//in hours
extern int len;

extern double V;			// (volum: L)
extern double c0; // (nutrient inflow concentration g / L)
extern double M; //Bacterial dry mass / cell(g)
extern double y;		// (stochiometric nutrient conversion to bacteria mass(g / L))
extern double Q;			// (flow rate L / h)
extern double q;		//flow through
extern double N0;// (initial nutrient mass g)

				 // Susceptible bacteria
extern double Vmax_1;		// (max growth rate h^. - 1)
extern double S0_1; //initial bacteria number
extern double Km_1; // (half Vmax conc g / mL)
extern double cBac_1;
extern double Q_1;

extern double Vmax_2;		// (max growth rate h^. - 1)
extern double S0_2;//3.33*pow(10, 11); //initial bacteria number
extern double Km_2; // (half Vmax conc g / mL)
extern double cBac_2;
extern double Q_2;

extern double a_1; //adsorption phage - 1 cell - 1 ml - 1 h - 1 NB!need proper units for use wth mass(dry mass conversion)
extern int b_1;				// phage released upon lysis
extern double T_1;			//average lysis time
extern double P0_1;	//initial phage number
extern double std_1;		//Standard deviation of bacterial lysis time
extern std::vector <double> distr_1;

extern double a_2; //adsorption phage - 1 cell - 1 ml - 1 h - 1 NB!need proper units for use wth mass(dry mass conversion)
extern int b_2;				// phage released upon lysis
extern double T_2;			//average lysis time
extern double P0_2;//pow(10, 3);	//initial phage number
extern double std_2;		//Standard deviation of bacterial lysis time
extern std::vector <double> distr_2;


//Definition of a variable class where all variables are held and modified. Makes solver simpler to use.

class Var2Bac2Phage {
public:
	double N, S1, S2, P1, P2;
	std::vector <double> I1;	//S1 infected by P1
	std::vector <double> I2;	//S1 infected by P2
	std::vector <double> I3;	//S2 infected by P2

	Var2Bac2Phage() : N(0), S1(0), S2(0), P1(0), P2(0) {
		for (int i = 0; i < len; ++i) {	
			I1.push_back(0);
			I2.push_back(0);
			I3.push_back(0);
		}
	}

	Var2Bac2Phage(double N, double S1, double S2, double P1, double P2, std::vector <double> &I1, std::vector <double> &I2, std::vector <double> &I3) : N(N), S1(S1), S2(S2), P1(P1), P2(P2), I1(I1), I2(I2), I3(I3) {
	};
};

Var2Bac2Phage operator+(Var2Bac2Phage lhs, const Var2Bac2Phage& rhs);
Var2Bac2Phage operator*(double a, Var2Bac2Phage rhs);
Var2Bac2Phage func2Bac2Phage(Var2Bac2Phage vars);


//Runge kutta 4th order solver. Template for easy changing of variable class. 

template <class T>
T rungeKutta2(T vars) {
	T k1 = func2Bac2Phage(vars);
	T k2 = func2Bac2Phage(vars + (h / 2)*k1);
	T k3 = func2Bac2Phage(vars + (h / 2)*k2);
	T k4 = func2Bac2Phage(vars + h*k3);

	T ny = vars + (h / 6)*(k1 + 2 * k2 + 2 * k3 + k4);
	ny.I1.back() = (h / 6)*(k1.I1.back() + 2 * k2.I1.back() + 2 * k3.I1.back() + k4.I1.back());
	ny.I2.back() = (h / 6)*(k1.I2.back() + 2 * k2.I2.back() + 2 * k3.I2.back() + k4.I2.back());
	ny.I3.back() = (h / 6)*(k1.I3.back() + 2 * k2.I3.back() + 2 * k3.I3.back() + k4.I3.back());
	return ny;
}