/*!\file Matpar.c
 * \brief: implementation of the Matpar object
 */

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../classes.h"
#include "../../shared/shared.h"

/*Matpar constructors and destructor*/
/*FUNCTION Matpar::Matpar() {{{*/
Matpar::Matpar(){
	return;
}
/*}}}*/
/*FUNCTION Matpar::Matpar(int matpar_mid,IoModel* iomodel){{{*/
Matpar::Matpar(int matpar_mid, IoModel* iomodel){

	bool ispdd;
	bool isefficientlayer;
	int  hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
	iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);

	this->mid = matpar_mid;
	iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
	iomodel->Constant(&this->rho_water,MaterialsRhoWaterEnum);
	iomodel->Constant(&this->rho_freshwater,MaterialsRhoFreshwaterEnum);
	iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum);
	iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
	iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum);
	iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum);
	iomodel->Constant(&this->beta,MaterialsBetaEnum);
	iomodel->Constant(&this->meltingpoint,MaterialsMeltingpointEnum);
	iomodel->Constant(&this->referencetemperature,ConstantsReferencetemperatureEnum);
	iomodel->Constant(&this->mixed_layer_capacity,MaterialsMixedLayerCapacityEnum);
	iomodel->Constant(&this->thermal_exchange_velocity,MaterialsThermalExchangeVelocityEnum);
	iomodel->Constant(&this->g,ConstantsGEnum);

	if(ispdd){
		iomodel->Constant(&this->desfac,SurfaceforcingsDesfacEnum);
		iomodel->Constant(&this->s0p,SurfaceforcingsS0pEnum);
	}

	if(hydrology_model==HydrologyshreveEnum){
		iomodel->Constant(&this->hydro_CR,HydrologyshreveCREnum);
		iomodel->Constant(&this->hydro_kn,HydrologyshreveKnEnum);
		iomodel->Constant(&this->hydro_n,HydrologyshreveNEnum);
		iomodel->Constant(&this->hydro_p,HydrologyshrevePEnum);
		iomodel->Constant(&this->hydro_q,HydrologyshreveQEnum);
	}
	else if(hydrology_model==HydrologydcEnum){
		iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
		iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);
		iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum);
		iomodel->Constant(&this->sediment_transmitivity,HydrologydcSedimentTransmitivityEnum);
		iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
		iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);

		if(isefficientlayer){
				iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
				iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
				iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum);
				iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum);
		}
	}
	else{
		_error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
	}

	/*gia: */
	iomodel->Constant(&this->lithosphere_shear_modulus,MaterialsLithosphereShearModulusEnum);
	iomodel->Constant(&this->lithosphere_density,MaterialsLithosphereDensityEnum);
	iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
	iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);

	this->inputs=NULL; /*not used here*/
}
/*}}}*/
/*FUNCTION Matpar::~Matpar() {{{*/
Matpar::~Matpar(){
	return;
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION Matpar::Echo {{{*/
void Matpar::Echo(void){

	_printf_("Matpar:\n");
	_printf_("   mid: " << mid << "\n");
	_printf_("   rho_ice: " << rho_ice << "\n");
	_printf_("   rho_water: " << rho_water << "\n");
	_printf_("   rho_freshwater: " << rho_freshwater << "\n");
	_printf_("   mu_water: " << mu_water << "\n");
	_printf_("   heatcapacity: " << heatcapacity << "\n");
	_printf_("   thermalconductivity: " << thermalconductivity << "\n");
	_printf_("   latentheat: " << latentheat << "\n");
	_printf_("   beta: " << beta << "\n");
	_printf_("   meltingpoint: " << meltingpoint << "\n");
	_printf_("   referencetemperature: " << referencetemperature << "\n");
	_printf_("   mixed_layer_capacity: " << mixed_layer_capacity << "\n");
	_printf_("   thermal_exchange_velocity: " << thermal_exchange_velocity << "\n");
	_printf_("   g: " << g << "\n");
	_printf_("   desfac: " << desfac << "\n");
	_printf_("   s0p: " << s0p << "\n");
	return;
}
/*}}}*/
/*FUNCTION Matpar::DeepEcho {{{*/
void Matpar::DeepEcho(void){

	this->Echo();
}		
/*}}}*/
/*FUNCTION Matpar::Id {{{*/
int    Matpar::Id(void){ return mid; }
/*}}}*/
/*FUNCTION Matpar::ObjectEnum{{{*/
int Matpar::ObjectEnum(void){

	return MatparEnum;

}
/*}}}*/
/*FUNCTION Matpar::copy {{{*/
Object* Matpar::copy() {
	return new Matpar(*this); 
}
/*}}}*/

/*Update virtual functions definitions:*/
/*FUNCTION Matpar::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
void   Matpar::InputUpdateFromVector(IssmDouble* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matpar::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
void   Matpar::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matpar::InputUpdateFromMatrixDakota(IssmDouble* vector, int name, int type) {{{*/
void  Matpar::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matpar::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/
void   Matpar::InputUpdateFromConstant(IssmDouble constant, int name){

	switch(name){
		case MaterialsRhoIceEnum:
			this->rho_ice=constant;
			break;
		case MaterialsRhoWaterEnum:
			this->rho_water=constant;
			break;
		case MaterialsRhoFreshwaterEnum:
			this->rho_freshwater=constant;
			break;
		case MaterialsMuWaterEnum:
			this->mu_water=constant;
			break;
		case MaterialsHeatcapacityEnum:
			this->heatcapacity=constant;
			break;
		case MaterialsThermalconductivityEnum:
			this->thermalconductivity=constant;
			break;
		case  MaterialsLatentheatEnum:
			this->latentheat=constant;
			break;
		case  MaterialsBetaEnum:
			this->beta=constant;
			break;
		case  MaterialsMeltingpointEnum:
			this->meltingpoint=constant;
			break;
		case  ConstantsReferencetemperatureEnum:
			this->referencetemperature=constant;
			break;
		case  MaterialsMixedLayerCapacityEnum:
			this->mixed_layer_capacity=constant;
			break;
		case  MaterialsThermalExchangeVelocityEnum:
			this->thermalconductivity=constant;
			break;
		case  ConstantsGEnum:
			this->g=constant;
			break;
  	        case  SurfaceforcingsDesfacEnum:
			this->desfac=constant;
			break;
		case SurfaceforcingsS0pEnum:
			this->s0p=constant;
			break;
		default: 
			break;
	}

}
/*}}}*/
/*FUNCTION Matpar::InputUpdateFromConstant(int constant, int name) {{{*/
void   Matpar::InputUpdateFromConstant(int constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matpar::InputUpdateFromConstant(bool constant, int name) {{{*/
void   Matpar::InputUpdateFromConstant(bool constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matpar::InputUpdateFromSolution{{{*/
void   Matpar::InputUpdateFromSolution(IssmDouble* solution){
	/*Nothing updated yet*/
}
/*}}}*/

/*Matpar management: */
/*FUNCTION Matpar::InputDuplicate{{{*/
void  Matpar::InputDuplicate(int original_enum,int new_enum){

	   /*Call inputs method*/
	   if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);

}
/*}}}*/
/*FUNCTION Matpar::Configure {{{*/
void  Matpar::Configure(Elements* elementsin){

	/*nothing done yet!*/

}
/*}}}*/
/*FUNCTION Matpar::GetBeta {{{*/
IssmDouble Matpar::GetBeta(){
	return beta;
}
/*}}}*/
/*FUNCTION Matpar::GetG {{{*/
IssmDouble Matpar::GetG(){
	return g;
}
/*}}}*/
/*FUNCTION Matpar::GetHeatCapacity {{{*/
IssmDouble Matpar::GetHeatCapacity(){
	return heatcapacity;
}
/*}}}*/
/*FUNCTION Matpar::GetLatentHeat {{{*/
IssmDouble Matpar::GetLatentHeat(){
	return latentheat;
}
/*}}}*/
/*FUNCTION Matpar::GetMeltingPoint {{{*/
IssmDouble Matpar::GetMeltingPoint(){
	return meltingpoint;
}
/*}}}*/
/*FUNCTION Matpar::GetReferenceTemperature {{{*/
IssmDouble Matpar::GetReferenceTemperature(){
	return referencetemperature;
}
/*}}}*/
/*FUNCTION Matpar::GetMixedLayerCapacity {{{*/
IssmDouble Matpar::GetMixedLayerCapacity(){
	return mixed_layer_capacity;
}
/*}}}*/
/*FUNCTION Matpar::GetRhoIce {{{*/
IssmDouble Matpar::GetRhoIce(){

	return rho_ice;
}
/*}}}*/
/*FUNCTION Matpar::GetRhoWater {{{*/
IssmDouble Matpar::GetRhoWater(){
	return rho_water;
}
/*}}}*/
/*FUNCTION Matpar::GetRhoFreshwater {{{*/
IssmDouble Matpar::GetRhoFreshwater(){
	return rho_freshwater;
}
/*}}}*/
/*FUNCTION Matpar::GetMuWater {{{*/
IssmDouble Matpar::GetMuWater(){
	return mu_water;
}
/*}}}*/
/*FUNCTION Matpar::GetDesFac {{{*/
IssmDouble Matpar::GetDesFac(){
	return desfac;
}
/*}}}*/
/*FUNCTION Matpar::GetS0p {{{*/
IssmDouble Matpar::GetS0p(){
	return s0p;
}
/*}}}*/
/*FUNCTION Matpar::GetThermalConductivity {{{*/
IssmDouble Matpar::GetThermalConductivity(){
	return thermalconductivity;
}
/*}}}*/
/*FUNCTION Matpar::GetThermalExchangeVelocity {{{*/
IssmDouble Matpar::GetThermalExchangeVelocity(){
	return thermal_exchange_velocity;
}
/*}}}*/
/*FUNCTION Matpar::GetHydrologyKn {{{*/
IssmDouble Matpar::GetHydrologyKn(){
	return hydro_kn;		 
}		 
/*}}}*/			 
/*FUNCTION Matpar::GetHydrologyP {{{*/			 
IssmDouble Matpar::GetHydrologyP(){		 
	return hydro_p;			 
}		 
/*}}}*/			 
/*FUNCTION Matqar::GetHydrologyQ {{{*/			 
IssmDouble Matpar::GetHydrologyQ(){		 
	return hydro_q;			 
}		 
/*}}}*/			 
/*FUNCTION Matpar::GetHydrologyCR {{{*/		 
IssmDouble Matpar::GetHydrologyCR(){		 
	return hydro_CR;		 
}		 
/*}}}*/			 
/*FUNCTION Matpar::GetHydrologyN {{{*/			 
IssmDouble Matpar::GetHydrologyN(){		 
	return hydro_n;			 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetSedimentStoring {{{*/
IssmDouble Matpar::GetSedimentStoring(){
	return this->rho_freshwater* this->g* this->sediment_porosity* this->sediment_thickness*
    ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity));		 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetEplStoring {{{*/
IssmDouble Matpar::GetEplStoring(){
	return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness*
    ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity));		 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetSedimentTransitivity {{{*/
IssmDouble Matpar::GetSedimentTransmitivity(){
	return sediment_transmitivity;		 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetSedimentThickness {{{*/
IssmDouble Matpar::GetSedimentThickness(){
	return sediment_thickness;		 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetEplTransitivity {{{*/
IssmDouble Matpar::GetEplTransmitivity(){
	return epl_transmitivity;		 
}		 
/*}}}*/			 
/*FUNCTION Matpar::TMeltingPoint {{{*/
IssmDouble Matpar::TMeltingPoint(IssmDouble pressure){
	return meltingpoint-beta*pressure;
}
/*}}}*/
/*FUNCTION Matpar::PureIceEnthalpy{{{*/
IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){
	return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
}
/*}}}*/
/*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{*/
IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
	if(enthalpy<PureIceEnthalpy(pressure)){
		return thermalconductivity/(rho_ice*heatcapacity);
	}
	else{
		return 1.045*1.e-4/rho_ice; // K0=1.045*1e-4 from Aschwanden 2012. TODO: fetch K0 from model
	}
}
/*}}}*/
/*FUNCTION Matpar::EnthalpyToThermal {{{*/
void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){

	/*Ouput*/
	IssmDouble temperature,waterfraction;

	if(enthalpy<PureIceEnthalpy(pressure)){
		temperature=referencetemperature+enthalpy/heatcapacity;
		waterfraction=0;
	}
	else{
		temperature=TMeltingPoint(pressure);
		waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
	}

	/*Assign output pointers:*/
	*pwaterfraction=waterfraction;
	*ptemperature=temperature;
}
/*}}}*/
/*FUNCTION Matpar::ThermalToEnthalpy {{{*/
void Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){

	/*Ouput*/
	IssmDouble enthalpy;

	if(temperature<TMeltingPoint(pressure)){
		enthalpy=heatcapacity*(temperature-referencetemperature);
	}
	else{
		enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
	}

	/*Assign output pointers:*/
	*penthalpy=enthalpy;
}
/*}}}*/

/*FUNCTION Matpar::GetLithosphereShearModulus {{{*/			 
IssmDouble Matpar::GetLithosphereShearModulus(){		 
	return lithosphere_shear_modulus;			 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetLithosphereDensity {{{*/			 
IssmDouble Matpar::GetLithosphereDensity(){		 
	return lithosphere_density;			 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetMantleDensity {{{*/			 
IssmDouble Matpar::GetMantleDensity(){		 
	return mantle_density;			 
}		 
/*}}}*/ 
/*FUNCTION Matpar::GetMantleShearModulus {{{*/			 
IssmDouble Matpar::GetMantleShearModulus(){		 
	return mantle_shear_modulus;			 
}		 
/*}}}*/ 
/*FUNCTION Matpar::IsInput{{{*/
bool Matpar::IsInput(int name){
	return false;
}
/*}}}*/
