/*!\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 <stdio.h>
#include <string.h>
#include "../objects.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
		
/*Matpar constructors and destructor*/
/*FUNCTION Matpar::Matpar() {{{1*/
Matpar::Matpar(){
	return;
}
/*}}}1*/
/*FUNCTION Matpar::Matpar(int matpar_mid,IoModel* iomodel){{{1*/
Matpar::Matpar(int matpar_mid, IoModel* iomodel){

	this->mid                       = matpar_mid;
	iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
	iomodel->Constant(&this->rho_water,MaterialsRhoWaterEnum);
	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);
	
	iomodel->Constant(&this->hydro_CR,HydrologyCREnum);
	iomodel->Constant(&this->kn,HydrologyKnEnum);
	iomodel->Constant(&this->hydro_n,HydrologyNEnum);
	iomodel->Constant(&this->hydro_p,HydrologyPEnum);
	iomodel->Constant(&this->hydro_q,HydrologyQEnum);
}
/*}}}1*/
/*FUNCTION Matpar::~Matpar() {{{1*/
Matpar::~Matpar(){
	return;
}
/*}}}1*/

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

	printf("Matpar:\n");
	printf("   mid: %i\n",mid);
	printf("   rho_ice: %g\n",rho_ice);
	printf("   rho_water: %g\n",rho_water);
	printf("   mu_water: %g\n",mu_water);
	printf("   heatcapacity: %g\n",heatcapacity);
	printf("   thermalconductivity: %g\n",thermalconductivity);
	printf("   latentheat: %g\n",latentheat);
	printf("   beta: %g\n",beta);
	printf("   meltingpoint: %g\n",meltingpoint);
	printf("   referencetemperature: %g\n",referencetemperature);
	printf("   mixed_layer_capacity: %g\n",mixed_layer_capacity);
	printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
	printf("   g: %g\n",g);
	return;
}
/*}}}1*/
/*FUNCTION Matpar::DeepEcho {{{1*/
void Matpar::DeepEcho(void){

	printf("Matpar:\n");
	printf("   mid: %i\n",mid);
	printf("   rho_ice: %g\n",rho_ice);
	printf("   rho_water: %g\n",rho_water);
	printf("   mu_water: %g\n",mu_water);
	printf("   heatcapacity: %g\n",heatcapacity);
	printf("   thermalconductivity: %g\n",thermalconductivity);
	printf("   latentheat: %g\n",latentheat);
	printf("   beta: %g\n",beta);
	printf("   meltingpoint: %g\n",meltingpoint);
	printf("   referencetemperature: %g\n",referencetemperature);
	printf("   mixed_layer_capacity: %g\n",mixed_layer_capacity);
	printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
	printf("   g: %g\n",g);
	return;
}		
/*}}}1*/
/*FUNCTION Matpar::Id {{{1*/
int    Matpar::Id(void){ return mid; }
/*}}}1*/
/*FUNCTION Matpar::MyRank {{{1*/
int    Matpar::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}1*/
/*FUNCTION Matpar::ObjectEnum{{{1*/
int Matpar::ObjectEnum(void){

	return MatparEnum;

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

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

	switch(name){
		case MaterialsRhoIceEnum:
			this->rho_ice=constant;
			break;

		case MaterialsRhoWaterEnum:
			this->rho_water=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;

		default: 
			break;
	}

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

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

	/*nothing done yet!*/

}
/*}}}*/
/*FUNCTION Matpar::GetBeta {{{1*/
double Matpar::GetBeta(){
	return beta;
}
/*}}}1*/
/*FUNCTION Matpar::GetG {{{1*/
double Matpar::GetG(){
	return g;
}
/*}}}1*/
/*FUNCTION Matpar::GetHeatCapacity {{{1*/
double Matpar::GetHeatCapacity(){
	return heatcapacity;
}
/*}}}1*/
/*FUNCTION Matpar::GetLatentHeat {{{1*/
double Matpar::GetLatentHeat(){
	return latentheat;
}
/*}}}1*/
/*FUNCTION Matpar::GetMeltingPoint {{{1*/
double Matpar::GetMeltingPoint(){
	return meltingpoint;
}
/*}}}1*/
/*FUNCTION Matpar::GetReferenceTemperature {{{1*/
double Matpar::GetReferenceTemperature(){
	return referencetemperature;
}
/*}}}1*/
/*FUNCTION Matpar::GetMixedLayerCapacity {{{1*/
double Matpar::GetMixedLayerCapacity(){
	return mixed_layer_capacity;
}
/*}}}1*/
/*FUNCTION Matpar::GetRhoIce {{{1*/
double Matpar::GetRhoIce(){
	
	return rho_ice;
}
/*}}}1*/
/*FUNCTION Matpar::GetRhoWater {{{1*/
double Matpar::GetRhoWater(){
	return rho_water;
}
/*}}}1*/
/*FUNCTION Matpar::GetMuWater {{{1*/
double Matpar::GetMuWater(){
	return mu_water;
}
/*}}}1*/
/*FUNCTION Matpar::GetThermalConductivity {{{1*/
double Matpar::GetThermalConductivity(){
	return thermalconductivity;
}
/*}}}1*/
/*FUNCTION Matpar::GetThermalExchangeVelocity {{{1*/
double Matpar::GetThermalExchangeVelocity(){
	return thermal_exchange_velocity;
}
/*}}}1*/
/*FUNCTION Matpar::GetKn {{{1*/		 
double Matpar::GetKn(){			 
	return kn;		 
}		 
/*}}}1*/			 
/*FUNCTION Matpar::GetHydrologyP {{{1*/			 
double Matpar::GetHydrologyP(){		 
	return hydro_p;			 
}		 
/*}}}1*/			 
/*FUNCTION Matqar::GetHydrologyQ {{{1*/			 
double Matpar::GetHydrologyQ(){		 
	return hydro_q;			 
}		 
/*}}}1*/			 
/*FUNCTION Matpar::GetHydrologyCR {{{1*/		 
double Matpar::GetHydrologyCR(){		 
	return hydro_CR;		 
}		 
/*}}}1*/			 
/*FUNCTION Matpar::GetHydrologyN {{{1*/			 
double Matpar::GetHydrologyN(){		 
	return hydro_n;			 
}		 
/*}}}1*/ 
/*FUNCTION Matpar::TMeltingPoint {{{1*/
double Matpar::TMeltingPoint(double pressure){
	return meltingpoint-beta*pressure;
}
/*}}}1*/
/*FUNCTION Matpar::PureIceEnthalpy{{{1*/
double Matpar::PureIceEnthalpy(double pressure){
	return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
}
/*}}}1*/
/*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{1*/
double Matpar::GetEnthalpyDiffusionParameter(double enthalpy,double pressure){
	if(enthalpy<PureIceEnthalpy(pressure)){
		return thermalconductivity/(rho_ice*heatcapacity);
	}
	else{
		return 0.1*thermalconductivity/(rho_ice*heatcapacity);
	}
}
/*}}}1*/
/*FUNCTION Matpar::EnthalpyToThermal {{{1*/
void Matpar::EnthalpyToThermal(double* ptemperature,double* pwaterfraction,double enthalpy,double pressure){

	/*Ouput*/
	double 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;
}
/*}}}1*/
/*FUNCTION Matpar::ThermalToEnthalpy {{{1*/
void Matpar::ThermalToEnthalpy(double * penthalpy,double temperature,double waterfraction,double pressure){

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

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