Commit 4941a7b8 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Now we generate the true plummer distribution.

parent 862d7cb0
...@@ -107,9 +107,9 @@ If you consider the ...@@ -107,9 +107,9 @@ If you consider the
\subsection{Plummer Model} \subsection{Plummer Model}
This is a hard test case in astrophysics problem, and it models a globular cluster of stars, which is highly non uniform. It is called the plummer distribution. To construct such distribution, first we construct a uniform points distribution on the unit sphere. Second, the radius is chosen according to the plummer distribution (double power law in astrophysics). We consider $u$ a random number between 0 and 1, then the associated radius is given by This is a hard test case in astrophysics problem, and it models a globular cluster of stars, which is highly non uniform. It is called the plummer distribution. To construct such distribution, first we construct a uniform points distribution on the unit sphere. Second, the radius is chosen according to the plummer distribution (double power law in astrophysics). We consider $u$ a random number between 0 and 1, then the associated radius is given by
\begin{equation*} \begin{equation*}
r = \sqrt{\frac{u^{2/3}}{u^{2/3}-1}} r = 1.0/\sqrt{u^{-2/3}-1},
\end{equation*} \end{equation*}
and the total mass is one. Then, $m_i = \frac{1}{Npt}$.
\begin{figure}[h] \begin{figure}[h]
\centering \centering
\begin{minipage}{0.45\textwidth}% \begin{minipage}{0.45\textwidth}%
...@@ -140,6 +140,8 @@ The corresponding potential is ...@@ -140,6 +140,8 @@ The corresponding potential is
\begin{equation} \begin{equation}
\Phi_P(r) = - \frac{G M}{\sqrt{r^2+a^2}} \Phi_P(r) = - \frac{G M}{\sqrt{r^2+a^2}}
\end{equation} \end{equation}
In N-body units, $G = M = 1$ and $a = 3\pi/16 \sim 0.589$
\subsection{Diagonal Model} \subsection{Diagonal Model}
%, shape end size=.5cm},decoration={shape start size=.5cm, shape end size=.125cm %, shape end size=.5cm},decoration={shape start size=.5cm, shape end size=.125cm
......
...@@ -97,7 +97,7 @@ namespace Param { ...@@ -97,7 +97,7 @@ namespace Param {
const FParameterNames Size const FParameterNames Size
= {{"-size"}, "Size of the geometry a:b:c - default values are 1.0:1.0:2.0"}; = {{"-size"}, "Size of the geometry a:b:c - default values are 1.0:1.0:2.0"};
const FParameterNames Radius const FParameterNames Radius
= {{"-radius"}, "used to specified the radius of the sphere an dthe plummer distribution or R - default value for R is 2.0"}; = {{"-radius"}, "used to specified the radius of the sphere and the plummer distribution or R - default value for R is 2.0"};
const FParameterNames Charge const FParameterNames Charge
= {{"-charge"}, "generate physical values between -1 and 1 otherwise generate between 0 and 1"}; = {{"-charge"}, "generate physical values between -1 and 1 otherwise generate between 0 and 1"};
const FParameterNames ZM const FParameterNames ZM
...@@ -250,6 +250,10 @@ int main(int argc, char ** argv){ ...@@ -250,6 +250,10 @@ int main(int argc, char ** argv){
std::cout << "End of writing" <<std::endl; std::cout << "End of writing" <<std::endl;
// Generate file for visualization // Generate file for visualization
// if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)){
// std::string outfilename(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.vtp"));
// driverExportData(outfilename, particles , NbPoints,loader.getNbRecordPerline() );
// }
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)) { if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)) {
std::string visufile(FParameters::getStr(argc, argv, FParameterDefinitions::OutputVisuFile.options, "output.vtp")); std::string visufile(FParameters::getStr(argc, argv, FParameterDefinitions::OutputVisuFile.options, "output.vtp"));
driverExportData(visufile, particles , NbPoints); driverExportData(visufile, particles , NbPoints);
......
...@@ -243,12 +243,12 @@ FReal plummerDist(FSize& cpt, const FReal &R) { ...@@ -243,12 +243,12 @@ FReal plummerDist(FSize& cpt, const FReal &R) {
} }
} }
/** /**
* \brief Build N points following the Plummer distribution * \brief Build N points following the Plummer distribution
* *
* First we construct N points uniformly distributed on the unit sphere. Then * First we construct N points uniformly distributed on the unit sphere. Then
* the radius in construct according to the Plummr distribution. * the radius in construct according to the Plummer distribution for
* a constant mass of 1/N
* *
* \tparam FReal Floating point type * \tparam FReal Floating point type
* *
...@@ -259,6 +259,32 @@ FReal plummerDist(FSize& cpt, const FReal &R) { ...@@ -259,6 +259,32 @@ FReal plummerDist(FSize& cpt, const FReal &R) {
template <class FReal> template <class FReal>
void unifRandomPlummer(const FSize N, const FReal R, FReal * points) { void unifRandomPlummer(const FSize N, const FReal R, FReal * points) {
unifRandomPointsOnSphere<FReal>(N, 1, points); unifRandomPointsOnSphere<FReal>(N, 1, points);
FReal mc = 1.0/static_cast<FReal>(N);
for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) {
FReal m = getRandom<FReal>();
FReal r = FMath::Sqrt( 1.0/(FMath::pow(m, -2.0/3.0) - 1.0)) ;
points[j] *= r;
points[j+1] *= r;
points[j+2] *= r;
points[j+3] = mc; // the mass
}
}
/**
* \brief Build N points following the Plummer like distribution
*
* First we construct N points uniformly distributed on the unit sphere. Then
* the radius in construct according to the Plummer like distribution.
*
* \tparam FReal Floating point type
*
* \param N the number of points following the Plummer distribution
* \param R the radius of the sphere that contains all the points
* \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
*/
template <class FReal>
void unifRandomPlummerLike(const FSize N, const FReal R, FReal * points) {
FReal a = 1.0 ;
unifRandomPointsOnSphere<FReal>(N, 1, points);
FReal r, rm = 0.0; FReal r, rm = 0.0;
FSize cpt = 0; FSize cpt = 0;
for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) { for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment