#include "../includes.h"
#include "../../inc/ics/flattenedmodel.h"
CFlattenedModel::CFlattenedModel()
{
m_name = "FlattenedModel";
m_toomreQ = 0.0;}
void CFlattenedModel::Activate( double toomreQ )
{
m_toomreQ = toomreQ;
}
double CFlattenedModel::PrepCumulativeMass( double r )
{
return 2.0 * M_PI * r * rho( r );
}double CFlattenedModel::dphi_dr(double r)
{
return (m_G * CumulativeMass( r )) / (r * r);}
void CFlattenedModel::RealizeModel()
{ double Epot = 0;
double Ekin = 0; gsl_rng * rnd = gsl_rng_alloc ( gsl_rng_mt19937 );
std::string file = StartBinary();
double mass = m_totalMass / m_count; 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) ) || ( radius < (m_rMax/100.0) ) );
double phi = 0 + gsl_rng_uniform_pos( rnd ) * ( 2 * M_PI + 0.0 );
double px = radius * cos( phi );
double py = radius * sin( phi );
double pz = 0;
double vx, vy, vz;
vx = vy = vz = 0;
double pot = Potential( radius );
Epot += pot;
double psi = -pot;
double v = sqrt( dphi_dr( radius ) * radius );
vx = - v * sin( phi );
vy = v * cos( phi );
vz = 0; if (m_toomreQ != 0.0)
{
double omega = v/radius; double kappa = sqrt(2.0)*omega; double sigmaRad = ( m_toomreQ * 3.36 * m_G * sigma(radius) ) / kappa;
double sigmaRadX = sigmaRad * sin( phi );
double sigmaRadY = sigmaRad * cos( phi );
vx = CSphericalModel::RandomGaussian( vx, sigmaRadX );
vy = CSphericalModel::RandomGaussian( vy, sigmaRadY ); double sigmaTan = ( m_toomreQ * 3.36 * m_G * sigma(radius) ) / (2.0 * omega );
double sigmaTanX = sigmaTan * sin( phi );
double sigmaTanY = sigmaTan * cos( phi );
vx = CSphericalModel::RandomGaussian( vx, sigmaTanX );
vy = CSphericalModel::RandomGaussian( vy, sigmaTanY ); vz = CSphericalModel::RandomGaussian( vz, 0.5 * sigmaTan );
}
Ekin += 0.5 * ( vx * vx + vy * vy + vz * vz ); BinaryWrite( mass );
BinaryWrite( px );
BinaryWrite( py );
BinaryWrite( pz );
BinaryWrite( vx );
BinaryWrite( vy );
BinaryWrite( vz );
}
StopBinary();
xmlNode * pRootNode = StartXml(); StopXml(pRootNode); gsl_rng_free( rnd );
std::cout << "Q = " << m_toomreQ << std::endl;
PrintFinish( Ekin, Epot );
std::cout << "Press enter key to finish...\a";
getchar();
}
double CFlattenedModel::rho( double r )
{
return sigma( r );
}