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

/*header files: */
/*{{{*/
#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"
/*}}}*/

/*TransientParam constructors and destructor*/
/*FUNCTION TransientParam::TransientParam(){{{*/
TransientParam::TransientParam(){
	return;
}
/*}}}*/
/*FUNCTION TransientParam::TransientParam(int enum_type,IssmDoubleMat value){{{*/
TransientParam::TransientParam(int in_enum_type,IssmDouble* in_values,IssmDouble* in_time,int in_N){

	_assert_(in_values && in_time);

	enum_type=in_enum_type;
	N=in_N;

	values=xNew<IssmDouble>(N);
	xMemCpy<IssmDouble>(values,in_values,N);

	timesteps=xNew<IssmDouble>(N);
	xMemCpy<IssmDouble>(timesteps,in_time,N);
}
/*}}}*/
/*FUNCTION TransientParam::~TransientParam(){{{*/
TransientParam::~TransientParam(){
	xDelete<IssmDouble>(values);
	xDelete<IssmDouble>(timesteps);
}
/*}}}*/

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

	_printf_("TransientParam:\n");
	_printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
	_printf_("   size: " << this->N << "\n");

}
/*}}}*/
/*FUNCTION TransientParam::DeepEcho{{{*/
void TransientParam::DeepEcho(void){

	_printf_("TransientParam:\n");
	_printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
	_printf_("   size: " << this->N << "\n");
	for(int i=0;i<this->N;i++){
		_printf_(   "time: " << this->timesteps[i] << " value: " << this->values[i] << "\n");
	}
}
/*}}}*/
/*FUNCTION TransientParam::Id{{{*/
int    TransientParam::Id(void){ return -1; }
/*}}}*/
/*FUNCTION TransientParam::ObjectEnum{{{*/
int TransientParam::ObjectEnum(void){

	return TransientParamEnum;

}
/*}}}*/
/*FUNCTION TransientParam::copy{{{*/
Object* TransientParam::copy() {

	return new TransientParam(this->enum_type,this->values,this->timesteps,this->N);

}
/*}}}*/

/*TransientParam virtual functions definitions: */
/*FUNCTION TransientParam::GetParameterValue(IssmDouble* pdouble,IssmDouble time){{{*/
void  TransientParam::GetParameterValue(IssmDouble* pdouble,IssmDouble time){

	IssmDouble output;
	bool   found;

	/*Ok, we have the time, go through the timesteps, and figure out which interval we 
	 *fall within. Then interpolate the values on this interval: */
	if(time<this->timesteps[0]){
		/*get values for the first time: */
		output=this->values[0];
		found=true;
	}
	else if(time>this->timesteps[this->N-1]){
		/*get values for the last time: */
		output=this->values[this->N-1];
		found=true;
	}
	else{
		/*Find which interval we fall within: */
		for(int i=0;i<this->N;i++){
			if(time==this->timesteps[i]){
				/*We are right on one step time: */
				output=this->values[i];
				found=true;
				break; //we are done with the time interpolation.
			}
			else{
				if(this->timesteps[i]<time && time<this->timesteps[i+1]){
					/*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
					IssmDouble deltat=this->timesteps[i+1]-this->timesteps[i];
					IssmDouble alpha=(time-this->timesteps[i])/deltat;
					output=(1.0-alpha)*this->values[i] + alpha*this->values[i+1];
					found=true;
					break;
				}
				else continue; //keep looking on the next interval
			}
		}
	}
	if(!found)_error_("did not find time interval on which to interpolate values");
	*pdouble=output;
}
/*}}}*/
