#include "EquationCalculator.h"
#include "SolutionPrinter.h"
#include<math.h>

EquationCalculator::EquationCalculator()
{
	if (econParams.EUfunc != NULL)
		econParams.EUfunc = NULL;

	if (econParams.EVfunc != NULL)
		econParams.EVfunc = NULL;
}

EquationCalculator::~EquationCalculator() { }

void EquationCalculator::SetLoopState(OptimizationParams * oParams, int j, int k, int l, int m)
{
	state.j = j;
	state.k = k;
	state.l = l;
	state.m = m;

	econParams.w1 = calc_w1(oParams);
	econParams.w2 = calc_w2(oParams);

	mexCommand command = oParams->GetCommand();

	if (command == mexCommand::enum_foc_int ||
		command == mexCommand::enum_solve_loop_E1_0 ||
		command == mexCommand::enum_solve_loop_E1_1 ||
		command == mexCommand::enum_solve_loop_E1_0_nk ||
		command == mexCommand::enum_solve_loop_E1_1_nk ||
		command == mexCommand::enum_foc_calc_single)
	{
		if (econParams.EUfunc != NULL)
			delete(econParams.EUfunc);
		econParams.EUfunc = oParams->Get_EUfunc(j, k);

		if (econParams.EVfunc != NULL)
			delete(econParams.EVfunc);
		econParams.EVfunc = oParams->Get_EVfunc(j, k);
	}

	state.level = (command == mexCommand::enum_foc_calc_single) ?
		oParams->Get_print_level() : printLevel::enum_printLevel_unknown;

	if (!econParams.isSet)
		SetEconParams(oParams);
}

void EquationCalculator::SetNormSwitch(int norm_switch)
{
	state.norm_switch = norm_switch;
}

void EquationCalculator::SetEconParams(OptimizationParams * oParams)
{
	econParams.isSet = true;

	int i = oParams->Get_i();
	int tt = oParams->Get_tt();

	int kt_external = oParams->Get_tt();

	econParams.Aprev = oParams->Get_A_grid_val(i);

	//Extract Parameters:
	mexCellItem *MU_0 = oParams->GetStructParamCell("mu_0");
	mexCellItem *MU_l1 = oParams->GetStructParamCell("mu_l1");
	mexCellItem *MU_l2 = oParams->GetStructParamCell("mu_l2");
	mexCellItem *MU_t1 = oParams->GetStructParamCell("mu_t1");
	mexCellItem *MU_t2 = oParams->GetStructParamCell("mu_t2");
	mexCellItem *KT = oParams->GetStructParamCell("KT");
	mexCellItem *CHI = oParams->GetStructParamCell("chi");
	mexCellItem *TAUMU = oParams->GetStructParamCell("taumu");
	mexCellItem *TAU_GAMMA = oParams->GetStructParamCell("tau_gamma");
	mexCellItem *TAU_B = oParams->GetStructParamCell("tau_b");

	econParams.eta_cp = oParams->GetStructParamDouble("eta_cp");
	econParams.r = oParams->GetStructParamDouble("r");
	econParams.beta = oParams->GetStructParamDouble("beta");
	econParams.amin = oParams->GetStructParamDouble("amin");

	econParams.psi_l1 = oParams->GetStructParamDouble("psi_l1");
	econParams.psi_l2 = oParams->GetStructParamDouble("psi_l2");
	econParams.rhou_l = oParams->GetStructParamDouble("rhou_l");

	econParams.psi_t1 = oParams->GetStructParamDouble("psi_t1");
	econParams.psi_t2 = oParams->GetStructParamDouble("psi_t2");
	econParams.rhou_t = oParams->GetStructParamDouble("rhou_t");

	econParams.L_bar = oParams->GetStructParamInt("L_bar");

	int kt = kt_external <=0? KT->getInt(tt): KT->getInt(kt_external);
	kt -= 1; 

	econParams.kttp1 = (kt_external <= 0 ? KT->getInt(tt+1) : KT->getInt(kt_external+1));
	econParams.kttp1 -= 1; 

	econParams.mu_0 = MU_0->getDouble(kt);
	econParams.mu_l1 = MU_l1->getDouble(kt);
	econParams.mu_l2 = MU_l2->getDouble(kt);
	econParams.mu_t1 = MU_t1->getDouble(kt);
	econParams.mu_t2 = MU_t2->getDouble(kt);
	econParams.chi = CHI->getDouble(kt);
	econParams.taumu = TAUMU->getDouble(kt);
	econParams.tau_gamma = TAU_GAMMA->getDouble(kt);
	econParams.tau_b = TAU_B->getDouble(kt);

	if (KT->checkIndexWithinRange(tt + 1))  
	{
		int kt_next = KT->getInt(tt + 1) - 1;
		econParams.mu_0tag = MU_0->getDouble(kt_next);
	}

	econParams.mu_l1_tild = econParams.mu_l1*(1 / econParams.psi_l1 - 1);
	econParams.mu_l2_tild = econParams.mu_l2*(1 / econParams.psi_l2 - 1);

	econParams.mu_t1_tild = econParams.mu_t1*(1 / econParams.psi_t1 - 1);
	econParams.mu_t2_tild = econParams.mu_t2*(1 / econParams.psi_t2 - 1);
}

#pragma region Private Value Calc Region

static double CalcH(EquationCalculator * eqCalc, double L, double T)
{
	return eqCalc->econParams.L_bar - L - T;
}
static double CalcY(EquationCalculator * eqCalc, double H1, double H2)
{
	double result = eqCalc->econParams.tau_b + eqCalc->econParams.w1*H1 + eqCalc->econParams.w2*H2;
	return result <= 0 ? EPS_ZERO : result;
}
static double CalcAtY(EquationCalculator * eqCalc, double Y)
{
	double result = eqCalc->econParams.chi * pow(Y, (1 - eqCalc->econParams.taumu)) - eqCalc->econParams.tau_gamma;
	return result;
}
static double CalcMY(EquationCalculator * eqCalc, double Y)
{
	double result = eqCalc->econParams.chi * (1 - eqCalc->econParams.taumu) * pow(Y, (-eqCalc->econParams.taumu));
	return result;
}
static double CalcA(EquationCalculator * eqCalc, double atY, double C)
{
	double result = (1 + eqCalc->econParams.r)*(eqCalc->econParams.Aprev + atY - C);
	return result;
}


static double CalcC_Generic(EquationCalculator * eqCalc, double mY, double w, double MUL, double mu_l_tild, double L, double psi)
{
	double result = pow((eqCalc->econParams.mu_0 * mY * w / (MUL * mu_l_tild * pow(L, (-1 / psi)))), (eqCalc->econParams.eta_cp));
	return result <= 0 ? EPS_ZERO : result;
}
static double CalcCfromL1(EquationCalculator * eqCalc, double mY, double MUL, double L1)
{
	double w = eqCalc->econParams.w1;
	double mu_l_tild = eqCalc->econParams.mu_l1_tild;
	double psi = eqCalc->econParams.psi_l1;

	return CalcC_Generic(eqCalc, mY, w, MUL, mu_l_tild, L1, psi);
}
static double CalcCfromL2(EquationCalculator * eqCalc, double mY, double MUL, double L2)
{
	double w = eqCalc->econParams.w2;
	double mu_l_tild = eqCalc->econParams.mu_l2_tild;
	double psi = eqCalc->econParams.psi_l2;

	return CalcC_Generic(eqCalc, mY, w, MUL, mu_l_tild, L2, psi);
}
static double CalcCfromA(EquationCalculator * eqCalc, double A, double atY)
{
	return eqCalc->econParams.Aprev + atY - A / (1 + eqCalc->econParams.r);
}

static double CalcVPart_Generic(double rhou_param, double mu_param1, double mu_param2, double psi_param1, double psi_param2, double param1, double param2)
{
	double result = (1 / (1 - rhou_param))*pow(mu_param1*pow(param1, (1 - 1 / psi_param1)) + mu_param2*pow(param2, (1 - 1 / psi_param2)), (1 - rhou_param));
	return result;
}
static double CalcVfromL(EquationCalculator * eqCalc, double L1, double L2)
{
	double rhou_param = eqCalc->econParams.rhou_l;
	double mu_param1 = eqCalc->econParams.mu_l1;
	double mu_param2 = eqCalc->econParams.mu_l2;
	double psi_param1 = eqCalc->econParams.psi_l1;
	double psi_param2 = eqCalc->econParams.psi_l2;

	return CalcVPart_Generic(rhou_param, mu_param1, mu_param2, psi_param1, psi_param2, L1, L2);
}
static double CalcVfromT(EquationCalculator * eqCalc, double T1, double T2)
{
	double rhou_param = eqCalc->econParams.rhou_t;
	double mu_param1 = eqCalc->econParams.mu_t1;
	double mu_param2 = eqCalc->econParams.mu_t2;
	double psi_param1 = eqCalc->econParams.psi_t1;
	double psi_param2 = eqCalc->econParams.psi_t2;

	return CalcVPart_Generic(rhou_param, mu_param1, mu_param2, psi_param1, psi_param2, T1, T2);
}
static double CalcVFromAll(EquationCalculator * eqCalc, double C, double L1, double  L2, double  T1, double  T2, double EV)
{
	return eqCalc->econParams.mu_0 * (1 / (1 - 1 / eqCalc->econParams.eta_cp)) * pow(C , (1 - 1 / eqCalc->econParams.eta_cp))
			- (1 / (1 - eqCalc->econParams.rhou_l)) * pow((eqCalc->econParams.mu_l1 * pow(L1, (1 - 1 / eqCalc->econParams.psi_l1)) + eqCalc->econParams.mu_l2 * pow(L2,  (1 - 1 / eqCalc->econParams.psi_l2))) , (1 - eqCalc->econParams.rhou_l))
			- (1 / (1 - eqCalc->econParams.rhou_t)) * pow((eqCalc->econParams.mu_t1 * pow(T1, (1 - 1 / eqCalc->econParams.psi_t1)) + eqCalc->econParams.mu_t2 * pow(T2,  (1 - 1 / eqCalc->econParams.psi_t2))) , (1 - eqCalc->econParams.rhou_t))
			+ eqCalc->econParams.beta * EV;
}

static double CalcVFromAll_nk(EquationCalculator * eqCalc, double C, double L1, double  L2, double EV)
{
	return eqCalc->econParams.mu_0 * (1 / (1 - 1 / eqCalc->econParams.eta_cp)) * pow(C, (1 - 1 / eqCalc->econParams.eta_cp))
		- (1 / (1 - eqCalc->econParams.rhou_l)) * pow((eqCalc->econParams.mu_l1 * pow(L1, (1 - 1 / eqCalc->econParams.psi_l1)) + eqCalc->econParams.mu_l2 * pow(L2, (1 - 1 / eqCalc->econParams.psi_l2))), (1 - eqCalc->econParams.rhou_l))
		+ eqCalc->econParams.beta * EV;
}

static double CalcParam2_Generic(EquationCalculator * eqCalc, double param, double mu_param1_tild, double mu_param2_tild, double psi_param1, double psi_param2)
{
	double result = pow((eqCalc->econParams.w1 / eqCalc->econParams.w2)*(mu_param2_tild / mu_param1_tild)* pow(param, (1 / psi_param1)), psi_param2);
	return result <= 0 ? EPS_ZERO : result;
}
static double CalcL2fromL1(EquationCalculator * eqCalc, double L1)
{
	double mu_param1_tild = eqCalc->econParams.mu_l1_tild;
	double mu_param2_tild = eqCalc->econParams.mu_l2_tild;
	double psi_param1 = eqCalc->econParams.psi_l1;
	double psi_param2 = eqCalc->econParams.psi_l2;

	return CalcParam2_Generic(eqCalc, L1, mu_param1_tild, mu_param2_tild, psi_param1, psi_param2);
}
static double CalcT2fromT1(EquationCalculator * eqCalc, double T1)
{
	double mu_param1_tild = eqCalc->econParams.mu_t1_tild;
	double mu_param2_tild = eqCalc->econParams.mu_t2_tild;
	double psi_param1 = eqCalc->econParams.psi_t1;
	double psi_param2 = eqCalc->econParams.psi_t2;

	return CalcParam2_Generic(eqCalc, T1, mu_param1_tild, mu_param2_tild, psi_param1, psi_param2);
}
static double CalcT2fromL2(EquationCalculator * eqCalc, double L2)
{
	return eqCalc->econParams.L_bar - L2;
}
static double CalcT1fromL1(EquationCalculator * eqCalc, double L1)
{
	return eqCalc->econParams.L_bar - L1;
}


static double CalcMU_Generic(double param1, double param2, double mu_param1, double mu_param2, double psi_param1, double psi_param2, double rhou_param)
{
	double result = pow(mu_param1*pow(param1, (1 - 1 / psi_param1)) + mu_param2*pow(param2, (1 - 1 / psi_param2)), (-rhou_param));
	return result <= 0 ? EPS_ZERO : result;
}
static double CalcMUL(EquationCalculator * eqCalc, double L1, double L2)
{
	double mu_param1 = eqCalc->econParams.mu_l1;
	double mu_param2 = eqCalc->econParams.mu_l2;
	double psi_param1 = eqCalc->econParams.psi_l1;
	double psi_param2 = eqCalc->econParams.psi_l2;
	double rhou_param = eqCalc->econParams.rhou_l;

	return CalcMU_Generic(L1, L2, mu_param1, mu_param2, psi_param1, psi_param2, rhou_param);
}
static double CalcMUT(EquationCalculator * eqCalc, double T1, double T2)
{
	double mu_param1 = eqCalc->econParams.mu_t1;
	double mu_param2 = eqCalc->econParams.mu_t2;
	double psi_param1 = eqCalc->econParams.psi_t1;
	double psi_param2 = eqCalc->econParams.psi_t2;
	double rhou_param = eqCalc->econParams.rhou_t;

	return CalcMU_Generic(T1, T2, mu_param1, mu_param2, psi_param1, psi_param2, rhou_param);
}
static double CalcF_Generic(double coeff, double val, double neg_inv_exp)
{
	return coeff * pow(val, (-1 / neg_inv_exp));
}
static double CalcFxfromT1(EquationCalculator * eqCalc, double MUT, double T1)
{
	return CalcF_Generic(MUT * eqCalc->econParams.mu_t1_tild, T1, eqCalc->econParams.psi_t1);
}
static double CalcFxfromT2(EquationCalculator * eqCalc, double MUT, double T2)
{
	return CalcF_Generic(MUT * eqCalc->econParams.mu_t2_tild, T2, eqCalc->econParams.psi_t2);
}
static double CalcFxfromL1(EquationCalculator * eqCalc, double MUL, double L1)
{
	return CalcF_Generic(MUL * eqCalc->econParams.mu_l1_tild, L1, eqCalc->econParams.psi_l1);
}
static double CalcFxfromL2(EquationCalculator * eqCalc, double MUL, double L2)
{
	return CalcF_Generic(MUL * eqCalc->econParams.mu_l2_tild, L2, eqCalc->econParams.psi_l2);
}
static double CalcFxfromC(EquationCalculator * eqCalc, double C)
{
	return CalcF_Generic(eqCalc->econParams.mu_0, C, eqCalc->econParams.eta_cp);
}

static double CalcLamParam_Generic(double tild_nom, double tild_denum, double w_nom, double w_denum, double lam, double psi_exp, double param, double psi_nom, double psi_denum)
{
	double result = pow(((tild_nom / tild_denum) * (w_nom / w_denum + lam)), (psi_exp)) * pow(param, (psi_nom / psi_denum));
	return result <= 0 ? EPS_ZERO : result;
}
static double CalcL1fromL2Lam1(EquationCalculator * eqCalc, double L2, double Lam1)
{
	double tild_nom = eqCalc->econParams.mu_l2_tild;
	double tild_denum = eqCalc->econParams.mu_l1_tild;
	double w_nom = eqCalc->econParams.w1;
	double w_denum = eqCalc->econParams.w2;
	double psi_exp = (-eqCalc->econParams.psi_l1); 
	double psi_nom = eqCalc->econParams.psi_l1;
	double psi_denum = eqCalc->econParams.psi_l2;

	return CalcLamParam_Generic(tild_nom, tild_denum, w_nom, w_denum, Lam1, psi_exp, L2, psi_nom, psi_denum);
}
static double CalcT2fromT1Lam1(EquationCalculator * eqCalc, double T1, double lam1)
{
	double tild_nom = eqCalc->econParams.mu_t2_tild;
	double tild_denum = eqCalc->econParams.mu_t1_tild;
	double w_nom = eqCalc->econParams.w1;
	double w_denum = eqCalc->econParams.w2;
	double psi_exp = (eqCalc->econParams.psi_t2);
	double psi_nom = eqCalc->econParams.psi_t2;
	double psi_denum = eqCalc->econParams.psi_t1;

	return CalcLamParam_Generic(tild_nom, tild_denum, w_nom, w_denum, lam1, psi_exp, T1, psi_nom, psi_denum);
}
static double CalcL2fromL1Lam2(EquationCalculator * eqCalc, double L1, double Lam2)
{
	double tild_nom = eqCalc->econParams.mu_l1_tild;
	double tild_denum = eqCalc->econParams.mu_l2_tild;
	double w_nom = eqCalc->econParams.w2;
	double w_denum = eqCalc->econParams.w1;
	double psi_exp = (-eqCalc->econParams.psi_l2); 
	double psi_nom = eqCalc->econParams.psi_l2;
	double psi_denum = eqCalc->econParams.psi_l1;

	return CalcLamParam_Generic(tild_nom, tild_denum, w_nom, w_denum, Lam2, psi_exp, L1, psi_nom, psi_denum);
}
static double CalcT1fromT2Lam2(EquationCalculator * eqCalc, double T2, double Lam2)
{
	double tild_nom = eqCalc->econParams.mu_t1_tild;
	double tild_denum = eqCalc->econParams.mu_t2_tild;
	double w_nom = eqCalc->econParams.w2;
	double w_denum = eqCalc->econParams.w1;
	double psi_exp = (eqCalc->econParams.psi_t1);
	double psi_nom = eqCalc->econParams.psi_t1;
	double psi_denum = eqCalc->econParams.psi_t2;

	return CalcLamParam_Generic(tild_nom, tild_denum, w_nom, w_denum, Lam2, psi_exp, T2, psi_nom, psi_denum);
}

#pragma endregion

#pragma region Calc Interface

int EquationCalculator::GetKttp()
{
	return econParams.kttp1;
}

double EquationCalculator::calc_w1(OptimizationParams * oParams)
{
	int tt = oParams->Get_tt();
	double lw_grid_1_val = oParams->Get_lw_grid_1_val(state.j, tt);
	double u1_grid_val = oParams->Get_u1_grid_val(state.l, tt);
	double expVal = lw_grid_1_val + u1_grid_val;
	return exp(expVal);
}
double EquationCalculator::calc_w2(OptimizationParams * oParams)
{
	int tt = oParams->Get_tt();
	double lw_grid_2_val = oParams->Get_lw_grid_2_val(state.k, tt);
	double u2_grid_val = oParams->Get_u2_grid_val(state.m, tt);
	double expVal = lw_grid_2_val + u2_grid_val;
	return exp(expVal);
}

double EquationCalculator::calc_L2(double L1)
{
	return CalcL2fromL1(this, L1);
}

double EquationCalculator::calc_T2(double T1)
{
	return CalcT2fromT1(this, T1);
}

DoublePair EquationCalculator::calc_A(double L1, double T1, double L2, double T2)
{
	if (L1 <= 0) L1 = EPS_ZERO;
	if (L2 <= 0) L2 = EPS_ZERO;
	if (T1 <= 0) T1 = EPS_ZERO;
	if (T2 <= 0) T2 = EPS_ZERO;

	double H1 = CalcH(this, L1, T1);
	double H2 = CalcH(this, L2, T2);
	double Y = CalcY(this, H1, H2);
	double atY = CalcAtY(this, Y);
	double mY = CalcMY(this, Y);
	double MUL = CalcMUL(this, L1, L2);
	double C = CalcCfromL1(this, mY, MUL, L1);
	double A = CalcA(this, atY, C);

	return { this, A , C };
}

double EquationCalculator::calc_V(double L1, double T1, double L2, double T2, double Csol, double Asol)
{
	double EV = econParams.EVfunc->Evaluate(Asol);
	return CalcVFromAll(this, Csol, L1, L2, T1, T2, EV);
}

double EquationCalculator::calc_U_from_C(double C)
{
	return econParams.mu_0 * pow(C, (-1 / econParams.eta_cp));
}

#pragma endregion

#pragma region Full Solution Calculation Region

SolutionSet EquationCalculator::CalcFullSolutionSetStd(DoublePair LT1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::std_cor;

	result.L1 = LT1.x0;
	result.T1 = LT1.x1;
	result.L2 = CalcL2fromL1(this, result.L1);
	result.T2 = CalcT2fromT1(this, result.T1);
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);
	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.MUL = CalcMUL(this, result.L1, result.L2);
	result.C = CalcCfromL1(this, result.mY, result.MUL, result.L1);
	result.A = CalcA(this, result.atY, result.C);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA(DoublePair LT1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A;

	result.L1 = LT1.x0;
	result.T1 = LT1.x1;

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (LT1.x0 < 0 || LT1.x1 < 0)
	{
		result.L1 = EPS_ZERO;
		result.T1 = EPS_ZERO;
		result.V = -DBL_MAX;
		return result;
	}

	result.L2 = CalcL2fromL1(this, result.L1);
	result.T2 = CalcT2fromT1(this, result.T1);
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);
	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin; 
	result.C = CalcCfromA(this, result.A, result.atY);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 2);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetH2(DoublePair L1LAM2)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_h2;

	result.L1 = L1LAM2.x0;
	result.LAM2 = L1LAM2.x1;
	result.L2 = CalcL2fromL1Lam2(this, result.L1, result.LAM2);
	result.T2 = CalcT2fromL2(this, result.L2);

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (result.T2 <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.T1 = CalcT1fromT2Lam2(this, result.T2, result.LAM2);
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	PrintSolutionPart(&result, &econParams, state.level, 2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.MUL = CalcMUL(this, result.L1, result.L2);
	result.C = CalcCfromL1(this, result.mY, result.MUL, result.L1);
	result.A = CalcA(this, result.atY , result.C);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetH1(DoublePair L2LAM1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_h1;

	result.L2 = L2LAM1.x0;
	result.LAM1 = L2LAM1.x1;
	result.L1 = CalcL1fromL2Lam1(this, result.L2, result.LAM1);
	result.T1 = CalcT1fromL1(this, result.L2);

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (result.T1 <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.T2 = CalcT2fromT1Lam1(this, result.T1, result.LAM1);
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	PrintSolutionPart(&result, &econParams, state.level, 2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.MUL = CalcMUL(this, result.L1, result.L2);
	result.C = CalcCfromL2(this, result.mY, result.MUL, result.L2);
	result.A = CalcA(this, result.atY , result.C);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_H2(DoublePair L1LAM2)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A_h2;

	result.L1 = L1LAM2.x0;
	result.LAM2 = L1LAM2.x1;
	result.L2 = CalcL2fromL1Lam2(this, result.L1, result.LAM2);
	result.T2 = CalcT2fromL2(this, result.L2);

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (result.T2 <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.T1 = CalcT1fromT2Lam2(this, result.T2, result.LAM2);
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	PrintSolutionPart(&result, &econParams, state.level, 2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin;
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_H1(DoublePair L2LAM1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A_h1;

	result.L2 = L2LAM1.x0;
	result.LAM1 = L2LAM1.x1;
	result.L1 = CalcL1fromL2Lam1(this, result.L2, result.LAM1);
	result.T1 = CalcT1fromL1(this, result.L2);

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (result.T1 <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.T2 = CalcT2fromT1Lam1(this, result.T1, result.LAM1);
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	PrintSolutionPart(&result, &econParams, state.level, 2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin; 
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetH1_H2(DoubleTrio L1L2A)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_h1_h2;

	result.L1 = L1L2A.x0;
	result.L2 = L1L2A.x1;
	result.A = L1L2A.x2;
	result.T1 = CalcT1fromL1(this, result.L1);
	result.T2 = CalcT2fromL2(this, result.L2);

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (result.T1 <= 0 || result.T2 <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.H1 = 0;
	result.H2 = 0;
	result.Y = CalcY(this, result.H1, result.H2);

	PrintSolutionPart(&result, &econParams, state.level, 2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	if (result.C <= 0 || (result.A < econParams.amin))
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 4);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_H1_H2(DoublePair L1L2)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A_h1_h2;

	result.L1 = L1L2.x0;
	result.L2 = L1L2.x1;
	result.T1 = CalcT1fromL1(this, result.L1);
	result.T2 = CalcT2fromL2(this, result.L2);

	PrintSolutionPart(&result, &econParams, state.level, 1);

	if (result.T1 <= 0 || result.T2 <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.H1 = 0;
	result.H2 = 0;
	result.Y = CalcY(this, result.H1, result.H2); 

	PrintSolutionPart(&result, &econParams, state.level, 2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin;
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);
	result.V = CalcVFromAll(this, result.C, result.L1, result.L2, result.T1, result.T2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetStd_nk(DoubleSingle LT1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::std_cor;

	result.L1 = LT1.x0;
	result.T1 = 0;
	result.L2 = CalcL2fromL1(this, result.L1);
	result.T2 = 0;
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);
	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.MUL = CalcMUL(this, result.L1, result.L2);
	result.C = CalcCfromL1(this, result.mY, result.MUL, result.L1);
	result.A = CalcA(this, result.atY, result.C);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_nk(DoubleSingle LT1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A;

	result.L1 = LT1.x0;
	result.T1 = 0;

	if (LT1.x0 < 0 )
	{
		result.L1 = EPS_ZERO;
		result.V = -DBL_MAX;
		return result;
	}

	result.L2 = CalcL2fromL1(this, result.L1);
	result.T2 = 0;
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);
	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin; 
	result.C = CalcCfromA(this, result.A, result.atY);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetH2_nk(DoubleSingle L1LAM2)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_h2;

	result.L1 = L1LAM2.x0;
	result.L2 = econParams.L_bar;
	result.T2 = 0;

	result.T1 = 0;
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.MUL = CalcMUL(this, result.L1, result.L2);
	result.C = CalcCfromL1(this, result.mY, result.MUL, result.L1);
	result.A = CalcA(this, result.atY, result.C);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetH1_nk(DoubleSingle L2LAM1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_h1;

	result.L2 = L2LAM1.x0;
	result.L1 = econParams.L_bar;
	result.T1 = 0;

	result.T2 = 0;
	result.H1 = 0;
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.MUL = CalcMUL(this, result.L1, result.L2);
	result.C = CalcCfromL2(this, result.mY, result.MUL, result.L2);
	result.A = CalcA(this, result.atY, result.C);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_H2_nk(DoubleSingle L1LAM2)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A_h2;

	result.L1 = L1LAM2.x0;
	result.L2 = econParams.L_bar;
	result.T2 = 0;

	result.T1 = 0;
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin;
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_H1_nk(DoubleSingle L2LAM1)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A_h1;

	result.L2 = L2LAM1.x0;
	result.L1 = econParams.L_bar;
	result.T1 = 0;

	result.T2 = 0;
	result.H1 = CalcH(this, result.L1, result.T1);
	result.H2 = CalcH(this, result.L2, result.T2);
	result.Y = CalcY(this, result.H1, result.H2);

	if (result.Y <= 0)
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin; 
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 3);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetH1_H2_nk(DoubleSingle L1L2A)
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_h1_h2;

	result.A = L1L2A.x0;
	result.T1 = 0;
	result.T2 = 0;

	result.H1 = 0;
	result.H2 = 0;
	result.Y = CalcY(this, result.H1, result.H2);

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);

	if (result.C <= 0 || (result.A < econParams.amin))
	{
		result.V = -DBL_MAX;
		return result;
	}

	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	PrintSolutionPart(&result, &econParams, state.level, 4);

	return result;
}

SolutionSet EquationCalculator::CalcFullSolutionSetA_H1_H2_nk()
{
	SolutionSet result;
	result.CorIndex = sol_corner::cor_A_h1_h2;

	result.L1 = econParams.L_bar;
	result.L2 = econParams.L_bar;
	result.T1 = 0;
	result.T2 = 0;

	result.H1 = 0;
	result.H2 = 0;
	result.Y = CalcY(this, result.H1, result.H2);

	result.atY = CalcAtY(this, result.Y);
	result.mY = CalcMY(this, result.Y);
	result.A = econParams.amin;
	result.EU = econParams.EUfunc->Evaluate(result.A);
	result.EV = econParams.EVfunc->Evaluate(result.A);
	result.C = CalcCfromA(this, result.A, result.atY);
	result.V = CalcVFromAll_nk(this, result.C, result.L1, result.L2, result.EV);

	return result;
}

#pragma endregion

#pragma region Foc_Int Region

double EquationCalculator::foc_int_terminal_nk(DoubleSingle LT1)
{
	EquationCalculator * eqCalc = LT1.eqCalc;
	double L1 = LT1.x0;
	double L2 = CalcL2fromL1(eqCalc, L1);
	double H1 = CalcH(eqCalc, L1, 0); 
	double H2 = CalcH(eqCalc, L2, 0); 

	double Y = eqCalc->econParams.w1 * H1 + eqCalc->econParams.w2*H2;

	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double C = atY + eqCalc->econParams.Aprev;

	if (C < 0) return 10 ^ 6;

	double MUL = CalcMUL(eqCalc, L1, L2);

	double c_part = mY * eqCalc->econParams.w1 * CalcFxfromC(eqCalc, C);
	double l_part = CalcFxfromL1(eqCalc, MUL, L1);

	double F = (l_part - c_part) / (c_part);
	
	return F;
}

DoublePair EquationCalculator::foc_int(DoublePair LT)
{
	EquationCalculator * eqCalc = LT.eqCalc;
	DoublePair F = { eqCalc, 0, 0 };

	double L1 = LT.x0;
	double T1 = LT.x1;

	double L2 = CalcL2fromL1(eqCalc, L1);
	double T2 = CalcT2fromT1(eqCalc, T1);
	double H1 = CalcH(eqCalc, L1, T1);
	double H2 = CalcH(eqCalc, L2, T2);
	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);
	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);
	
	double C = CalcCfromL1(eqCalc, mY, MUL, L1);
	double A = CalcA(eqCalc, atY, C);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C); 
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;
	
	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? abs(EU_part) : abs(c_part);

	double l_part = CalcFxfromL1(eqCalc, MUL, L1);
	double t_part = CalcFxfromT1(eqCalc, MUT, T1);

	double Fx1_num = (l_part - t_part);
	double Fx1_denom = eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 ? t_part : l_part;

	F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoublePair EquationCalculator::foc_int_std(DoublePair LT1)
{
	EquationCalculator * eqCalc = LT1.eqCalc;
	double L1 = LT1.x0;
	double T1 = LT1.x1;

	double L2 = CalcL2fromL1(eqCalc, L1);
	double T2 = CalcT2fromT1(eqCalc, T1);

	double H1 = CalcH(eqCalc, L1, T1);
	double H2 = CalcH(eqCalc, L2, T2);

	double Y = CalcY(eqCalc, H1, H2);

	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double C = CalcCfromL1(eqCalc, mY, MUL, L1);
	double A = CalcA(eqCalc, atY, C);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C);
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? EU_part : c_part;

	double l_part = CalcFxfromL1(eqCalc, MUL, L1);
	double t_part = CalcFxfromT1(eqCalc, MUT, T1);

	double Fx1_num = (l_part - t_part);
	double Fx1_denom = eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 ? t_part : l_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoublePair EquationCalculator::foc_cor_h1(DoublePair L2LAM1)
{
	EquationCalculator * eqCalc = L2LAM1.eqCalc;
	double L2 = L2LAM1.x0;
	double LAM1 = L2LAM1.x1;

	double L1 = CalcL1fromL2Lam1(eqCalc, L2, LAM1);
	
	double T1 = eqCalc->econParams.L_bar - L1;
	if(T1 <= 0) T1 = pow(10 , -9);

	double T2 = CalcT2fromT1Lam1(eqCalc, T1, LAM1);

	double H2 = CalcH(eqCalc, L2, T2);
	double H1 = 0;

	double Y = CalcY(eqCalc, H1, H2); 

	double atY = CalcAtY(eqCalc, Y); 
	double mY = CalcMY(eqCalc, Y); 


	double MUL = CalcMUL(eqCalc, L1, L2); 
	double MUT = CalcMUT(eqCalc, T1, T2); 

	double C = CalcCfromL2(eqCalc, mY, MUL, L2); 

	double A = CalcA(eqCalc, atY, C); 

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C); 
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? EU_part : c_part;

	double l_part = CalcFxfromL2(eqCalc, MUL, L2);
	double t_part = CalcFxfromT2(eqCalc, MUT, T2);

	double Fx1_num = (l_part - t_part);
	double Fx1_denom = eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 ? t_part : l_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoublePair EquationCalculator::foc_cor_h2(DoublePair L1LAM2)
{
	EquationCalculator * eqCalc = L1LAM2.eqCalc;
	double L1 = L1LAM2.x0;
	double LAM2 = L1LAM2.x1;

	double L2 = CalcL2fromL1Lam2(eqCalc, L1, LAM2);

	double 	T2 = eqCalc->econParams.L_bar - L2;
	if (T2 <= 0) T2 = EPS_ZERO;

	double T1 = CalcT1fromT2Lam2(eqCalc, T2, LAM2);

	double H1 = CalcH(eqCalc, L1, T1);
	double H2 = 0;

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double C = CalcCfromL1(eqCalc, mY, MUL, L1);
	double A = CalcA(eqCalc, atY, C);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C); 
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? EU_part : c_part;

	double l_part = CalcFxfromL1(eqCalc, MUL, L1);
	double t_part = CalcFxfromT1(eqCalc, MUT, T1);

	double Fx1_num = (l_part - t_part);
	double Fx1_denom = eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 ? t_part : l_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoubleTrio EquationCalculator::foc_cor_h1_h2(DoubleTrio L1L2A)
{
	EquationCalculator * eqCalc = L1L2A.eqCalc;
	double L1 = L1L2A.x0;
	double L2 = L1L2A.x1;
	double A  = L1L2A.x2;

	double T1 = eqCalc->econParams.L_bar - L1;
	if (T1 <= 0) T1 = EPS_ZERO;

	double T2 = eqCalc->econParams.L_bar - L2;
	if (T2 <= 0) T2 = EPS_ZERO;

	double Y = eqCalc->econParams.tau_b;
	if (Y <= 0) T1 = EPS_ZERO;

	double atY = CalcAtY(eqCalc, Y);

	double C = atY + eqCalc->econParams.Aprev - A / (1 + eqCalc->econParams.r);
	if (C <= 0) C = EPS_ZERO;

	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double l1_part = CalcFxfromL1(eqCalc, MUL, L1);
	double l2_part = CalcFxfromL2(eqCalc, MUL, L2);

	double t1_part = CalcFxfromT1(eqCalc, MUT, T1);
	double t2_part = CalcFxfromT2(eqCalc, MUT, T2);

	double c_part = CalcFxfromC(eqCalc, C);
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = l1_part - t1_part;
	double Fx0_denom = (eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 || eqCalc->state.norm_switch == 13 || eqCalc->state.norm_switch == 123) ? t1_part : l1_part;
	
	double Fx1_num = l2_part - t2_part;
	double Fx1_denom = (eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 || eqCalc->state.norm_switch == 23 || eqCalc->state.norm_switch == 123) ? t2_part : l2_part;
	
	double Fx2_num = c_part - EU_part;
	double Fx2_denom = (eqCalc->state.norm_switch == 3 || eqCalc->state.norm_switch == 13 || eqCalc->state.norm_switch == 23 || eqCalc->state.norm_switch == 123) ? EU_part : c_part;

	DoubleTrio F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom, Fx2_num / Fx2_denom };

	return F;
}

DoublePair EquationCalculator::foc_cor_A(DoublePair LT)
{
	EquationCalculator * eqCalc = LT.eqCalc;
	double L1 = LT.x0;
	double T1 = LT.x1;

	double L2 = CalcL2fromL1(eqCalc, L1);

	double T2 = CalcT2fromT1(eqCalc, T1);

	double H1 = CalcH(eqCalc, L1, T1);
	double H2 = CalcH(eqCalc, L2, T2);

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);
	double k = (1,1);
	double C = CalcCfromA(eqCalc, eqCalc->econParams.amin, atY);
	if (C <= 0) C = EPS_ZERO;

	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double l_part = CalcFxfromL1(eqCalc, MUL, L1);
	double c_part = mY * eqCalc->econParams.w1 * CalcFxfromC(eqCalc, C);
	double t_part = CalcFxfromT1(eqCalc, MUT, T1);

	double Fx0_num = l_part - c_part;
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? l_part : c_part;

	double Fx1_num = l_part - t_part;
	double Fx1_denom = eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 ? l_part : t_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoublePair EquationCalculator::foc_cor_A_h1(DoublePair L2LAM1)
{
	EquationCalculator * eqCalc = L2LAM1.eqCalc;
	double L2 = L2LAM1.x0;
	double LAM1 = L2LAM1.x1;

	double L1 = CalcL1fromL2Lam1(eqCalc, L2, LAM1);

	double T1 = eqCalc->econParams.L_bar - L1;
	if (T1 <= 0) T1 = EPS_ZERO;

	double T2 = CalcT2fromT1Lam1(eqCalc, T1, LAM1);

	double H1 = 0;
	double H2 = CalcH(eqCalc, L2, T2);

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double C = CalcCfromA(eqCalc, eqCalc->econParams.amin, atY);
	if (C <= 0) C = EPS_ZERO;

	double l1_part = CalcFxfromL1(eqCalc, MUL, L1);
	double l2_part = CalcFxfromL2(eqCalc, MUL, L2);
	double c_part = mY * eqCalc->econParams.w2 * CalcFxfromC(eqCalc, C);
	double t_part = CalcFxfromT1(eqCalc, MUT, T1);

	double Fx0_num = l2_part - c_part;
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? l2_part : c_part;

	double Fx1_num = l1_part - t_part;
	double Fx1_denom = eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12 ? t_part : l1_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoublePair EquationCalculator::foc_cor_A_h2(DoublePair L1LAM2)
{
	EquationCalculator * eqCalc = L1LAM2.eqCalc;

	double L1 = L1LAM2.x0;
	double LAM2 = L1LAM2.x1;

	double L2 = CalcL2fromL1Lam2(eqCalc, L1, LAM2);
	
	double T2 = eqCalc->econParams.L_bar - L2;
	if (T2 <= 0) T2 = EPS_ZERO;

	double T1 = CalcT1fromT2Lam2(eqCalc, T2, LAM2);

	double H1 = CalcH(eqCalc, L1, T1);
	double H2 = 0;

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double C = CalcCfromA(eqCalc, eqCalc->econParams.amin, atY);
	if (C <= 0) C = EPS_ZERO;

	double c_part = mY * eqCalc->econParams.w1 * CalcFxfromC(eqCalc, C);
	double l1_part = CalcFxfromL1(eqCalc, MUL, L1);
	double t_part = CalcFxfromT1(eqCalc, MUT, T1);

	double Fx0_num = l1_part - c_part;
	double Fx0_denom = (eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12) ? l1_part : c_part;

	double Fx1_num = l1_part - t_part;
	double Fx1_denom = (eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12) ? t_part : l1_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoublePair EquationCalculator::foc_cor_A_h1_h2(DoublePair L1L2)
{
	EquationCalculator * eqCalc = L1L2.eqCalc;
	double L1 = L1L2.x0;
	double L2 = L1L2.x1;

	double T1 = eqCalc->econParams.L_bar - L1;
	if (T1 <= 0) T1 = EPS_ZERO;

	double T2 = eqCalc->econParams.L_bar - L2;
	if (T2 <= 0) T2 = EPS_ZERO;

	double H2 = 0;
	double H1 = 0;


	double MUL = CalcMUL(eqCalc, L1, L2);
	double MUT = CalcMUT(eqCalc, T1, T2);

	double l1_part = CalcFxfromL1(eqCalc, MUL, L1);
	double l2_part = CalcFxfromL2(eqCalc, MUL, L2);
	double t1_part = CalcFxfromT1(eqCalc, MUT, T1);
	double t2_part = CalcFxfromT2(eqCalc, MUT, T2);

	double Fx0_num = l1_part - t1_part;
	double Fx0_denom = (eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12) ? t1_part : l1_part;

	double Fx1_num = l2_part - t2_part;
	double Fx1_denom = (eqCalc->state.norm_switch == 2 || eqCalc->state.norm_switch == 12) ? t2_part : l2_part;

	DoublePair F = { eqCalc, Fx0_num / Fx0_denom , Fx1_num / Fx1_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_int_std_nk(DoubleSingle LT1)
{
	EquationCalculator * eqCalc = LT1.eqCalc;
	double L1 = LT1.x0;

	double L2 = CalcL2fromL1(eqCalc, L1);

	double H1 = CalcH(eqCalc, L1, 0);
	double H2 = CalcH(eqCalc, L2, 0);

	double Y = CalcY(eqCalc, H1, H2);
	
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);

	double C = CalcCfromL1(eqCalc, mY, MUL, L1);
	double A = CalcA(eqCalc, atY, C);
	
	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C);
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? EU_part : c_part;

	DoubleSingle F = { eqCalc, Fx0_num / Fx0_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_cor_A_nk(DoubleSingle LT)
{
	EquationCalculator * eqCalc = LT.eqCalc;
	double L1 = LT.x0;

	double L2 = CalcL2fromL1(eqCalc, L1);

	double H1 = CalcH(eqCalc, L1, 0);
	double H2 = CalcH(eqCalc, L2, 0);

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);
	double k = (1, 1);
	double C = CalcCfromA(eqCalc, eqCalc->econParams.amin, atY);
	if (C <= 0) C = EPS_ZERO;

	double MUL = CalcMUL(eqCalc, L1, L2);

	double l_part = CalcFxfromL1(eqCalc, MUL, L1);
	double c_part = mY * eqCalc->econParams.w1 * CalcFxfromC(eqCalc, C);

	double Fx0_num = l_part - c_part;
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? l_part : c_part;

	DoubleSingle F = { eqCalc,  Fx0_num / Fx0_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_cor_h2_nk(DoubleSingle L1LAM2)
{
	EquationCalculator * eqCalc = L1LAM2.eqCalc;
	double L1 = L1LAM2.x0;

	double L2 = eqCalc->econParams.L_bar;

	double H1 = CalcH(eqCalc, L1, 0);
	double H2 = 0;

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);

	double C = CalcCfromL1(eqCalc, mY, MUL, L1);
	double A = CalcA(eqCalc, atY, C);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C);

	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? EU_part : c_part;

	DoubleSingle F = { eqCalc, Fx0_num / Fx0_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_cor_h1_nk(DoubleSingle L2LAM1)
{
	EquationCalculator * eqCalc = L2LAM1.eqCalc;
	double L2 = L2LAM1.x0;

	double L1 = eqCalc->econParams.L_bar;

	double H2 = CalcH(eqCalc, L2, 0);
	double H1 = 0;

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);
	double C = CalcCfromL2(eqCalc, mY, MUL, L2);

	double A = CalcA(eqCalc, atY, C);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C);
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = (c_part - EU_part);
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? EU_part : c_part;

	DoubleSingle F = { eqCalc, Fx0_num / Fx0_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_cor_A_h2_nk(DoubleSingle L1LAM2)
{
	EquationCalculator * eqCalc = L1LAM2.eqCalc;

	double L1 = L1LAM2.x0;

	double L2 = eqCalc->econParams.L_bar;

	double H1 = CalcH(eqCalc, L1, 0);
	double H2 = 0;

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);

	double C = CalcCfromA(eqCalc, eqCalc->econParams.amin, atY);

	double c_part = mY * eqCalc->econParams.w1 * CalcFxfromC(eqCalc, C);
	double l1_part = CalcFxfromL1(eqCalc, MUL, L1);

	double Fx0_num = l1_part - c_part;
	double Fx0_denom = (eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12) ? l1_part : c_part;

	DoubleSingle F = { eqCalc,  Fx0_num / Fx0_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_cor_A_h1_nk(DoubleSingle L2LAM1)
{
	EquationCalculator * eqCalc = L2LAM1.eqCalc;
	double L2 = L2LAM1.x0;

	double L1 = eqCalc->econParams.L_bar;

	double H1 = 0;
	double H2 = CalcH(eqCalc, L2, 0);

	double Y = CalcY(eqCalc, H1, H2);
	double atY = CalcAtY(eqCalc, Y);
	double mY = CalcMY(eqCalc, Y);

	double MUL = CalcMUL(eqCalc, L1, L2);

	double C = CalcCfromA(eqCalc, eqCalc->econParams.amin, atY);

	double l2_part = CalcFxfromL2(eqCalc, MUL, L2);
	double c_part = mY * eqCalc->econParams.w2 * CalcFxfromC(eqCalc, C);

	double Fx0_num = l2_part - c_part;
	double Fx0_denom = eqCalc->state.norm_switch == 1 || eqCalc->state.norm_switch == 12 ? l2_part : c_part;

	DoubleSingle F = { eqCalc, Fx0_num / Fx0_denom };

	return F;
}

DoubleSingle EquationCalculator::foc_cor_h1_h2_nk(DoubleSingle L1L2A)
{
	EquationCalculator * eqCalc = L1L2A.eqCalc;
	double A = L1L2A.x0;
	double Y = eqCalc->econParams.tau_b;
	double atY = CalcAtY(eqCalc, Y);

	double C = CalcCfromA(eqCalc, A,atY);

	double EU = eqCalc->econParams.EUfunc->Evaluate(A);

	double c_part = CalcFxfromC(eqCalc, C);
	double EU_part = (1 + eqCalc->econParams.r) * eqCalc->econParams.beta * EU;

	double Fx0_num = c_part - EU_part;
	double Fx0_denom = (eqCalc->state.norm_switch == 1) ? EU_part : c_part;

	DoubleSingle F = { eqCalc, Fx0_num / Fx0_denom };

	return F;
}


#pragma endregion

