#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(); 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 ); gsl_integration_workspace_free( m_wspc );
}
void CSphericalModel::Generate()
{
CInitialConditions::Generate();
TabulateCumulativeMass();
TabulatePotential();
Tabulate_rho_phi();
TabulateDistributionFunction();
Tabulate_rm();
RealizeModel();
}double __PrepCumulativeMass( double r, void * obj )
{
return reinterpret_cast<CSphericalModel*>(obj)->PrepCumulativeMass( r );
}double __dphi_dr(double r, void * obj)
{
return reinterpret_cast<CSphericalModel*>(obj)->dphi_dr( r );
}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; 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; double res, abserr;
gsl_function F;
F.function = &__PrepCumulativeMass;
F.params = this; 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;
} m_mr = gsl_spline_alloc(gsl_interp_cspline, m_funcTabs);
gsl_spline_init( m_mr, radii, mass, m_funcTabs );
delete [] radii;
delete [] mass; 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) {
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)
{ double res, abserr;
gsl_function F;
F.function = &__dphi_dr;
F.params = this;
if (r > 0)
{ 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 - 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]; double max = TIMES_RADII * m_rMax; 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;
} m_Phir = gsl_spline_alloc(gsl_interp_cspline, m_funcTabs);
gsl_spline_init( m_Phir, radii, potential, m_funcTabs );
delete [] radii;
delete [] potential; 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]; double max = TIMES_RADII * m_rMax; double dStep = 1.0 / m_funcTabs;
double d = dStep; int start = 0;
bool bSmall = false;
for(unsigned int i = 0; i < m_funcTabs; i++ )
{
double r = log(d) * max;
phi[i] = - Potential(r); if ((phi[i] > 1e-6) && !bSmall)
{
bSmall = true; start = i;
}
rhos[i] = rho(r);
d += dStep;
} 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; if (m_bUsePlot)
{
m_plotOut.StartPlot( "phi", "rho" );
double max = TIMES_RADII * m_rMax; 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();
} if (m_bUsePlot)
{
m_plotOut.StartPlot( "phi", "d_phi-d_rho" );
double max = TIMES_RADII * m_rMax; 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; 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]; double max = TIMES_RADII * m_rMax; double dStep = 1.0 / m_funcTabs;
double d = dStep; int start = 0;
bool bSmall = false;
for(unsigned int i = 0; i < m_funcTabs; i++ )
{
double r = log(d) * max;
m_E = -Potential( r ); 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; if ((pf[i] > 1e-6) && !bSmall)
{
bSmall = true; start = i;
}
d += dStep;
} 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; if (m_bUsePlot)
{
m_plotOut.StartPlot( "undif-phi", "undif-f" );
double max = TIMES_RADII * m_rMax; 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();
} if (m_bUsePlot)
{
m_plotOut.StartPlot( "phi", "f" );
double max = TIMES_RADII * m_rMax; 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]; const double max = TIMES_RADII * m_rMax; double d = 1.0;
const double dStep = 1.0 / m_funcTabs; int start = 0;
bool bSmall = false;
for(unsigned int i = 0; i < m_funcTabs; i++)
{
double r = log(d) * max; double res, abserr;
gsl_function F;
F.function = &__PrepCumulativeMass;
F.params = this; 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;
} m_Rm = gsl_spline_alloc(gsl_interp_akima, m_funcTabs);
gsl_spline_init( m_Rm, M, fr, m_funcTabs );
delete [] fr;
delete [] M; 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()
{ double Epot = 0;
double Ekin = 0; 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; radius = rm( mx );
}
while ( radius > m_rCut );
double theta = acos( -1.0 + gsl_rng_uniform_pos( rnd ) * ( 1.0 + 1.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 );
double fmax = 1.1 * f( psi );
double f0 = 0.0; double f1 = 1.0; double v2;
while ( f1 > f0 )
{ v2 = 1.1 * vmax2; while ( v2 > vmax2 ) {
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 );
f1 = fmax * gsl_rng_uniform_pos( rnd ); }
Ekin += 0.5 * ( v2 ); BinaryWrite( mass );
BinaryWrite( px );
BinaryWrite( py );
BinaryWrite( pz );
BinaryWrite( vx );
BinaryWrite( vy );
BinaryWrite( vz );
}
StopBinary();
xmlNode * pRootNode = StartXml(); StopXml(pRootNode); 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;
}