#include "ihacres_auto_cal.h"

#include <stdlib.h> // random numbers
#include <time.h>	// random numbers

#include <math.h>	// function "pow()"

#include <fstream>  // used for textfile output (test only)

using std::ofstream;
//---------------------------------------------------------

Cihacres_auto_cal::Cihacres_auto_cal()
{
	//-----------------------------------------------------
	// 1. Info...
	Set_Name("IHACRES Automatic Calibration");

	Set_Author("copyrights (c) 2007 Stefan Liersch");

	Set_Description("Automatic Calibration Tool for the Model IHACRES \n \n \n"
		"Reference: \n \n"
		"Jakeman, A.J. / Hornberger, G.M. (1993). \n"
		"How Much Complexity Is Warranted in a Rainfall-Runoff Model? \n"
		"Water Resources Research, (29), NO. 8 (2637-2649) \n \n"   
		"Kokkonen, T. S. et al. (2003). \n"
		"Predicting daily flows in ungauged catchments:"
		"model regionalization from catchment descriptors"
		"at the Coweeta Hydrologic Laboratory, North Carolina \n "
		"Hydrological Processes (17), 2219-2238 \n \n"
		"Croke, B. F. W., W. S. Merritt, et al. (2004).\n"
		"A dynamic model for predicting hydrologic response"
		"to land cover changes in gauged and"
		"ungauged catchments. \n"
		"Journal Of Hydrology 291(1-2): 115-131."
		);

	//-----------------------------------------------------
	// Create First Module Dialog
	_CreateDialog1();
	//-----------------------------------------------------
}
//---------------------------------------------------------

Cihacres_auto_cal::~Cihacres_auto_cal()
{
	//delete[] T_NLP;
}
//---------------------------------------------------------

bool Cihacres_auto_cal::On_Execute()
{
ofstream file("IHAC_TEST.txt");
ofstream file2("occ.txt");
ofstream error("error.txt");
	int				i,j,k;
	CSG_Table		*pTable; //, *pTable_out, *pTable_parms;
	int*			TopIndex1;
	int*			TopIndex2;
	int*			TopIndex3;
	double*			TopVal1;
	double*			TopVal2;
	double*			TopVal3;
	//int				steps;
	int				storconf;	// storage configuration
								// 0 = single,
								// 1 = two parallel,
								// 2 = two in series
	double			dev_eRainDis; // deviation of estimated effective rainfall
								  // to observed discharge in [%]
	double			dev_eRainPCP; // deviation of estimated effective rainfall
								  // to observed precipitation in [%]
	double			NSE;		// Nash-Sutcliffe Efficiency

	double			l_temp = 0;
	double			p_temp = 1;

	//---------------------------------------------------------
	// Assign parameters from First Module Dialog
	pTable				= Parameters("TABLE")			->asTable();
	// Field numbers
	T_in.dateField		= Parameters("DATE_Field")		->asInt();
	T_in.dischargeField	= Parameters("DISCHARGE_Field")	->asInt();
	T_in.pcpField		= Parameters("PCP_Field")		->asInt();
	T_in.tmpField		= Parameters("TMP_Field")		->asInt();
	T_in.bTMP			= Parameters("USE_TMP")			->asBool();

	area				= Parameters("AREA")			->asDouble();
	dev_eRainDis		= Parameters("DEV_ERAINDIS")	->asDouble();
	dev_eRainPCP		= Parameters("DEV_ERAINPCP")	->asDouble();
	//steps				= Parameters("STEP")			->asInt();
	storconf			= Parameters("STORAGE")			->asInt();
	IHAC_version		= Parameters("IHACVERS")		->asInt();
	T_in.first = T_in.last = 0;
	//---------------------------------------------------------

	//---------------------------------------------------------
	// After pressing OK in the first Dialog Dialog2 appears.
	// The user selects the first and the last date
//	if (_CreateDialog2(bTMP,bnl_only,storconf,pTable,dateField,dischargeField))
	if (_CreateDialog2(storconf, pTable, T_in, IHAC_version))
	{
error << "after Dialog2" << endl;
error << "IHAC_version: " << IHAC_version << endl;
		//---------------------------------------------------------
		// Searching the first and the last record of the time range
		ihac.AssignFirstLastRec(*pTable, T_in.first, T_in.last, T_in.date1, T_in.date2, T_in.dateField);
		T_in.size = T_in.last - T_in.first + 1;
error << "after AssignFirstLlast" << endl;
		//---------------------------------------------------------
		// Initialize arrays
		_InitArrays(IHAC_version, T_in, T_SM);
error << "after _InitArrays" << endl;
		//_InitArrays_nsim(nsim);
		//---------------------------------------------------------
		// Initialize global variables
		_ZeroVariables();
		//---------------------------------------------------------
		// Initialize variables of T_Cal (Calibration Parameters)
		//_Init_T_Cal(steps);
		//-----------------------------------------------------
		// Initialize array p_bValidSim (valid simulations)
error << "T_CP.ncombi " << T_CP.ncombi << endl;
		p_bValidSim = new bool[T_CP.ncombi];
		for (i = 0; i < T_CP.ncombi; i++)
			p_bValidSim[i] = false;
		//-----------------------------------------------------
		// Assign parameter sets for non-linear module

		T_NLP = new T_NonLinear_Parms[T_CP.ncombi];
		//T_NLP = _Set_PreCalibration1_Parms(T_CP); //, T_NLP);
		_Set_PreCalibration1_Parms(T_CP);
error << "after PreCalibration1" << endl;
		//---------------------------------------------------------
		// Read input table
		// Assign selected time range to vector date, 
		// discharge[], pcp[], and tmp[]
		for (j = 0, k = T_in.first; j < T_in.size, k < T_in.last + 1; j++, k++)
		{
			// put the date string into the vector
			T_in.date[j].append(pTable->Get_Record(k)->asString(T_in.dateField));
			T_in.p_discharge[j]	= pTable->Get_Record(k)->asDouble(T_in.dischargeField);
			T_in.p_pcp[j]		= pTable->Get_Record(k)->asDouble(T_in.pcpField);
			T_in.p_tmp[j]		= pTable->Get_Record(k)->asDouble(T_in.tmpField);
		}
		//---------------------------------------------------------
		// Convert streamflow from m3/s to mm/day
		ihac.ConvM3Sec2MMDay(T_in.p_discharge,T_in.p_dischargeMM,area,T_in.size);
		// calculate sum of observed discharge in [mm]
		// this value is used to calibrate the non-linear module parameters
		sum_obsDisMM = ihac.SumVector(T_in.p_dischargeMM, T_in.size);

		/*
		//---------------------------------------------------------
		// SNOW MODULE
		//---------------------------------------------------------
		// Calculate daily MeltRate
		if (IHAC_version == 2)
		{
			for (int i = 0; i < sizeAll; i++)
			{
				T_SM.p_SnowStorage[i] = 0.0;
				T_SM.p_MeltRate[i] = 0.0;
			}
			T_SM.T_Rain = -1.0;
			T_SM.T_Melt = 1.0;
			T_SM.DD_FAC = 3.5;
			ihac.SnowModule(T_SM.p_SnowStorage,T_SM.p_MeltRate,T_in.tmp,T_in.pcp, T_SM.T_Rain,T_SM.T_Melt,T_SM.DD_FAC,T_in.size);
		}
		//---------------------------------------------------------
		*/

		//---------------------------------------------------------
		// Perform non-linear module calibration
		_CalcNonLinearModule(T_NLP,T_in.bTMP, IHAC_version, T_CP.ncombi, dev_eRainDis, dev_eRainPCP, bSnowModule);
		//---------------------------------------------------------
error << "after _CalcNonLinearModule" << endl;
		//---------------------------------------------------------
		// Select valid simulations
		file2 << "nValidSims: " << nValidSims << endl;
		if (nValidSims > 0)
		{
			T_NLP_valid = new T_NonLinear_Parms[nValidSims];
error << "after T_NLP_valid = new" << endl;
			//T_NLP_valid = _Select_Valid_Sims(T_CP.ncombi, T_NLP, T_NLP_valid);
			_Select_Valid_Sims(T_CP.ncombi, T_NLP); 
error << "_Select_Valid_Sims" << endl;
		//---------------------------------------------------------
			//---------------------------------------------------------
			// count number of occurences of parameter values in valid simulations
			T_occ = new T_nOccurence[3]; // 3 Parameters
			// initialize T_occ (assign values to T_occ parameters)
			//T_occ = _Init_Tocc(T_occ, T_CP);
			_Init_Tocc(T_CP);
			// select n_best in order to refine parameter ranges
			TopIndex1 = new int[T_occ[0].n_best];
			TopIndex2 = new int[T_occ[1].n_best];
			TopIndex3 = new int[T_occ[2].n_best];
/*
			TopIndex1 = _FindTopIndex(T_occ[0].n_occ,T_occ[0].nsteps+1,TopIndex1,T_occ[0].n_best,nValidSims);
			TopIndex2 = _FindTopIndex(T_occ[1].n_occ,T_occ[1].nsteps+1,TopIndex2,T_occ[1].n_best,nValidSims);
			TopIndex3 = _FindTopIndex(T_occ[2].n_occ,T_occ[2].nsteps+1,TopIndex3,T_occ[2].n_best,nValidSims);
*/
			TopIndex1 = _FindTopIndex(T_occ[0].n_occ,T_occ[0].nsteps+1,T_occ[0].n_best,nValidSims);
			TopIndex2 = _FindTopIndex(T_occ[1].n_occ,T_occ[1].nsteps+1,T_occ[1].n_best,nValidSims);
			TopIndex3 = _FindTopIndex(T_occ[2].n_occ,T_occ[2].nsteps+1,T_occ[2].n_best,nValidSims);

error << "after TopIndex=" << endl;
error << "TopIndex1= " << TopIndex1 << endl;
error << "TopIndex2= " << TopIndex2 << endl;
error << "TopIndex3= " << TopIndex3 << endl;
error << "TopIndex1= " << &TopIndex1 << endl;
error << "TopIndex2= " << &TopIndex2 << endl;
error << "TopIndex3= " << &TopIndex3 << endl;
			TopVal1 = new double[T_occ[0].n_best];
			TopVal2 = new double[T_occ[1].n_best];
			TopVal3 = new double[T_occ[2].n_best];

			TopVal1 = _AssignTopVal(T_occ[0],TopIndex1);
			TopVal2 = _AssignTopVal(T_occ[1],TopIndex2);
			TopVal3 = _AssignTopVal(T_occ[2],TopIndex3);
error << " after AssignTopVal =" << endl;
			// Assign new parameter bounds/ranges
			T_occ[0].new_lb = _lowest(TopVal1, T_occ[0].n_best);
			T_occ[0].new_ub = _largest(TopVal1, T_occ[0].n_best);
			T_occ[1].new_lb = _lowest(TopVal2, T_occ[1].n_best);
			T_occ[1].new_ub = _largest(TopVal2, T_occ[1].n_best);
			T_occ[2].new_lb = _lowest(TopVal3, T_occ[2].n_best);
			T_occ[2].new_ub = _largest(TopVal3, T_occ[2].n_best);
	
			file2 << "lower bound, upper bound" << endl;
			file2 << "Tw: ," << T_occ[0].new_lb << " , " << T_occ[0].new_ub << " , ";
			file2 << "n_steps: " << T_occ[0].nsteps << " , ";
			file2 << "n_best: " << T_occ[0].n_best << endl;
			file2 << "f : ," <<T_occ[1].new_lb << " , " << T_occ[1].new_ub << " , ";
			file2 << "n_steps: " << T_occ[1].nsteps << " , ";
			file2 << "n_best: " << T_occ[1].n_best << endl;
			file2 << "c : ," <<T_occ[2].new_lb << " , " << T_occ[2].new_ub << " , ";
			file2 << "n_steps: " << T_occ[2].nsteps << " , ";
			file2 << "n_best: " << T_occ[2].n_best << endl;
			
			file2 << "Tw, f, c" << endl;
			for (i = 0; i < T_occ[0].n_best; i++)
			{
				file2 << TopIndex1[i] << ", ";
				if (i < T_occ[1].n_best)
					file2 << TopIndex2[i] << ", ";
				if (i < T_occ[2].n_best)
					file2 << TopIndex3[i] << endl;
			}

			file2 << "Tw,n_Tw,f,n_f,c,n_c" << endl;
			for (i = 0; i < T_occ[0].nsteps+1; i++)
			{
				file2 << T_occ[0].val[i] << " , ";
				file2 << T_occ[0].n_occ[i] << " , ";
				if (i < T_occ[1].nsteps+1)
				{
					file2 << T_occ[1].val[i] << " , ";
					file2 << T_occ[1].n_occ[i] << " , ";
				}
				if (i < T_occ[2].nsteps+1)
				{
					file2 << T_occ[2].val[i] << " , ";
					file2 << T_occ[2].n_occ[i] << endl;
				}
			}

			file << "simulation,Tw,c,f,l,p,sum eRainGT,Dif_eRainDis,T_NLP[j].Sum_eRainMM" << endl;
			for (i = 0; i < nValidSims; i++)
			{
				file << i << ", ";
				file << T_NLP_valid[i].Tw  <<",";
				file << T_NLP_valid[i].c  <<",";
				file << T_NLP_valid[i].f  <<",";
				file << T_NLP_valid[i].l  <<",";
				file << T_NLP_valid[i].p  <<",";
				file << T_NLP_valid[i].sum_eRainGTpcp  <<",";
				file << T_NLP_valid[i].Dif_eRainDis  <<",";
				file << T_NLP_valid[i].Sum_eRainMM << endl;
			}

			delete[] T_NLP_valid;

		} // end if (nValidSims > 0)
		//---------------------------------------------------------

	} // end if (_CreateDialog2)

	delete[] T_NLP;
	delete[] p_bValidSim;


	// **************************************************************
	// Perform Monte Carlo simulation of the non-linear module
	// **************************************************************
	//int nsimMC_NLM = 10000;
	p_bValidSim = new bool[nsimMC_NLM];
	for (i = 0; i < nsimMC_NLM; i++) p_bValidSim[i] = false;
	T_NLP = new T_NonLinear_Parms[nsimMC_NLM];
	_CalcNonLinearModuleMC(T_NLP,T_occ,T_in.bTMP, IHAC_version, nsimMC_NLM, dev_eRainDis, dev_eRainPCP);
error << "after NonLinearModuleMC" << endl;
	file2 << endl;
	file2 << "Monte Carlo Simulation ***********************" << endl;
	file2 << "nValidSims: " << nValidSims << endl;
	file2 << "dev_eRainDis: " << dev_eRainDis << endl;
	file2 << "dev_eRainPCP: " << dev_eRainPCP << endl;
	file2 << "Tw, f, c,.l.p sum_eRainGTpcp, Dif_eRainDis" << endl;
	if (nValidSims > 0)
	{
		T_NLP_valid = new T_NonLinear_Parms[nValidSims];
		//T_NLP_valid = _Select_Valid_Sims(nsim, T_NLP, T_NLP_valid);
		_Select_Valid_Sims(nsimMC_NLM, T_NLP);

		for (i = 0; i < nValidSims; i++)
		{
			file2 << T_NLP_valid[i].Tw << ", " << T_NLP_valid[i].f << ", " << T_NLP_valid[i].c << ", ";
			file2 << T_NLP_valid[i].l << ", "<< T_NLP_valid[i].p << ", ";
			file2 << T_NLP_valid[i].sum_eRainGTpcp << ", " << T_NLP_valid[i].Dif_eRainDis << endl;
		}

		// ------------------------------------------------------------
		// Select best non-linear MC-Simulation parameter settings
		// using Top 10 parameter settings according to the lowest
		// difference of simulated effective rainfall and observed streamflow
		// Be aware that nTop can be reduced to the number of valid simulations!
		nTop = 10;
		Top_NLM = _Assign_TopNLM(nTop, T_NLP_valid);

		file2 << endl << endl;
		file2 << "Top Simulations" << endl;
		file2 << "Tw, f, c,l,p" << endl;
		for (i = 0; i < Top_NLM.nTop; i++)
		{
			file2 << Top_NLM.Tw[i] << ", ";
			file2 << Top_NLM.f[i] << ", ";
			file2 << Top_NLM.c[i] << ", ";
			file2 << Top_NLM.l[i] << ", ";
			file2 << Top_NLM.p[i] << endl;
		}

		// initialize simulated effective rainfall and streamflow arrays
		T_out = new T_Output_Arrays[Top_NLM.nTop];
		for (i = 0; i < Top_NLM.nTop; i++)
		{
			T_out[i].eR = new double[T_in.size];
			T_out[i].simStreamflow = new double[T_in.size];
			T_out[i].vq = new double[nsimMC_LM];
			T_out[i].vs = new double[nsimMC_LM];
			T_out[i].tq = new double[nsimMC_LM];
			T_out[i].ts = new double[nsimMC_LM];
		
			// calculate excess rainfall time series for each linear module calibration
			// if temperature data are available
			// calculate Wetness Time Constant
			if (T_in.bTMP)
			if (IHAC_version == 3)
				ihac.CalcWetnessTimeConst_Redesign(T_in.p_tmp,p_Tw, Top_NLM.Tw[i], Top_NLM.f[i], T_in.size);
			else
				ihac.CalcWetnessTimeConst(T_in.p_tmp, p_Tw, Top_NLM.Tw[i], Top_NLM.f[i], T_in.size);

			// calculate the Catchment Wetness Index
			if (IHAC_version == 3)
				ihac.CalcWetnessIndex_Redesign(p_Tw, T_in.p_pcp, p_WI, T_in.size);
			else
				ihac.CalcWetnessIndex(p_Tw,T_in.p_pcp, p_WI, Top_NLM.c[i], T_in.size);
			// calculate effective rainfall
			switch(IHAC_version)
			{
			case 0: // JAKEMAN & HORNBERGER (1993)
				Top_NLM.Sum_eRainMM[i] = ihac.CalcExcessRain(T_in.p_pcp, T_in.p_tmp, p_WI, T_out[i].eR, Top_NLM.sum_eRainGTpcp[i], T_in.size, bSnowModule);
				break;
			case 1: // Croke et al. (2005) Redesign
				Top_NLM.Sum_eRainMM[i] = ihac.CalcExcessRain_Redesign(T_in.p_pcp, T_in.p_tmp, p_WI, T_out[i].eR, Top_NLM.sum_eRainGTpcp[i], T_in.size,
																	Top_NLM.c[i], Top_NLM.l[i], Top_NLM.p[i], bSnowModule);
				break;
			} // end switch
		}	// end for i
ofstream fl("excessR.txt");
		for (i = 0; i < T_in.size; i++)
		{
			for (j = 0; j < Top_NLM.nTop; j++)
			{
				if (j < Top_NLM.nTop-1)
					fl << T_out[j].eR[i] << ", ";
				else
					fl << T_out[j].eR[i] << endl;
			}
		}
		delete[] T_NLP_valid;
		// ------------------------------------------------------------

		// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
		// Perform linear module calibration
		// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
		// number of simulations per nTop non-linear parameter set
		//int nsimMC_LM = 10000;
		ofstream file3("streamflow.txt");

		// initialize random function
		srand((unsigned) time(NULL)); // using time.h

		switch (storconf)
		{
			case 0: // single storage
				file3 << "j,i,NSE, a, b" << endl;
				// ----------------------------------
				// initializing storage parameters
				T_LM = new T_Linear_Parms[nsimMC_LM];
				for (i = 0; i < nsimMC_LM; i++)
				{
					// assign storage values
					T_LM[i].nStorages = 1;
					T_LM[i].a	= new double[1];
					T_LM[i].b	= new double[1];
					T_LM[i].a[0]= _Random_double(T_LB.a_lb, T_LB.a_ub);
					if (b_freebee)
						T_LM[i].b[0]= _Random_double(T_LB.b_lb, T_LB.b_ub);
					if (!b_freebee)
						T_LM[i].b[0]= 1 + T_LM[i].a[0];
				}
				// ----------------------------------
				for (j = 0; j < Top_NLM.nTop && Set_Progress(j,Top_NLM.nTop); j++)
				{

					for (i = 0; i < nsimMC_LM; i++)
					{
						//---------------------------------------------------------
						// simulate streamflow
						ihac.SimStreamflowSingle(T_out[j].eR,T_in.p_discharge[0],T_out[j].simStreamflow,1,T_LM[i].a[0],T_LM[i].b[0],T_in.size);
						//---------------------------------------------------------
						// calculate Nash-Sutcliffe efficiency
						NSE = ihac.CalcEfficiency(T_in.p_discharge,T_out[j].simStreamflow, T_in.size);
						//---------------------------------------------------------
						file3 << j << ", " << i << ", " << NSE << ", ";
						file3 << T_LM[i].a[0] << ", " << T_LM[i].b[0] << endl;

					}
				} // end for (j...)

				break;
			case 1: // two storages (parallel configuration)
				file3 << "j,i,NSE, aq, bq, as, bs, vq, vs" << endl;
				// ----------------------------------
				// initializing storage parameters
				T_LM = new T_Linear_Parms[nsimMC_LM];
				for (i = 0; i < nsimMC_LM; i++)
				{
					// assign storage values
					T_LM[i].nStorages = 2;
					T_LM[i].a	= new double[2];
					T_LM[i].b	= new double[2];
					// a[0] / b[0] = quick storage component
					T_LM[i].a[0]= _Random_double(T_LB.aq_lb, T_LB.aq_ub);
					// a[1] / b[1]= slow storage component
					T_LM[i].a[1]= _Random_double(T_LB.as_lb, T_LB.as_ub);
					T_LM[i].b[0]= _Random_double(T_LB.bq_lb, T_LB.bq_ub);
					if (b_freebee)
						T_LM[i].b[1]= _Random_double(T_LB.bs_lb, T_LB.bs_ub);
					if (!b_freebee)
						// here the value of b slow is calculated depending on the settings
						// of the other storage parameters
						T_LM[i].b[1]= (1 - (T_LM[i].b[0] / (1 + T_LM[i].a[0]))) * (1 + T_LM[i].a[1]);
				}
error << "after random LM" << endl;
				// ----------------------------------
				for (j = 0; j < Top_NLM.nTop && Set_Progress(j,Top_NLM.nTop); j++)
				{

					for (i = 0; i < nsimMC_LM; i++)
					{
						//---------------------------------------------------------
						// simulate streamflow
						ihac.SimStreamflow2Parallel(T_out[j].eR, T_out[j].simStreamflow, T_in.p_discharge[0],
													T_LM[i].a[0], T_LM[i].a[1], T_LM[i].b[0], T_LM[i].b[1], 
													T_out[j].vq[i], T_out[j].vs[i],
													IHAC_version, T_in.size, 1);
						//---------------------------------------------------------
						// calculate Nash-Sutcliffe efficiency
						NSE = ihac.CalcEfficiency(T_in.p_discharge,T_out[j].simStreamflow, T_in.size);
						//---------------------------------------------------------
						file3 << j << ", " << i << ", " << NSE << ", ";
						file3 << T_LM[i].a[0] << ", " << T_LM[i].b[0] << ", ";
						file3 << T_LM[i].a[1] << ", " << T_LM[i].b[1] << ", ";
						file3 << T_out[j].vq[i] << ", " << T_out[j].vs[i] << endl;

					}
				} // end for (j...)
				break;
		} // end switch (storconf)

		delete[] T_LM;
		delete[] T_out;


		// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


		// ------------------------------------------------------------

		
	}
	// **************************************************************


	delete[]	T_NLP;
	delete[]	p_bValidSim;

	delete[]	p_Tw;
	delete[]	p_WI;
	delete[]	p_ExcessRain;
	delete[]	p_SimDischarge;

	return(true);
} // end Function On_Execute
//---------------------------------------------------------
/*
void Cihacres_auto_cal::_Init_T_Cal(int steps)
{
	T_CP.nsteps1 = 50; //steps; //50; // Tw
	T_CP.nsteps2 = 30; //steps; //50; // f
	T_CP.nsteps3 = 30; //steps; //50; // c
	T_CP.ncombi = (T_CP.nsteps1 + 1) * (T_CP.nsteps2 + 1) * (T_CP.nsteps3 + 1);

	T_CP.lb1 = 0.0;
	T_CP.ub1 = 150;

	T_CP.lb2 = 0.01;
	T_CP.ub2 = 1.0;

	T_CP.lb3 = 0.001;
	T_CP.ub3 = 0.01;
}
*/

//---------------------------------------------------------------------

void Cihacres_auto_cal::_ZeroVariables()
{
	sum_obsDisMM		= 0.0;
	sum_simDisMM		= 0.0;
	sum_pcpMM			= 0.0;
	RunoffCoeff			= 0.0;
}
//---------------------------------------------------------

void Cihacres_auto_cal::_InitArrays(int IHAC_version, T_InputTable &T_in, T_SnowModule &T_SM)
{
	int size = T_in.size;

	T_in.date.resize(size);
	T_in.p_discharge = new double[size];
	T_in.p_dischargeMM = new double[size];
	T_in.p_pcp = new double[size];
	T_in.p_tmp = new double[size];
	p_Tw = new double[size];
	p_WI = new double[size];
	p_ExcessRain = new double[size];
	p_SimDischarge = new double[size];
	if (IHAC_version == 2) // snow module
	{
		T_SM.p_SnowStorage = new double[size];
		T_SM.p_MeltRate = new double[size];
	}
}
//---------------------------------------------------------------------

//T_NonLinear_Parms *Cihacres_auto_cal::_Set_PreCalibration1_Parms(T_Cal T_CP) //, T_NonLinear_Parms *T_nlp)
void Cihacres_auto_cal::_Set_PreCalibration1_Parms(T_Cal T_CP) 
{
	int		i, j, k;
	double	parm1, parm2, parm3;
	double	step1 = (T_CP.ub1 - T_CP.lb1) / T_CP.nsteps1;
	double	step2 = (T_CP.ub2 - T_CP.lb2) / T_CP.nsteps2;
	double	step3 = (T_CP.ub3 - T_CP.lb3) / T_CP.nsteps3;
	int		counter = 0;
	
	for (i = 0; i < T_CP.nsteps1+1; i++)
	{
		parm1 = T_CP.lb1 + i * step1;
		for (j = 0; j < T_CP.nsteps2+1; j++)
		{
			parm2 = T_CP.lb2 + j * step2;
			for (k = 0; k < T_CP.nsteps3+1; k++)
			{
				parm3 = T_CP.lb3 + k * step3;

				// Assign non-linear parameter settings
				T_NLP[counter].Tw	= parm1;
				T_NLP[counter].f	= parm2;
				T_NLP[counter].c	= parm3;
				T_NLP[counter].l	= 0.0;
				T_NLP[counter].p	= 1.0;
				counter++;
			} // end for(k...)

		} // end for(j...)

	} // end for(i...)

	//return(T_nlp);
}
//---------------------------------------------------------

//---------------------------------------------------------------------
//						Calculations of the Non-linear Module
//---------------------------------------------------------------------
void Cihacres_auto_cal::_CalcNonLinearModule(T_NonLinear_Parms* T_NLP, bool bTMP, int IHAC_version, int nsim, 
										double dev_eRainDis, double dev_eRainPCP, bool bSnowModule)
{
	int j;
	nValidSims = 0;
ofstream file("PreCal.txt");
file << "Tw, f, c, l, p, Dif_eRainDis, sum_eRainGTpcp" << endl;

	for (j = 0; j < nsim && Set_Progress(j,nsim); j++)
	{
		// if temperature data are available (bTMP = true), then adjust the rate
		// at which the catchment wetness declines in the absence of rainfall
		// to daily temperature
		if (bTMP)
		{
			if (IHAC_version == 3)
				ihac.CalcWetnessTimeConst_Redesign(T_in.p_tmp,p_Tw, T_NLP[j].Tw, T_NLP[j].f,T_in.size);
			else
				ihac.CalcWetnessTimeConst(T_in.p_tmp, p_Tw, T_NLP[j].Tw, T_NLP[j].f, T_in.size);
		}
		// calculate the catchment wetness index
		if (IHAC_version == 3)
			ihac.CalcWetnessIndex_Redesign(p_Tw, T_in.p_pcp, p_WI, T_in.size);
		else
			ihac.CalcWetnessIndex(p_Tw, T_in.p_pcp, p_WI, T_NLP[j].c, T_in.size);

		// calculate effective rainfall
		switch(IHAC_version)
		{
		case 0: // JAKEMAN & HORNBERGER (1993)
			T_NLP[j].Sum_eRainMM = ihac.CalcExcessRain(T_in.p_pcp, T_in.p_tmp, p_WI, p_ExcessRain, T_NLP[j].sum_eRainGTpcp, T_in.size,bSnowModule);
			break;
		case 3: // Croke et al. (2005) Redesign
			T_NLP[j].Sum_eRainMM = ihac.CalcExcessRain_Redesign(T_in.p_pcp, T_in.p_tmp, p_WI, p_ExcessRain, T_NLP[j].sum_eRainGTpcp, T_in.size,
													T_NLP[j].c, T_NLP[j].l, T_NLP[j].p, bSnowModule);
			break;

		} // end switch

		// calculate difference between discharge and excess rainfall
		T_NLP[j].Dif_eRainDis = sum_obsDisMM - T_NLP[j].Sum_eRainMM;
		if (T_NLP[j].Dif_eRainDis < 0.0) T_NLP[j].Dif_eRainDis *= -1.0;

		//---------------------------------------------------------
		// decide whether a simulation is valid or not
		//---------------------------------------------------------
		if ( T_NLP[j].sum_eRainGTpcp  / T_NLP[j].Sum_eRainMM * 100 <= dev_eRainPCP &&
			 T_NLP[j].Dif_eRainDis / sum_obsDisMM   * 100 <= dev_eRainDis   )
		{
			p_bValidSim[j] = true;
			nValidSims += 1;
			file << T_NLP[j].Tw << ", ";
			file << T_NLP[j].f << ", ";
			file << T_NLP[j].c << ", ";
			file << T_NLP[j].l << ", ";
			file << T_NLP[j].p << ", ";
			file << T_NLP[j].Dif_eRainDis << ", ";
			file << T_NLP[j].sum_eRainGTpcp << endl;
		}
/*
			file1 << T_NLP[j].Tw << ", ";
			file1 << T_NLP[j].f << ", ";
			file1 << T_NLP[j].c << ", ";
			file1 << T_NLP[j].Dif_eRainDis << ", ";
			file1 << T_NLP[j].sum_eRainGTpcp << endl;
*/
	} // end nsim calibration (non-linear module)
}

//---------------------------------------------------------------------
//			Calculations of the Non-linear Module (Monte Carlo)
//---------------------------------------------------------------------
void Cihacres_auto_cal::_CalcNonLinearModuleMC(T_NonLinear_Parms* T_nlp, T_nOccurence *T_occ, bool bTMP, int IHAC_version, int nsim, 
										double dev_eRainDis, double dev_eRainPCP)
{
	int j;
	nValidSims = 0;
	// initialize random function
	srand((unsigned) time(NULL)); // using time.h
	for (j = 0; j < nsim && Set_Progress(j,nsim); j++)
	{
		//T_nlp[j].Tw = _Random_double(T_occ[0].new_lb, T_occ[0].new_ub); // Tw
		T_nlp[j].Tw = _Random_double(1,15);
		//T_nlp[j].Tw = _Random_double(5, 25); // Tw
		//T_nlp[j].f	= _Random_double(T_occ[1].new_lb, T_occ[1].new_ub); // f
		T_nlp[j].f	= _Random_double(0,4); // f
		//T_nlp[j].f	= _Random_double(0.05, 0.3); // f
		T_nlp[j].c	= _Random_double(0.001, 0.01); // c
		//T_nlp[j].c	= _Random_double(T_occ[2].new_lb, T_occ[2].new_ub); // c

		T_nlp[j].l = 0;
		T_nlp[j].p = 1;

		// if temperature data are available (bTMP = true), then adjust the rate
		// at which the catchment wetness declines in the absence of rainfall
		// to daily temperature
		if (bTMP)
		{
			if (IHAC_version == 1)
				ihac.CalcWetnessTimeConst_Redesign(T_in.p_tmp,p_Tw, T_nlp[j].Tw, T_nlp[j].f,T_in.size);
			else
				ihac.CalcWetnessTimeConst(T_in.p_tmp, p_Tw, T_nlp[j].Tw, T_nlp[j].f, T_in.size);
		}
		// calculate the catchment wetness index
		if (IHAC_version == 1)
			ihac.CalcWetnessIndex_Redesign(p_Tw, T_in.p_pcp, p_WI, T_in.size);
		else
			ihac.CalcWetnessIndex(p_Tw, T_in.p_pcp, p_WI, T_nlp[j].c, T_in.size);

		// calculate effective rainfall
		switch(IHAC_version)
		{
		case 0: // JAKEMAN & HORNBERGER (1993)
			T_nlp[j].Sum_eRainMM = ihac.CalcExcessRain(T_in.p_pcp, T_in.p_tmp, p_WI, p_ExcessRain, T_nlp[j].sum_eRainGTpcp, T_in.size, bSnowModule);
			break;
		case 1: // Croke et al. (2005) Redesign
			T_nlp[j].Sum_eRainMM = ihac.CalcExcessRain_Redesign(T_in.p_pcp, T_in.p_tmp, p_WI, p_ExcessRain, T_nlp[j].sum_eRainGTpcp, T_in.size,
													T_nlp[j].c, T_nlp[j].l, T_nlp[j].p, bSnowModule);
			break;

		} // end switch

		// calculate difference between discharge and excess rainfall
		T_nlp[j].Dif_eRainDis = sum_obsDisMM - T_nlp[j].Sum_eRainMM;
		if (T_nlp[j].Dif_eRainDis < 0.0) T_nlp[j].Dif_eRainDis *= -1.0;

		//---------------------------------------------------------
		// decide whether a simulation is valid or not
		//---------------------------------------------------------
		if ( T_nlp[j].sum_eRainGTpcp  / T_nlp[j].Sum_eRainMM * 100 <= dev_eRainPCP &&
			 T_nlp[j].Dif_eRainDis / sum_obsDisMM   * 100 <= dev_eRainDis   )
		{
			p_bValidSim[j] = true;
			nValidSims += 1;
		}

	} // end nsim calibration (non-linear module)
}

//---------------------------------------------------------------------
// Count number of true values in bool array
//---------------------------------------------------------------------
int Cihacres_auto_cal::_CountNumberTrue(bool *array, int size)
{
	int counter = 0;
	// count number of true values
	for (int i = 0; i < size; i++)
	{
		if (array[i]) counter += 1;
	}
	return (counter);
}

//---------------------------------------------------------------------
// Select valid simulation parameter sets in a structure
//---------------------------------------------------------------------
//T_NonLinear_Parms *Cihacres_auto_cal::_Select_Valid_Sims(int nsim,
//														 T_NonLinear_Parms* T_NLP,
//														 T_NonLinear_Parms* valSims)
void Cihacres_auto_cal::_Select_Valid_Sims(int nsim, T_NonLinear_Parms* T_NLP)
{
	int j,k;
	k = 0;
	for (j = 0; j < nsim; j++)
	{
		if (p_bValidSim[j])
		{
			T_NLP_valid[k].c				= T_NLP[j].c;
			T_NLP_valid[k].Tw				= T_NLP[j].Tw;
			T_NLP_valid[k].f				= T_NLP[j].f;
			T_NLP_valid[k].l				= T_NLP[j].l;
			T_NLP_valid[k].p				= T_NLP[j].p;
			T_NLP_valid[k].Dif_eRainDis		= T_NLP[j].Dif_eRainDis;
			T_NLP_valid[k].sum_eRainGTpcp	= T_NLP[j].sum_eRainGTpcp;
			T_NLP_valid[k].Sum_eRainMM		= T_NLP[j].Sum_eRainMM;
			k++;
		}
	}
	//return(valSims);
}

//---------------------------------------------------------------------
// Initialize struct T_nOccurence
//---------------------------------------------------------------------
//T_nOccurence* Cihacres_auto_cal::_Init_Tocc(T_nOccurence *T_occ, T_Cal T_CP)
void Cihacres_auto_cal::_Init_Tocc(T_Cal T_CP)
{
	double	step;
	int		i,j;
	int		counter;

	// Parameter 1 = Tw
	T_occ[0].lb = T_CP.lb1;
	T_occ[0].ub = T_CP.ub1;
	T_occ[0].nsteps = T_CP.nsteps1;
	step = (T_occ[0].ub - T_occ[0].lb) / T_occ[0].nsteps;
	T_occ[0].val = new double[T_occ[0].nsteps+1];
	T_occ[0].n_occ = new int[T_occ[0].nsteps+1];
	T_occ[0].n_classes = _Calc_num_selection(T_occ[0]);
	T_occ[0].n_best = _round_i(T_occ[0].nsteps / T_occ[0].n_classes);
	// assign parameter setting to val[]	
	for (i = 0; i < T_occ[0].nsteps+1; i++)
	{
		T_occ[0].val[i] = T_occ[0].lb + i * step;
		counter = 0;
		// count number of occurences
		for (j = 0; j < nValidSims; j++)
		{
			if (T_NLP_valid[j].Tw == T_occ[0].val[i])
				counter += 1;
		}
		T_occ[0].n_occ[i] = counter;
	}

	// Parameter 2 = f
	T_occ[1].lb = T_CP.lb2;
	T_occ[1].ub = T_CP.ub2;
	T_occ[1].nsteps = T_CP.nsteps2;
	step = (T_occ[1].ub - T_occ[1].lb) / T_occ[1].nsteps;
	T_occ[1].val = new double[T_occ[1].nsteps+1];
	T_occ[1].n_occ = new int[T_occ[1].nsteps+1];
	T_occ[1].n_classes = _Calc_num_selection(T_occ[1]);
	T_occ[1].n_best = _round_i(T_occ[1].nsteps / T_occ[1].n_classes);
	for (i = 0; i < T_occ[1].nsteps+1; i++)
	{
		T_occ[1].val[i] = T_occ[1].lb + i * step;
		counter = 0;
		// count number of occurences
		for (j = 0; j < nValidSims; j++)
		{
			if (T_NLP_valid[j].f == T_occ[1].val[i])
				counter += 1;
		}
		T_occ[1].n_occ[i] = counter;
	}

	// Parameter 3 = c
	T_occ[2].lb = T_CP.lb3;
	T_occ[2].ub = T_CP.ub3;
	T_occ[2].nsteps = T_CP.nsteps3;
	step = (T_occ[2].ub - T_occ[2].lb) / T_occ[2].nsteps;
	T_occ[2].val = new double[T_occ[2].nsteps+1];
	T_occ[2].n_occ = new int[T_occ[2].nsteps+1];
	T_occ[2].n_classes = _Calc_num_selection(T_occ[2]);
	T_occ[2].n_best = _round_i(T_occ[2].nsteps / T_occ[2].n_classes);
	for (i = 0; i < T_occ[2].nsteps+1; i++)
	{
		T_occ[2].val[i] = T_occ[2].lb + i * step;
		counter = 0;
		// count number of occurences
		for (j = 0; j < nValidSims; j++)
		{
			if (T_NLP_valid[j].c == T_occ[2].val[i])
				counter += 1;
		}
		T_occ[2].n_occ[i] = counter;
	}

	//return (T_occ);
}

//---------------------------------------------------------------------
// Select best parameter settings in order to refine parameter ranges
//---------------------------------------------------------------------
void Cihacres_auto_cal::_Select_Best(T_nOccurence* T_occ)
{
	int i,j;
	for (i = 0; i < 3; i++)
	{
		for (j = 0; j < T_occ[i].nsteps+1; j++)
		{

		}
	}

}

//---------------------------------------------------------------------
// Calculate the number of best parameter settings
//---------------------------------------------------------------------
int Cihacres_auto_cal::_Calc_num_selection(T_nOccurence T_occ)
{
	// number of classes after STURGES
	// num Classes = 1 + 3.32 * log nValues
	return(_round_i(1 + 3.32 * log10(double(T_occ.nsteps + 1))));
}

//---------------------------------------------------------------------
// Find Parameter Indicees of best simulations
//---------------------------------------------------------------------
int* Cihacres_auto_cal::_FindTopIndex(int* array, int size, int top, int nSim)
{
	int largest			= 0;
	int n_occ			= 0;
	int temp_largest	= nSim + 1;
	int dif				= top;
	int setIndex		= 0;
	int start			= 0;
	int begin_i			= 0;
	int* TopIndex;
	TopIndex = new int[size];

	for (int i = 0; i < top; i++)
	{
		// find largest value in array below temp_largest
		largest = _largest_below(array, size, temp_largest);
		// counter number of occurence of largest value
		// this must be taken into account, because multiple
		// similar largest values are possible
		n_occ = _count_n_occ(array,size,largest);
		dif -= n_occ;

		if (dif >= 0) setIndex = n_occ;
		if (dif < 0)
			setIndex = n_occ - dif;

		begin_i = i;
		if (begin_i < top)
		{
			for (int j = begin_i; j < setIndex+begin_i; j++)
			{
				TopIndex[j] = _Get_Index(array, size, largest, start);
				start = TopIndex[j] + 1;
			}
		}
		i += setIndex-1;
		start = 0;
		temp_largest = largest;
		if (dif <= 0) i = top; // exit for (int i)-loop
	}
	return(TopIndex);
}
/*
int* Cihacres_auto_cal::_FindTopIndex(int* array, int size, int *TopIndex, int top, int nSim)
{
	int largest			= 0;
	int n_occ			= 0;
	int temp_largest	= nSim + 1;
	int dif				= top;
	int setIndex		= 0;
	int start			= 0;
	int begin_i			= 0;

	for (int i = 0; i < top; i++)
	{
		// find largest value in array below temp_largest
		largest = _largest_below(array, size, temp_largest);
		// counter number of occurence of largest value
		// this must be taken into account, because multiple
		// similar largest values are possible
		n_occ = _count_n_occ(array,size,largest);
		dif -= n_occ;

		if (dif >= 0) setIndex = n_occ;
		if (dif < 0)
			setIndex = n_occ - dif;

		begin_i = i;
		if (begin_i < top)
		{
			for (int j = begin_i; j < setIndex+begin_i; j++)
			{
				TopIndex[j] = _Get_Index(array, size, largest, start);
				start = TopIndex[j] + 1;
			}
		}
		i += setIndex-1;
		start = 0;
		temp_largest = largest;
		if (dif <= 0) i = top; // exit for (int i)-loop
	}
	return(TopIndex);
}
*/
//---------------------------------------------------------------------
// Assign top values
//---------------------------------------------------------------------
double* Cihacres_auto_cal::_AssignTopVal(T_nOccurence T_occ, int* topindex)
{
	double* out;
	out = new double[T_occ.n_best];

	for (int i = 0; i < T_occ.n_best; i++)
	{
		out[i] = T_occ.val[topindex[i]];
	}

	return(out);
}

//---------------------------------------------------------------------
// Assign top parameter settings from non-linear module calibration
//---------------------------------------------------------------------
T_Top_NLM_Sims Cihacres_auto_cal::_Assign_TopNLM(int nTop, T_NonLinear_Parms* NLM)
{
	int				i;
	int				ntop;
	int				id;
	double*			dif;
	double			except = -999999.9;
	T_Top_NLM_Sims	TopNLM;

	if (nValidSims >= nTop) ntop = nTop;
	else ntop = nValidSims;
	
	// assign Dif_eRainDis values of valid simulations to array dif
	dif = new double[nValidSims];
	for (i = 0; i < nValidSims; i++)
		dif[i] = NLM[i].Dif_eRainDis;

	// initialize struct T_Top_NLM_Sims
	TopNLM.nTop = ntop;
	TopNLM.Tw	= new double[ntop];
	TopNLM.f	= new double[ntop];
	TopNLM.c	= new double[ntop];
	TopNLM.l	= new double[ntop];
	TopNLM.p	= new double[ntop];
	TopNLM.Dif_eRainDis = new double[ntop];
	TopNLM.sum_eRainGTpcp = new double[ntop];
	TopNLM.Sum_eRainMM = new double[ntop];

	// Identify nTop simulations and
	// assign top parameter settings to TopNLM
	for (i = 0; i < ntop; i++)
	{
		id = _lowest_after(dif, nValidSims, except);
		TopNLM.Tw[i]	= NLM[id].Tw;
		TopNLM.f[i]		= NLM[id].f;
		TopNLM.c[i]		= NLM[id].c;
		TopNLM.l[i]		= NLM[id].l;
		TopNLM.p[i]		= NLM[id].p;
		TopNLM.Dif_eRainDis[i] = NLM[id].Dif_eRainDis;
		TopNLM.sum_eRainGTpcp[i] = NLM[id].sum_eRainGTpcp;
		TopNLM.Sum_eRainMM[i] = NLM[id].Sum_eRainMM;
		except			= dif[id];
	}

	return (TopNLM);
}

/////////////////////////////////////////////////////////////////////////////
//
//									HELPER
//
/////////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------
// round double value
//---------------------------------------------------------------------
int Cihacres_auto_cal::_round_i(double d)
{
  return d < 0 ? d -0.5 : d + 0.5;
}
//---------------------------------------------------------------------
double Cihacres_auto_cal::_largest(double* array, int size)
{
	double largest = array[0];;

	for (int i = 1; i < size; i++)
		if (largest < array[i]) largest = array[i];

	return largest;
}
//---------------------------------------------------------------------
double Cihacres_auto_cal::_lowest(double* array, int size)
{
	double lowest = array[0];
	for (int i = 1; i < size; i++)
		if (array[i] < lowest) lowest = array[i];

	return lowest;
}
//---------------------------------------------------------------------

int Cihacres_auto_cal::_lowest_after(double *array, int size, double exception)
{
	int		i;
	int		id = 0;
	double	tempMIN = array[0];

	for (i = 1; i < size; i++)
	{
		if (array[i] > exception && array[i] < tempMIN)
		{
			id = i;
			tempMIN = array[i];
		}
	}
	return (id);
}
//---------------------------------------------------------------------

int Cihacres_auto_cal::_largest_below(int* array, int size, int exception)
{
	int second_largest = 0;
	if (array[0] < exception) second_largest = array[0];

	//int second_largest = exception;

	for (int i = 1; i < size; i++)
		if (array[i] < exception && array[i] > second_largest)
			second_largest = array[i];

	return second_largest;
}
//---------------------------------------------------------------------
int Cihacres_auto_cal::_count_n_occ(int* array, int size, int val)
{
	int counter = 0;
	// count number of occurence of current largest value
	for (int i = 0; i < size; i++)
	{
		if (array[i] == val) counter++;
	}
	return(counter);
}
//---------------------------------------------------------------------
int Cihacres_auto_cal::_Get_Index(int* array, int size, int value, int start)
{
	for (int i = start; i < size; i++)
		if (array[i] == value)
			return(i);
	return (-1);
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
//                       Calculate random number
//---------------------------------------------------------------------
double Cihacres_auto_cal::_Random_double(double lb, double ub)
{
	// lb = lower boud
	// ub = upper bound
	double			random;

	random = rand()/(1.+RAND_MAX);
	random = lb + (ub - lb) * random;

	return random;
}


/////////////////////////////////////////////////////////////////////////////
//
//									DIALOGS
//
/////////////////////////////////////////////////////////////////////////////
void Cihacres_auto_cal::_CreateDialog1()
{
	CSG_Parameter	*pNode;
	CSG_String s;

	pNode = Parameters.Add_Table(
		NULL	, "TABLE"	, _TL("Table"),
		"",
		PARAMETER_INPUT
	);

	Parameters.Add_Table_Field(
		pNode	, "DATE_Field"	, _TL("Date Column"),
		"Select the column containing the Date"
	);

	Parameters.Add_Table_Field(
		pNode	, "DISCHARGE_Field"	, _TL("Streamflow Column"),
		"Select the Column containing Discharge Values"
	);

	Parameters.Add_Table_Field(
		pNode	, "PCP_Field"	, _TL("Precipitation Column"),
		"Select the Column containing precipitation Values"
	);

	Parameters.Add_Table_Field(
		pNode	, "TMP_Field"	, _TL("Temperature Column"),
		"Select the Column containing Temperature Values"
	);

	Parameters.Add_Value(
		pNode,	"USE_TMP",	_TL("Using temperature data?"),
		_TL("If checked, then temperature data are used."),
		PARAMETER_TYPE_Bool, true
	);

	pNode = Parameters.Add_Value(
		NULL,	"AREA",	_TL("Area of the Watershed in [km2]"),
		_TL("Area of the Watershed in [km2] used for unit conversion"),
		PARAMETER_TYPE_Double,
		100.0, 0.00001, true, 0.0, false
	);

	Parameters.Add_Value(
		pNode,	"DEV_ERAINDIS",	_TL("Maximum Deviation: ExcessRain-Streamflow [%]"),
		_TL("Maximum Deviation between estimated excess rainfall and observed streamflow volumes"),
		PARAMETER_TYPE_Double,
		5.0, 0.0, true, 50.0, true
	);

	Parameters.Add_Value(
		pNode,	"DEV_ERAINPCP",	_TL("Maximum Overestimation of ExcessRain [%]"),
		_TL("Maximum Sum of overestimated ExcessRain to total ExcessRain in [%]"),
		PARAMETER_TYPE_Double,
		5.0, 0.0, true, 50.0, true
	);
/*
	Parameters.Add_Value(
		pNode,	"STEP",	_TL("Parameter steps"),
		_TL(""),
		PARAMETER_TYPE_Int,
		50, 0, true, 500, true
	);
*/
	s.Printf("Node3", 3);
	pNode = Parameters.Add_Node(NULL,s,"Storage Configuration","");

	Parameters.Add_Choice(
		pNode	, "STORAGE"		, "Storage",
		"",
		CSG_String::Format("%s|%s|%s|",
			_TL("Single Storage"),			//  0
			_TL("Two Parallel Storages"),	//  1 
			_TL("Two Storages in Series")	//  2 
		)
	);

	s.Printf("Node4", 4);
	pNode = Parameters.Add_Node(NULL,s,"IHACRES Version","");

	Parameters.Add_Choice(
		pNode	, "IHACVERS"		, "IHACRES Version",
		"",
		CSG_String::Format("%s|%s|%s|%s|",
			_TL("Jakeman & Hornberger (1993)"),	//  0
			_TL("Kokkonen et al. (2003)"),		//  1 
			_TL("IHACRES v1.0 + SNOW MODULE"),	//  2 
			_TL("Croke et al. (2005)")			//	3
		)
	);

	Parameters.Add_Table_Output(
		NULL	, "TABLEout"	, _TL("Table"),
		""
	);
	Parameters.Add_Table_Output(
		NULL	, "TABLEparms"	, _TL("Table"),
		""
	);
}
///////////////////////////////////////////////////////////////////////
// DIALOG 2
///////////////////////////////////////////////////////////////////////
bool Cihacres_auto_cal::_CreateDialog2(int storconf, CSG_Table *pTable, T_InputTable &T_in, int IHAC_vers)
{
	CSG_String		s;
	CSG_Parameters	P;  // used to add Parameters in the second dialog
	CSG_Parameter	*pNode, *pNode1;


	//	Dialog design
	P.Set_Name(_TL("Choose Time Range"));

	s.Printf("Node1", 1);
	pNode = P.Add_Node(NULL,s,"Time Range","");

	s.Printf("FDAY" , 1-1);
	P.Add_String(pNode,s,_TL("First Day"),_TL(""),
				 pTable->Get_Record(0)->asString(T_in.dateField));

	s.Printf("LDAY" , 1-2);
	P.Add_String(pNode,s,_TL("Last Day"),_TL(""),
				 pTable->Get_Record(pTable->Get_Record_Count()-1)->asString(T_in.dateField));

	// Parameter Tw
	s.Printf("Node2", 2);
	pNode = P.Add_Node(NULL,s,"Non-Linear Module Parameters","");

	s.Printf("Node2-1", 2-1);
	pNode1 = P.Add_Node(pNode,s,"Parameter Tw","");

	P.Add_Value(
		pNode1,"TW_LB","Tw (lower bound)",
		_TL(""),
		PARAMETER_TYPE_Double,
		0.0, 0.0, true, 500.0, true
	);
	P.Add_Value(
		pNode1,"TW_UB", "Tw (upper bound)",
		_TL(""),
		PARAMETER_TYPE_Double,
		50.0, 0.0, true, 500.0, true
	);
	P.Add_Value(
		pNode1,"TW_Steps",_TL("Tw steps"),
		_TL(""),
		PARAMETER_TYPE_Int,
		10, 1, true, 500, true
	);

	if (T_in.bTMP)
	{
		// Parameter f
		s.Printf("Node2-2", 2-2);
		pNode1 = P.Add_Node(pNode,s,"Parameter f","");

		P.Add_Value(
			pNode1,"F_LB",_TL("f (lower bound)"),
			_TL(""),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 10.0, true
		);
		P.Add_Value(
			pNode1,"F_UB",_TL("f (upper bound)"),
			_TL(""),
			PARAMETER_TYPE_Double,
			5.0, 0.0, true, 20.0, true
		);
		P.Add_Value(
			pNode1,"F_Steps",_TL("f steps"),
			_TL(""),
			PARAMETER_TYPE_Int,
			10, 1, true, 500, true
		);
	}

	// Parameter c
	s.Printf("Node2-3", 2-3);
	pNode1 = P.Add_Node(pNode,s,"Parameter c","");

	P.Add_Value(
		pNode1,"C_LB",_TL("c (lower bound)"),
		_TL(""),
		PARAMETER_TYPE_Double,
		0.001, 0.0000001, true, 0.1, true
	);
	P.Add_Value(
		pNode1,"C_UB",_TL("c (upper bound)"),
		_TL(""),
		PARAMETER_TYPE_Double,
		0.01, 0.0000001, true, 0.1, true
	);
	P.Add_Value(
		pNode1,"C_Steps",_TL("c steps"),
		_TL(""),
		PARAMETER_TYPE_Int,
		10, 1, true, 500, true
	);

	if (IHAC_vers == 3) // Croke et al. 2005
	{
		s.Printf("Node2-32", 2-32);
		pNode1 = P.Add_Node(pNode,s,"Parameter l","");

		P.Add_Value(
			pNode1,"L_LB",_TL("l (lower bound)"),
			_TL(""),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
		P.Add_Value(
			pNode1,"L_UB",_TL("l (upper bound)"),
			_TL(""),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
		
		s.Printf("Node2-33", 2-33);
		pNode1 = P.Add_Node(pNode,s,"Parameter p","");

		P.Add_Value(
			pNode1,"P_LB",_TL("p (lower bound)"),
			_TL(""),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
		P.Add_Value(
			pNode1,"P_UB",_TL("p (upper bound)"),
			_TL(""),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
	}

	// number of Monte-Carlo-Simulation (non-linear module)
	s.Printf("Node2-4", 2-4);
	pNode1 = P.Add_Node(pNode,s,"Number of Monte-Carlo-Simulations","");

	P.Add_Value(
		pNode1,"NSIM_MC_NLM","Number of MC-Simulations",
		_TL("Number of Monte-Carlo-Simulations for the non-linear module"),
		PARAMETER_TYPE_Int,
		10000, 1, true, 10000000, true
	);

	// Linear Module Parameters
	s.Printf("Node3", 3);
	pNode = P.Add_Node(NULL,s,"Linear Module Parameters","");

		switch(storconf)
		{
		case 0: // single storage
			s.Printf("Node3-1", 3-1);
			pNode1 = P.Add_Node(pNode,s,"Parameter (a)","");

			P.Add_Value(
				pNode1,	"AFAC_lb",	_TL("(a) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.5, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AFAC_ub",	_TL("(a) [upper bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.9, -0.99, true, -0.01, true
				);

			s.Printf("Node3-2", 3-2);
			pNode1 = P.Add_Node(pNode,s,"Parameter (b)","");

			P.Add_Value(
				pNode1,	"BFAC_lb",	_TL("(b) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				0.0, 0.001, true, 1.0, true
			);
			P.Add_Value(
				pNode1,	"BFAC_ub",	_TL("(b) [upper bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				1.0, 0.001, true, 1.0, true
			);

			// number of Monte-Carlo-Simulation (non-linear module)
			s.Printf("Node3-4", 3-4);
			pNode1 = P.Add_Node(pNode,s,"Number of Monte-Carlo-Simulations","");

			P.Add_Value(
				pNode1,"NSIM_MC_LM","Number of MC-Simulations (non-linear module)",
				_TL("Number of Monte-Carlo-Simulations for the linear module"),
				PARAMETER_TYPE_Int,
				10000, 1, true, 10000000, true
			);
			break;

		case 1: // two parallel storages
			// Parameter a
			s.Printf("Node3-1", 3-1);
			pNode1 = P.Add_Node(pNode,s,"Parameter (a)","");

			P.Add_Value(
				pNode1,	"AQ_lb",	_TL("a(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.5, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AQ_ub",	_TL("a(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.9, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AS_lb",	_TL("a(s) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.5, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AS_ub",	_TL("a(s) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.9, -0.99, true, -0.01, true
			);

			// Parameter b
			s.Printf("Node3-2", 3-2);
			pNode1 = P.Add_Node(pNode,s,"Parameter (b)","");

			P.Add_Value(
				pNode1,	"BQ_lb",	_TL("b(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				0.0, 0.0, true, 1.0, true
			);
			P.Add_Value(
				pNode1,	"BQ_ub",	_TL("b(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				1.0, 0.0, true, 1.0, true
			);
			P.Add_Value(
				pNode1,	"BS_lb",	_TL("b(s) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				0.0, 0.0, true, 1.0, true
			);
			P.Add_Value(
				pNode1,	"BS_ub",	_TL("b(s) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				1.0, 0.0, true, 1.0, true
			);
			P.Add_Value(
				pNode1,"B_FREE","b-value independent",
				_TL("Usually the value of b is depending on the settings of a. If this box is checked b can be chosen independendly"),
				PARAMETER_TYPE_Bool, true
			);
			// number of Monte-Carlo-Simulation (non-linear module)
			s.Printf("Node3-4", 3-4);
			pNode1 = P.Add_Node(pNode,s,"Number of Monte-Carlo-Simulations","");

			P.Add_Value(
				pNode1,"NSIM_MC_LM","Number of MC-Simulations (linear module)",
				_TL("Number of Monte-Carlo-Simulations for the linear module"),
				PARAMETER_TYPE_Int,
				10000, 1, true, 10000000, true
			);
			break;
		case 2: // two storages in series
			break;
		} // end switch (storconf)


	// End of second dialog design


	if( SG_UI_Dlg_Parameters(&P, _TL("Choose Time Range")) )
	{
		// assign data from second dialog

		T_in.date1		= P(CSG_String::Format("FDAY",
							T_in.dateField).c_str())->asString();
		T_in.date2		= P(CSG_String::Format("LDAY",
							T_in.dischargeField).c_str())->asString();

		// assign Non-Linear Module Parameters
		T_CP.lb1		= P("TW_LB")		->asDouble();
		T_CP.ub1		= P("TW_UB")		->asDouble();
		T_CP.nsteps1	= P("TW_Steps")		->asInt();
		T_CP.lb2		= P("F_LB")			->asDouble();
		T_CP.ub2		= P("F_UB")			->asDouble();
		T_CP.nsteps2	= P("F_Steps")		->asInt();	
		T_CP.lb3		= P("C_LB")			->asDouble();
		T_CP.ub3		= P("C_UB")			->asDouble();
		T_CP.nsteps3	= P("C_Steps")		->asInt();
		T_CP.ncombi		= (T_CP.nsteps1 + 1) * (T_CP.nsteps2 + 1) * (T_CP.nsteps3 + 1);
		nsimMC_NLM		= P("NSIM_MC_NLM")	->asInt();
		if (IHAC_vers == 3)
		{
			// Hier Parameter l und p einlesen
		}

		// assign Linear Module Parameters
		nsimMC_LM		= P("NSIM_MC_LM")	->asInt();
		switch (storconf)
		{
			case 0: // single storage
				T_LB.a_lb		= P("AFAC_lb")		->asDouble();
				T_LB.a_ub		= P("AFAC_ub")		->asDouble();
				T_LB.b_lb		= P("BFAC_lb")		->asDouble();
				T_LB.b_ub		= P("BFAC_ub")		->asDouble();
				break;

			case 1: // two storages (parallel configuration)
				T_LB.aq_lb	= P("AQ_lb")		->asDouble();
				T_LB.aq_ub	= P("AQ_ub")		->asDouble();
				T_LB.bq_lb	= P("BQ_lb")		->asDouble();
				T_LB.bq_ub	= P("BQ_ub")		->asDouble();
				T_LB.as_lb	= P("AS_lb")		->asDouble();
				T_LB.as_ub	= P("AS_ub")		->asDouble();
				T_LB.bs_lb	= P("BS_lb")		->asDouble();
				T_LB.bs_ub	= P("BS_ub")		->asDouble();
				break;

			case 2: // two storages (seriel configuration)
				break;
		} // end switch (storconf)




		return(true);
	} else {
		return(false);
	}
}
