/*************************************************************************************
 *
 *  sphericalmodel.cpp:  suitable for deriving spherical models
 *
 *  Copyright (c) 2007 Jakub Schwarzmeier
 *  University of West Bohemia, Department of General Physics
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or (at
 *  your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful, but
 *  WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 ************************************************************************/

#include "../includes.h"

#include "../../../inc/ics/sphericalmodel.h"
#include "../../../../amon/inc/simulation/bh.h"

CSphericalModel::CSphericalModel(void)
{
	m_name = "SphericalModel";

	m_E = 0;

	m_Rm = NULL;
	m_mr = NULL;
	m_Phir = NULL;

	m_Rhophi = NULL;
	m_Fp = NULL;

	m_accRm = gsl_interp_accel_alloc();
	m_accMr = gsl_interp_accel_alloc();
	m_accPhir = gsl_interp_accel_alloc();

	m_accRhophi = gsl_interp_accel_alloc();
	m_accFp = gsl_interp_accel_alloc();

	//  Allocate work space for integration.
	m_wspc = gsl_integration_workspace_alloc(m_funcTabs); 
}

CSphericalModel::~CSphericalModel(void)
{
	gsl_interp_accel_free (m_accMr);
	if (m_mr != NULL)  gsl_spline_free( m_mr );

	gsl_interp_accel_free (m_accRm);
	if (m_Rm != NULL)  gsl_spline_free( m_Rm );

	gsl_interp_accel_free (m_accPhir);
	if (m_Phir != NULL)  gsl_spline_free( m_Phir );

	gsl_interp_accel_free (m_accRhophi);
	if (m_Rhophi != NULL)  gsl_spline_free( m_Rhophi );

	gsl_interp_accel_free (m_accFp);
	if (m_Fp != NULL)  gsl_spline_free( m_Fp );

	//  Free work space allocated for integration.
	gsl_integration_workspace_free( m_wspc );
}


void CSphericalModel::Generate()
{
	CInitialConditions::Generate();

	TabulateCumulativeMass();
	TabulatePotential();
	Tabulate_rho_phi();
	TabulateDistributionFunction();
	Tabulate_rm();

	RealizeModel();
}

//--  Bridge to close gap between non-OOP GSL routines and OOP-functions.

//  Calles objective PrepCumulativeMass function.
double __PrepCumulativeMass( double r, void * obj )
{
	return reinterpret_cast<CSphericalModel*>(obj)->PrepCumulativeMass( r );
}

//  Calles objective dphi_dr function.
double __dphi_dr(double r, void * obj)
{
	return reinterpret_cast<CSphericalModel*>(obj)->dphi_dr( r );
} 

//  Calles objective energy function.
double __Energy(double r, void * obj)
{
	return reinterpret_cast<CSphericalModel*>(obj)->energy( r );
}

// ***********************************************************************

double CSphericalModel::PrepCumulativeMass( double r )
{
	return 4 * M_PI * r * r * rho( r );
}

void CSphericalModel::TabulateCumulativeMass()
{
	double * radii = new double[m_funcTabs];
	double * mass = new double[m_funcTabs];

	const double max = TIMES_RADII * m_rMax;	//  out to 10x radii
	double d = 1.0;
	const double dStep = 1.0 / m_funcTabs;
	for(unsigned int i = 0; i < m_funcTabs; i++)
	{ 
		double r = log(d) * max;

		//  Prepare variables for GSL integration.
		double res, abserr; 
		gsl_function F;
		F.function = &__PrepCumulativeMass;
		F.params = this;

		//  Integrate CumulativeMass.
		gsl_integration_qag(&F, 0.0, r, 0.0, 1.0e-8, m_funcTabs, 
			GSL_INTEG_GAUSS31, m_wspc, &res, &abserr); 
		
		radii[i] = r;
		mass[i] = res;

		d -= dStep;
	}

	//  Construct spline m(r).
	m_mr = gsl_spline_alloc(gsl_interp_cspline, m_funcTabs); 
	gsl_spline_init( m_mr, radii, mass, m_funcTabs );

	delete [] radii;
	delete [] mass;

	//  Plot m(r) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "r", "m" );

		double plotStep = TIMES_RADII * m_rMax / m_plotTabs;
		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = i * plotStep;

			double mr = CumulativeMass( r );
			m_plotOut.Plot( r, mr );

		}

		m_plotOut.StopPlot();
	}
}

double CSphericalModel::CumulativeMass( double r )
{
	return gsl_spline_eval (m_mr, r, m_accMr);
}

// ***********************************************************************

double CSphericalModel::dphi_dr(double r)
{
	double val;
	if (r == 0.0)		//  avoid 1/0
	{
		val = 0.0;
	}
	else
	{
		val = (m_G/(r*r)) * CumulativeMass( r );
		if (val > 1.0) val = 1.0;
		if (val < 0.0) val = 0.0;
	}

	return val;
}

double CSphericalModel::PrePotential(double r)
{
	//  Prepare variables for GSL integration.
	double res, abserr; 
	gsl_function F;
	F.function = &__dphi_dr;
	F.params = this;

	if (r > 0)
	{
		//  Perform integrations of dphi_dr that gives us phi(r).
		gsl_integration_qag(&F, 0, r, 0.0, 1.0e-4,
			m_funcTabs, GSL_INTEG_GAUSS31, m_wspc, &res, &abserr); 
	}
	else
	{
		res = 0;
	}

	double pot = 1.0/*central potential*/ - res;

	if (pot < 0.0)	pot = 0.0;

	return -m_G * pot;
}

void CSphericalModel::TabulatePotential()
{
	double * radii = new double[m_funcTabs];
	double * potential = new double[m_funcTabs];

	//  Tabulate phi(r) in logarithmical scale.
	double max = TIMES_RADII * m_rMax;	//  out to 10x radii
	double d = 1.0;
	double dStep = 1.0 / m_funcTabs;
	for(unsigned int i = 0; i < m_funcTabs; i++)
	{ 
		double r = log(d) * max;

		radii[i] = r;
		potential[i] = PrePotential(r);

		d -= dStep;
	}

	//  Construct spline Phi(r)
	m_Phir  = gsl_spline_alloc(gsl_interp_cspline, m_funcTabs); 
	gsl_spline_init( m_Phir, radii, potential, m_funcTabs );

	delete [] radii;
	delete [] potential;

	//  Plot phi(r) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "r", "phi" );

		double plotStep = TIMES_RADII * m_rMax / m_plotTabs;
		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = i * plotStep;

			double phix = Potential( r );
			m_plotOut.Plot( r, phix );
		}

		m_plotOut.StopPlot();
	}

	return;
}

double CSphericalModel::Potential( double r )
{
	return gsl_spline_eval( m_Phir, r, m_accPhir);
}

// ***********************************************************************

void CSphericalModel::Tabulate_rho_phi()
{
	double * rhos = new double[m_funcTabs];
	double * phi = new double[m_funcTabs];

	//  Prepare
	double max = TIMES_RADII * m_rMax;		//  out to 10x radii
	double dStep = 1.0 / m_funcTabs;
	double d = dStep;

	//  Avoid too small numbers variables.
	int start = 0;
	bool bSmall = false;

	//  Tabulate rho( phi ).

	for(unsigned int i = 0; i < m_funcTabs; i++ )
	{ 
		double r = log(d) * max;

		phi[i] = - Potential(r);
		//  Avoid too small numbers in order to prevent difficulties 
		//  with spline interpolation.
		if ((phi[i] > 1e-6) && !bSmall) 	
		{ 
			bSmall = true; start = i; 
		}

		rhos[i] = rho(r);

		d += dStep;
	}



	//  Construct spline rho( phi )
	m_Rhophi  = gsl_spline_alloc(gsl_interp_cspline, m_funcTabs-start); 
	gsl_spline_init( m_Rhophi, &phi[start], &rhos[start], m_funcTabs-start );

	delete [] phi;
	delete [] rhos;

	//  Plot rho( phi ) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "phi", "rho" );

		double max = TIMES_RADII * m_rMax;		//  out to 10x radii
		double dStep = 1.0 / m_plotTabs;
		double d = dStep;

		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = log(d) * max;

			double phix = Potential( r );
			double rho = gsl_spline_eval (m_Rhophi, -phix, m_accRhophi);
			m_plotOut.Plot( phix, rho );

			d += dStep;
		}

		m_plotOut.StopPlot();
	}

	//  Plot d_rho/d_rho( phi ) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "phi", "d_phi-d_rho" );

		double max = TIMES_RADII * m_rMax;		//  out to 10x radii
		double dStep = 1.0 / m_plotTabs;
		double d = dStep;

		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = log(d) * max;

			double phix = Potential( r );
			double der = d_rho_d_phi ( -phix );
			m_plotOut.Plot( phix, der );

			d += dStep;
		}

		m_plotOut.StopPlot();
	}
}

double CSphericalModel::d_rho_d_phi( double p )
{
	return gsl_spline_eval_deriv (m_Rhophi, p, m_accRhophi);
}

// ***********************************************************************

double CSphericalModel::energy( double p )
{
    double subtract = m_E - p;	
	//  division with small numbers introduces errors
	if (subtract < 1e-6)  return 0;

    double res = d_rho_d_phi( p ) / sqrt( subtract );
    return res;
}

void CSphericalModel::TabulateDistributionFunction()
{
	double * pf = new double[m_funcTabs];
	double * ff = new double[m_funcTabs];

	//  Prepare
	double max = TIMES_RADII * m_rMax;		//  out to 10x radii
	double dStep = 1.0 / m_funcTabs;
	double d = dStep;

	//  Avoid too small numbers variables.
	int start = 0;
	bool bSmall = false;

	for(unsigned int i = 0; i < m_funcTabs; i++ )
	{ 
		double r = log(d) * max;

		m_E = -Potential( r );

		//  Prepare variables for GSL integration.
		double res, abserr; 
		gsl_function F;
		F.function = &__Energy;
		F.params = this;

		gsl_integration_qag(&F, 0, m_E, 0.0, 1.0e-4,
			   m_funcTabs, GSL_INTEG_GAUSS31, m_wspc, &res, &abserr); 

		pf[i] = m_E;
		ff[i] = res;

		//  Avoid too small numbers in order to prevent difficulties 
		//  with spline interpolation.
		if ((pf[i] > 1e-6) && !bSmall) 	
		{ 
			bSmall = true; start = i; 
		}

		d += dStep;
	}

	//  Construct spline f( p )
	m_Fp  = gsl_spline_alloc(gsl_interp_cspline, m_funcTabs-start); 
	gsl_spline_init( m_Fp, &pf[start], &ff[start], m_funcTabs-start );

	delete [] ff;
	delete [] pf;

	//  Plot undifferentiated f( phi ) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "undif-phi", "undif-f" );

		double max = TIMES_RADII * m_rMax;		//  out to 10x radii
		double dStep = 1.0 / m_plotTabs;
		double d = dStep;

		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = log(d) * max;

			double phix = Potential( r );
			double f = gsl_spline_eval (m_Fp, -phix, m_accFp);
			m_plotOut.Plot( phix, f );

			d += dStep;
		}

		m_plotOut.StopPlot();
	}

	//  Plot final f( phi ) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "phi", "f" );

		double max = TIMES_RADII * m_rMax;		//  out to 10x radii
		double dStep = 1.0 / m_plotTabs;
		double d = dStep;

		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = log(d) * max;

			double phi = Potential(r);
			double fff = f( -phi );

			m_plotOut.Plot( phi, fff );

			d += dStep;
		}

		m_plotOut.StopPlot();
	}

}

double CSphericalModel::f( double p )
{
	double dummy;

	if (p < 0.09)
	{
		dummy = 1.0e-7;
	}
	else
	{
		dummy = gsl_spline_eval_deriv (m_Fp, p, m_accFp);
		dummy /= ( sqrt( 8.0 ) * M_PI * M_PI );
	}

	if (dummy < 0) dummy = 1.0e-7;

	return dummy;
}

// ***********************************************************************


void CSphericalModel::Tabulate_rm()
{
	//

	double * fr = new double[m_funcTabs];
	double * M = new double[m_funcTabs];

	//  Inverse m(r) to receive r(m).
	const double max = TIMES_RADII * m_rMax;	//  out to 10x radii
	double d = 1.0;
	const double dStep = 1.0 / m_funcTabs;

	//  Avoid too small deviations in numbers.
	int start = 0;
	bool bSmall = false;

	for(unsigned int i = 0; i < m_funcTabs; i++)
	{ 
		double r = log(d) * max;

		//  Prepare variables for GSL integration.
		double res, abserr; 
		gsl_function F;
		F.function = &__PrepCumulativeMass;
		F.params = this;

		//  Integrate CumulativeMass.
		gsl_integration_qag(&F, 0.0, r, 0.0, 1.0e-8, m_funcTabs, 
			GSL_INTEG_GAUSS31, m_wspc, &res, &abserr); 
		
		M[i] = res;
		fr[i] = r;

		d -= dStep;
	}

	//  Construct spline r( m )
	m_Rm  = gsl_spline_alloc(gsl_interp_akima, m_funcTabs);
	gsl_spline_init( m_Rm, M, fr, m_funcTabs );

	delete [] fr;
	delete [] M;

	//  Plot m(r) if requested.
	if (m_bUsePlot)
	{
		m_plotOut.StartPlot( "m", "r" );

		double plotStep = TIMES_RADII * m_rMax / m_plotTabs;
		for (unsigned int i = 0; i < m_plotTabs; i++ )
		{
			double r = i * plotStep;
			if (r > m_rMax)	r = m_rMax;

			double mm = CumulativeMass ( r );
			double rr = rm ( mm );
			m_plotOut.Plot( mm, rr );

		}

		m_plotOut.StopPlot();
	}

    return;
}

double CSphericalModel::rm( double m )
{
	return gsl_spline_eval (m_Rm, m, m_accRm);
}

// ***********************************************************************

void CSphericalModel::RealizeModel()
{
	//  test energies
    double Epot = 0;
    double Ekin = 0;

	//  prepare uniform (normal, Gaussian) random number generator
	gsl_rng * rnd = gsl_rng_alloc ( gsl_rng_mt19937 );

        double mass = m_totalMass / m_count;

	std::string file = StartBinary();

    for ( unsigned int i = 0; i < m_count; i++ )
    {
		double mx, radius;
		do
		{
			mx = gsl_rng_uniform_pos( rnd ) * m_totalMass; //  CumulativeMass( TIMES_RADII * m_rMax );
			radius = rm( mx );
		}
		while ( radius > m_rCut );


		double theta = acos( -1.0 + gsl_rng_uniform_pos( rnd ) * ( 1.0 + 1.0 ) );
//		double theta = 0 + gsl_rng_uniform_pos( rnd ) * ( 2 * M_PI + 0.0 );
		double phi = 0 + gsl_rng_uniform_pos( rnd ) * ( 2 * M_PI + 0.0 );

        double px = radius * sin( theta ) * cos( phi );
        double py = radius * sin( theta ) * sin( phi );
        double pz = radius * cos( theta );

        double vx, vy, vz;
        vx = vy = vz = 0;

        double pot = Potential( radius );
        Epot += pot;
		double psi = -pot;

        double vmax2 = 2.0 * psi;
        double vmax = sqrt( vmax2 );	//  v_esc = sqrt(2*Phi)

        double fmax = 1.1 * f( psi );	//  1.1 get over the graf for rejection
//		double fmax2 = 1.1 * 1.56e-1* pow(fabs(pot),7/2);

        double f0 = 0.0; double f1 = 1.0;	//  probabilities, just some dummy initial values
											//  to get into the while loop
		double v2;
        while ( f1 > f0 )
        {
            // pick a velocity
            v2 = 1.1 * vmax2;	//  1.1 get into while loop
            while ( v2 > vmax2 )		//  rejection technique
										//  leave loop when v_e^2 <= 2 * Phi
            {
				vx = vmax * ( -1.0 + gsl_rng_uniform_pos( rnd ) * ( 1.0 + 1.0 ) );
				vy = vmax * ( -1.0 + gsl_rng_uniform_pos( rnd ) * ( 1.0 + 1.0 ) );
				vz = vmax * ( -1.0 + gsl_rng_uniform_pos( rnd ) * ( 1.0 + 1.0 ) );
                v2 = vx * vx + vy * vy + vz * vz;

            }

            f0 = f( psi - 0.5 * v2 );
//			double f0_2 = 1.56e-1* pow( (-pot - 0.5 * v2),7/2);

			f1 = fmax * gsl_rng_uniform_pos( rnd );	//  y-maximum for rejection on the graph
        }

		Ekin += 0.5 * ( v2 );

		//  write ICs binary file
		BinaryWrite( mass );

		BinaryWrite( px );
		BinaryWrite( py );
		BinaryWrite( pz );

		BinaryWrite( vx );
		BinaryWrite( vy );
		BinaryWrite( vz );
    }

	StopBinary();

	xmlNode * pRootNode = StartXml();
	//  insert something into XML tree if needed
	StopXml(pRootNode);

	//  free random number generator object
	gsl_rng_free( rnd );

	PrintFinish( Ekin, Epot );

	std::cout << "Press enter key to finish...\a";
	getchar();
}


// ***********************************************************************

void CSphericalModel::SetPotential( gsl_spline * phir )
{
	if (m_Phir != NULL)  gsl_spline_free( m_Phir );
	m_Phir = phir;
}

// ***********************************************************************

double CSphericalModel::log(double x)
{
	double res;
	if (x < 0.5)	
		res = pow(10.0, -x);
	else
        res = -log10(x);

	return res;
}

double CSphericalModel::RandomGaussian(double mean, double stddev)
{
	double radiussq;                    
	double value1, value2;

	do
	{
		value1 = 2.0 * (double)::rand() / RAND_MAX - 1.0;
		value2 = 2.0 * (double)::rand() / RAND_MAX - 1.0;
		radiussq = value1 * value1 + value2 * value2;

	} while ((radiussq >= 1.0) || ((radiussq - 0.000001) < 0.0));

	double factor = ::sqrt(-2.0 * ::log(radiussq) / radiussq);

	return value1 * factor * stddev + mean;
}