/*!\file: love_core.cpp
 * \brief: core of the LOVE numbers solution
 */

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"
void love_core(FemModel* femmodel){

	Vector<IssmDouble> *wg    = NULL;
	Vector<IssmDouble> *dwdtg = NULL;
	IssmDouble          *x    = NULL;
	IssmDouble          *y    = NULL;

	/*love parameters: */
	IssmDouble *frequencies = NULL;
	int         nfreq,dummy;
	int         sh_nmin,sh_nmax;
	IssmDouble  g0,r0,mu0;
	bool        allow_layer_deletion;
	bool        love_kernels;
	int         forcing_type;
	bool        verbosemod = (int)VerboseModule();

	/*parameters: */
	bool save_results;

	if(VerboseSolution()) _printf0_("   computing LOVE numbers\n");

	/*Recover some parameters: */
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);

	/*recover love number parameters: */
	femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
	femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
	femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
	femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
	femmodel->parameters->FindParam(&g0,LoveG0Enum);
	femmodel->parameters->FindParam(&r0,LoveR0Enum);
	femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
	femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
	femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
	femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);

	/*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
	Matlitho* matlitho=NULL;
	for (Object* & object: femmodel->materials->objects){
		Material* material=xDynamicCast<Material*>(object);
		if (material->ObjectEnum()==MatlithoEnum)matlitho=xDynamicCast<Matlitho*>(material);
	}
	_assert_(matlitho);

	/*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
	IssmDouble*  LoveKr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
	IssmDouble*  LoveHr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
	IssmDouble*  LoveLr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
	IssmDouble*  LoveKi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
	IssmDouble*  LoveHi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
	IssmDouble*  LoveLi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));

	/*Initialize love kernels (real and imaginary parts): */
	IssmDouble*  LoveKernelsReal = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
	IssmDouble*  LoveKernelsImag = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);

	/*call the main module: */
	FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //output
			nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
			matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
			matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->isburgers, matlitho->issolid //matlitho inputs
			);

	/*Add love matrices to results:*/
	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
	/*Only when love_kernels is set unity*/
	if (love_kernels==1) {
		femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsRealEnum,LoveKernelsReal,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
		femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
	}

	/*Free ressources:*/
	xDelete<IssmDouble>(frequencies);
	xDelete<IssmDouble>(LoveKr);
	xDelete<IssmDouble>(LoveHr);
	xDelete<IssmDouble>(LoveLr);
	xDelete<IssmDouble>(LoveKi);
	xDelete<IssmDouble>(LoveHi);
	xDelete<IssmDouble>(LoveLi);
	xDelete<IssmDouble>(LoveKernelsReal);
	xDelete<IssmDouble>(LoveKernelsImag);
}
