/*!\file: Friction.cpp
 * \brief: wrapper for all friction parameters
 */ 

#include "./Friction.h"
#include "../shared/shared.h"
#include "../include/typedefs.h"

/*--------------------------------------------------
	NewFriction
  --------------------------------------------------*/

Friction* NewFriction(void)
{
	/* create a new Friction object */

	return (Friction*)xmalloc(sizeof(Friction));
}


/*--------------------------------------------------
	FrictionInit
  --------------------------------------------------*/

int FrictionInit(Friction* friction)
{
	friction->element_type=NULL;
	friction->gravity=UNDEF;
	friction->rho_ice=UNDEF;
	friction->rho_water=UNDEF;
	friction->K=NULL; 
	friction->bed=NULL;
	friction->thickness=NULL;
	friction->velocities=NULL;
	friction->p=UNDEF;
	friction->q=UNDEF;
	
	return 1;
}

/*--------------------------------------------------
	DeleteFriction
  --------------------------------------------------*/

void DeleteFriction(Friction** pfriction)
{
	Friction* friction = *pfriction;

	/*Just  erase element_type: */
	xfree((void**)&friction->element_type);

	/*Erase entire structure: */
	xfree((void**)pfriction);
}


/*--------------------------------------------------
	FrictionGetAlpha2
  --------------------------------------------------*/
void  FrictionGetAlpha2(double* alpha2, Friction* friction){

	/*This routine calculates the basal friction coefficient 
	alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p*
	alpha2 is assumed to be double[3]: */

	/*diverse: */
	int     i;
	int     numgrids=3;
	double  Neff[numgrids];
	double  r,s;
	double  velocity_x[numgrids];
	double  velocity_y[numgrids];
	double  velocity_mag[numgrids];

	//compute r and q coefficients: */
	r=friction->q/friction->p;
	s=1./friction->p;
		
	//From bed and thickness, compute effective pressure when drag is viscous:
	for(i=0;i<numgrids;i++){

		Neff[i]=friction->gravity*(friction->rho_ice*friction->thickness[i]+friction->rho_water*friction->bed[i]);

		/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
		the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
		for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
		replace it by Neff=0 (ie, equival it to an ice shelf)*/
		if (Neff[i]<0)Neff[i]=0;

		//We need the velocity magnitude to evaluate the basal stress:
		velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2
		velocity_y[i]=*(friction->velocities+2*i+1);
		velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
	
		alpha2[i]=pow(friction->K[i],2)*pow(Neff[i],r)*pow(velocity_mag[i],(s-1));
	}
}

/*--------------------------------------------------
	FrictionGetAlphaComplement
  --------------------------------------------------*/
void  FrictionGetAlphaComplement(double* alpha_complement, Friction* friction){

	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
	 * alpha_complement= Neff ^r * vel ^s*/

	/*diverse: */
	int     i;
	int     numgrids=3;
	double  Neff[numgrids];
	double  r,s;
	double  velocity_x[numgrids];
	double  velocity_y[numgrids];
	double  velocity_mag[numgrids];

	//compute r and q coefficients: */
	r=friction->q/friction->p;
	s=1./friction->p;
		
	//From bed and thickness, compute effective pressure when drag is viscous:
	for(i=0;i<numgrids;i++){

		Neff[i]=friction->gravity*(friction->rho_ice*friction->thickness[i]+friction->rho_water*friction->bed[i]);

		/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
		the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
		for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
		replace it by Neff=0 (ie, equival it to an ice shelf)*/
		if (Neff[i]<0)Neff[i]=0;

		//We need the velocity magnitude to evaluate the basal stress:
		velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2
		velocity_y[i]=*(friction->velocities+2*i+1);
		velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
	
		alpha_complement[i]=pow(Neff[i],r)*pow(velocity_mag[i],(s-1));
	}
}
