/*
 *  SNOWPACK stand-alone
 *
 *  Copyright WSL Institute for Snow and Avalanche Research SLF, DAVOS, SWITZERLAND
*/
/*  This file is part of Snowpack.
    Snowpack is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Snowpack is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Snowpack.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <snowpack/Aggregate.h>

/************************************************************
 * static section                                           *
 ************************************************************/

const double Aggregate::limit_dry     = 0.001; ///< Distinguishes between dry and wet snow layers (1)
const double Aggregate::diff_theta_w  = 0.7;   ///< Maximum water difference for aggregation (% by volume)
const double Aggregate::diff_jul      = 1.0;   ///< Maximum  age difference for aggregation (d)
const double Aggregate::diff_dg       = 0.25;  ///< Maximum  grain size difference for aggregation (mm)
const double Aggregate::diff_dg_rel   = 0.2;   ///< Maximum  relative grain size difference for aggregation (mm)
const double Aggregate::diff_sp       = 0.2;   ///< Maximum  sphericity difference for aggregation (1)
const double Aggregate::diff_dd       = 0.2;   ///< Maximum  dendricity difference for aggregation (1)
const double Aggregate::min_l_element = 0.3;   ///< Minimum length of element to be kept separate (cm)

/**
 * @brief Eliminate the "empty" layers shift the remaining layers to form a compact snowpack
 * @param nL_ini Iinitial number of layers prior to aggregation
 * @param Pdata A vector of SnowProfileLayer
 */
void Aggregate::shift(const size_t& nL_ini, std::vector<SnowProfileLayer>& Pdata)
{
	for (size_t ll=1, l_new=1; ll<nL_ini; ll++) {
		if (Pdata[ll].height != Constants::undefined) {
			if (l_new != ll) {
				Pdata[l_new].layerDate = Pdata[ll].layerDate;
				Pdata[l_new].height = Pdata[ll].height;
				Pdata[l_new].rho = Pdata[ll].rho;
				Pdata[l_new].T = Pdata[ll].T;
				Pdata[l_new].gradT = Pdata[ll].gradT;
				Pdata[l_new].strain_rate = Pdata[ll].strain_rate;
				Pdata[l_new].theta_i = Pdata[ll].theta_i;
				Pdata[l_new].theta_w = Pdata[ll].theta_w;
				Pdata[l_new].grain_size = Pdata[ll].grain_size;
				Pdata[l_new].dendricity = Pdata[ll].dendricity;
				Pdata[l_new].sphericity = Pdata[ll].sphericity;
				Pdata[l_new].ogs = Pdata[ll].ogs;
				Pdata[l_new].bond_size = Pdata[ll].bond_size;
				Pdata[l_new].coordin_num = Pdata[ll].coordin_num;
				Pdata[l_new].marker = Pdata[ll].marker;
				Pdata[l_new].hard = Pdata[ll].hard;
			}
			l_new++;
		}
	}
}

/**
* @brief Decide whether to join two similar layers
* @param l_upper index of upper layer
* @param Pdata A vector of SnowProfileLayer
*/
bool Aggregate::joinSimilarLayers(const size_t& l_upper, std::vector<SnowProfileLayer>& Pdata)
{
	const unsigned int l_lower = l_upper-1;

	if ((Pdata[l_upper].theta_w < limit_dry) || (Pdata[l_lower].theta_w < limit_dry)) {
		if (fabs(Pdata[l_upper].theta_w - Pdata[l_lower].theta_w) > limit_dry)
			return false;
	} else {
		if (fabs(Pdata[l_upper].theta_w - Pdata[l_lower].theta_w) > diff_theta_w)
			return false;
	}

	if (fabs(Pdata[l_upper].layerDate.getJulianDate() - Pdata[l_lower].layerDate.getJulianDate()) > diff_jul)
		return false;

	if (Pdata[l_upper].marker != Pdata[l_lower].marker)
		return false;

	// do not combine layers which are of quite different hardness 020917;Fz
	if (fabs(Pdata[l_upper].hard - Pdata[l_lower].hard) > 1.0)
		return false;

	if ((Pdata[l_upper].dendricity == 0) && (Pdata[l_lower].dendricity == 0)){
		if ( fabs(Pdata[l_upper].sphericity - Pdata[l_lower].sphericity) > diff_sp)
			return false;

		if (fabs(Pdata[l_upper].grain_size - Pdata[l_lower].grain_size) > MAX(diff_dg, diff_dg_rel * Pdata[l_upper].grain_size))
			return false;
	} else {
		if (fabs(Pdata[l_upper].sphericity - Pdata[l_lower].sphericity) > diff_sp)
			return false;

		if (fabs(Pdata[l_upper].dendricity - Pdata[l_lower].dendricity) > diff_dd)
			return false;
	}

	return true;
}

/**
 * @brief Decide whether to merge thin layers
 * @param l_lower Index of upper layer
 * @param Pdata A vector of SnowProfileLayer
 */
bool Aggregate::mergeThinLayer(const size_t& l_lower, std::vector<SnowProfileLayer>& Pdata)
{
	const size_t l_upper = l_lower+1;

	// if a dry layer is involved
	if ((Pdata[l_lower].theta_w < limit_dry) || (Pdata[l_upper].theta_w < limit_dry)){
		// do not combine dry with moist layers
		if (fabs(Pdata[l_lower].theta_w - Pdata[l_upper].theta_w) > limit_dry)
			return false;

		// do not combine layers which are of quite different age
		if (fabs(Pdata[l_lower].layerDate.getJulianDate() - Pdata[l_upper].layerDate.getJulianDate()) > 2 * diff_jul)
			return false;

		// do not combine layers with different grain classes
		if (Pdata[l_lower].marker != Pdata[l_upper].marker) {
			return false;
		}
		// do not combine layers which are of quite different hardness
		if (fabs(Pdata[l_lower].hard - Pdata[l_upper].hard) > 1.5) {
			return false;
		}
	}
	// for two wet layers
	else {
		if ((Pdata[l_lower].grain_size < 0.75) || (Pdata[l_lower].sphericity < 0.5) ||
			(Pdata[l_upper].grain_size < 0.75) || (Pdata[l_upper].sphericity < 0.5)) {
			if (Pdata[l_lower].marker != Pdata[l_upper].marker) {
				return false;
			}
		}
		else {
			if (!((((Pdata[l_lower].marker > 9) &&  (Pdata[l_lower].marker < 13)) ||
				((Pdata[l_lower].marker > 19) &&  (Pdata[l_lower].marker < 23)) ) &&
				( ((Pdata[l_lower].marker > 9) &&  (Pdata[l_lower].marker < 13)) ||
				((Pdata[l_lower].marker > 19) &&  (Pdata[l_lower].marker < 23))))) {

				if (Pdata[l_lower].marker != Pdata[l_upper].marker)
					return false;
			}
		}
	}
	return true;
}

/**
 * @brief Aggregate snow profile layers and compute the grain class
 * @note Leave top element untouched
 * @param Pdata A vector of SnowProfileLayer
 */
size_t Aggregate::aggregate(std::vector<SnowProfileLayer>& Pdata)
{
	size_t nL_ini = Pdata.size();
	size_t nL = nL_ini;

	// Initialize number of layers and aggregate only if more than 5 layers
	if (nL > 5) {
		// First Run - aggregate similar layers
		// keep track of the coordinates and length of elements
		double L0_lower = (Pdata[nL_ini-2].height -  Pdata[nL_ini-3].height);
		for (size_t l_upper=nL_ini-2; l_upper > 0; l_upper--) {
			const double L0_upper = L0_lower;
			const size_t l_lower = l_upper-1;
			if (l_lower > 0) {
				L0_lower = (Pdata[l_lower].height -  Pdata[l_lower-1].height);
			} else {
				L0_lower = Pdata[l_lower].height;
			}
			// if two layers are similar combine them; keep SH though
			if ((Pdata[l_lower].marker != 3) && (Pdata[l_upper].marker != 3)) {
				if (joinSimilarLayers(l_upper, Pdata)) {
					nL--;
					Pdata[l_lower].average(L0_lower, L0_upper, Pdata[l_upper]);
					Pdata[l_upper].height = Constants::undefined;
					L0_lower += L0_upper;
				}
			}
		}

		shift(nL_ini, Pdata);
		nL_ini = nL;

		// Second Run - aggregate remaining very thin layers
		if (nL_ini > 2) {
			bool flag = false;
			double L0_lower = (Pdata[nL_ini-2].height -  Pdata[nL_ini-3].height);
			for(size_t l_upper = nL_ini-2; l_upper > 0; l_upper--) {
				const size_t l_lower = l_upper-1;
				const double L0_upper = L0_lower;
				if (l_lower > 0) {
					L0_lower = (Pdata[l_lower].height -  Pdata[l_lower-1].height);
				} else {
					L0_lower = Pdata[l_lower].height;
				}
				if ((Pdata[l_lower].marker != 3) && (Pdata[l_upper].marker != 3)) {
					// trick to try to join with upper or lower level -> use flag to mark thin layer
					if (flag || (L0_lower < (sqrt(Pdata[nL_ini-1].height-Pdata[l_lower].height)/4.))
						         || (L0_lower < min_l_element)) {
						// if two layers are similar or one layer is very very small combine them
						if (mergeThinLayer(l_upper, Pdata)
							    || (L0_lower < min_l_element) || (L0_upper < min_l_element)) {
							nL--;
							Pdata[l_lower].average(L0_lower, L0_upper, Pdata[l_upper]);
							Pdata[l_upper].height = Constants::undefined;
							L0_lower += L0_upper;
							flag = false;
						} else {
							flag = true;
						}
					} else {
						flag = false;
					}
				} // if not surface hoar layer
				else {
					flag = false;
				}
			}  // for all elements
			shift(nL_ini, Pdata);
		} // if nL_ini > 2
	} // if more than 5 layers

	// Update snow type
	for(size_t ll=0; ll<nL; ll++) {
		Pdata[ll].type = ElementData::snowType(Pdata[ll].dendricity, Pdata[ll].sphericity,
                                               Pdata[ll].grain_size, Pdata[ll].marker, Pdata[ll].theta_w/100.,
                                               ElementData::snowResidualWaterContent(Pdata[ll].theta_i/100.));
	}
	return (nL);
}
