FGenerateDistribution.hpp

#ifndef FGENERATEDISTRIBUTION_HPP
#define FGENERATEDISTRIBUTION_HPP


#include <cstdlib>
#include <ctime>
#include <iostream>
#include <fstream>
#include <algorithm>
//
#include "Utils/FMath.hpp"
#include "Utils/FParameters.hpp"

/**  return a random number between 0 and 1 */
//double getRandom() {
//	return drand48();
//} ;
void initRandom() {
	srand( static_cast<unsigned int>(time(0))) ;
} ;
FReal getRandom() {
	//	return static_cast<FReal>(rdand48());
	return static_cast<FReal>(rand()/FReal(RAND_MAX));
} ;
//!  \fn   unifRandonPointsOnUnitSphere(const int N , FReal * points)

//! \brief Generate N points uniformly distributed on the unit sphere

//!
//! \param N the number of points uniformly randomly sample on the unit sphere
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//! \example  generateDistributions.cpp
void unifRandonPointsOnUnitSphere(const int N , FReal * points) {
	FReal u, v, theta, phi, sinPhi ;
	//
	initRandom() ;
	int j = 0 ;
	for (int i = 0 ; i< N ; ++i, j+=4)  {
		//
		u = getRandom() ;  v = getRandom() ;
		theta  = FMath::FTwoPi*u ;
		phi     = FMath::ACos(2*v-1);
		sinPhi = FMath::Sin(phi);
		//
		points[j]	  =	FMath::Cos(theta)*sinPhi ;
		points[j+1] =	FMath::Sin(theta)*sinPhi ;
		points[j+2] =	2*v-1 ;
		//
	}
};
//!  \fn  nonunifRandonPointsOnElipsoid(const int N , const FReal &a, const FReal &b, const FReal &c, FReal * points)

//! \brief  Generate N points non uniformly distributed on the ellipsoid of  aspect ratio a:b:c

//!
//! \param N the number of points
//! \param a  the x  semi-axe length
//! \param b  the y  semi-axe length
//! \param c  the z  semi-axe length
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
void nonunifRandonPointsOnElipsoid(const int N , const FReal &a, const FReal &b, const FReal &c, FReal * points) {
	//
	FReal u, v , cosu ;
	int j =0 ;
	for (int i = 0 ; i< N ; ++i, j+=4)  {
		u = getRandom() ;  v = getRandom() ;
		u  = FMath::FPi*u - FMath::FPiDiv2;   v   = FMath::FTwoPi*v - FMath::FPi;
		cosu = FMath::Cos(u) ;
		points[j]	   = a*cosu*FMath::Cos(v)  ;
		points[j+1]  = b*cosu*FMath::Sin(v)  ;
		points[j+2]  = c*FMath::Sin(u)  ;
	}
};
//!  \fn  nonunifRandonPointsOnElipsoid(const int N , const FReal &a, const FReal &c, FReal * points)

//! \brief  Generate N points uniformly distributed on the ellipsoid of  aspect ratio a:a:c

//!
//! \param N the number of points
//! \param a  the x  semi-axe length
//! \param c  the z  semi-axe length
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
void unifRandonPointsOnProlate(const int N , const FReal &a, const FReal &c, FReal * points){
	//
	FReal u, w,v ,ksi ;
	FReal e = (a*a*a*a)/(c*c*c*c) ;
	bool isgood = false;
	int j =0 , cpt =0 ;
	//
	for (int i = 0 ; i< N ; ++i, j+=4)  {
		// Select a random point on the prolate
		do {
			cpt++	;
			u = getRandom() ;  v = getRandom() ;
			u  = 2.0*u - 1.0;   v   = FMath::FTwoPi*v;
			w =FMath::Sqrt(1-u*u) ;
			points[j]	   = a*w*FMath::Cos(v)  ;
			points[j+1]  = a*w*FMath::Sin(v)  ;
			points[j+2]  = c*u ;
			// Accept the position ?
			ksi = a*getRandom()  ;
			//			std::cout << "Gradf  "<<  points[j]*points[j] + points[j+1] *points[j+1]  +e*points[j+2] *points[j+2]  << std::endl;
			isgood = (points[j]*points[j] + points[j+1] *points[j+1]  +e*points[j+2] *points[j+2]  < ksi*ksi );
		} while (isgood);
	}
	std::cout.precision(4);
	std::cout << "Total tested points: "<< cpt << " % of rejected points: "<<100*static_cast<FReal>(cpt-N)/cpt << " %" <<std::endl;

} ;

//!  \fn  unifRandonPointsOnSphere(const int N , const FReal R, FReal * points)

//! \brief Generate N points uniformly distributed on the sphere of radius R

//!
//! \param N the number of points uniformly randomly sample on the sphere
//! \param R the radius of the sphere
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
void unifRandonPointsOnSphere(const int N , const FReal R, FReal * points) {
	//
	unifRandonPointsOnUnitSphere(N , points) ;
	int j =0 ;
	for (int i = 0 ; i< N ; ++i, j+=4)  {
		points[j]	   *= R ;
		points[j+1]  *= R ;
		points[j+2]  *= R ;
	}
};
//!  \fn FReal plummerDist(int & cpt, const FReal &R)

//! \brief   Radial Plummer distribution

//!
//! \param cpt : counter to know how many random selections we need to obtain a radius less than R
//! \param R    : Radius of the sphere that contains the particles
//! @return Return the radius according to the Plummer distribution either double type or float type
//!
FReal plummerDist(int & cpt, const FReal &R) {
	//
	FReal radius ,u ;
	do  {
		//radius =  getRandom() ;
		u        = FMath::pow (getRandom() , 2.0/3.0) ;
		radius = FMath::Sqrt (u/(1.0-u));
		cpt++;
		if(radius  <=R){
			//			std::cout << radius << "    "  <<std::endl;
			return static_cast<FReal>(radius);
		}
	} while (true);
}
//! \fn void unifRandonPlummer(const int N , const FReal R, const FReal M, FReal * points)

//! \brief  Build N points following the Plummer distribution

//! First we construct N points uniformly distributed on the unit sphere. Then the radius in construct according to the Plummr distribution.
//!
//! \param N the number of points following the Plummer distribution
//! \param R the radius of the sphere that contains all the points
//! \param M the total mass of all the particles inside the Sphere or radius R
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
void unifRandonPlummer(const int N , const FReal R, const FReal M, FReal * points) {
	//
	unifRandonPointsOnUnitSphere(N , points) ;
	//
	FReal r , rm= 0.0;
	//	FReal Coeff =  3.0*M/(4.0*FMath::FPi*R*R*R) ;
	//am1 = 0 ;//1/FMath::pow(1+R*R,2.5);
	int cpt = 0 ;
	for (int i = 0,j=0 ; i< N ; ++i, j+=4)  {
		// u \in []
		r = plummerDist(cpt,R) ;
		rm = std::max(rm, r);
		points[j]	   *= r ;
		points[j+1]  *= r ;
		points[j+2]  *= r ;
	}

	std::cout << "Total tested points: "<< cpt << " % of rejected points: "
			<<100*static_cast<FReal>(cpt-N)/cpt << " %" <<std::endl;

} ;
//! \fn void exportCVS(std::ofstream& file, const int N, const FReal * particles )

//! \brief  Export particles in CVS Format
//!
//! Export particles in CVS Format as follow
//!      x ,  y  , z , physicalValue
//! It is useful to plot the distribution with paraView
//!
void exportCVS(std::ofstream& file, const int N, const FReal * particles ){
	int j = 0;
	file << " x ,  y , z, q " <<std::endl;
	for(int i = 0 ; i< N; ++i, j+=4){
		file <<    particles[j]    << " , "    <<   particles[j+1]    << " , "   <<   particles[j+2]    << " , "   <<   particles[j+3]   <<std::endl;
	}
}
//
void exportCOSMOS(std::ofstream& file, const int N, const FReal * particles ){
	int j = 0;
	file << " x ,  y , z, q " <<std::endl;
	for(int i = 0 ; i< N; ++i, j+=4){
		file <<    particles[j]    << "  "    <<   particles[j+1]    << "  "   <<   particles[j+2]    << "  0.0 0.0 0.0  "   <<   particles[j+3]   <<"  " << i << std::endl;
	}
}
//
void exportVTK(std::ofstream& VTKfile, const int N, const FReal * particles ){
	int j = 0;
	//---------------------------
	// print generic information
	//---------------------------
	VTKfile << "# vtk DataFile Version 3.0" << "\n";
	VTKfile << "#  Generated bt exportVTK" << "\n";

	VTKfile << "ASCII" << "\n";
	VTKfile << "DATASET POLYDATA" << "\n";
	//
	//---------------------------------
	// print nodes ordered by their TAG
	//---------------------------------
	VTKfile << "POINTS " << N << "  float" << "\n";
	//
	for(int i = 0 ; i< N; ++i, j+=4){
		VTKfile <<    particles[j]    << "  "    <<   particles[j+1]    << "   "   <<   particles[j+2]      <<std::endl;
	}
	// ------------------------------------------
	VTKfile << "\n";
	VTKfile << "VERTICES  " <<  N << " " << 2*N << "\n";
	for(int i = 0 ; i< N; ++i, j+=4){
		VTKfile <<    "  1 "    << " "    <<i<<std::endl;
	}
	VTKfile << "POINT_DATA  " <<  N << "\n";
	VTKfile << "SCALARS PhysicalValue  float 1" << "\n"
			<< "LOOKUP_TABLE default" << "\n" ;
	j = 0 ;
	for(int i = 0 ; i< N; ++i, j+=4){
		VTKfile <<    particles[j+3]    << " "    <<std::endl;
	}
	VTKfile << "\n";
};

void exportVTKxml(std::ofstream& VTKfile, const int N, const FReal * particles ){
	int j = 0;

	VTKfile << "<?xml version=\"1.0\"?>" <<std::endl
			<< "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<<std::endl
			<< "<PolyData>"<<std::endl
			<< "<Piece NumberOfPoints=\" " << N << " \"  NumberOfVerts=\" "<<N <<" \" NumberOfLines=\" 0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">"<<std::endl
			<< "<Points>"<<std::endl
			<< "<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"> "<<std::endl ;
	j = 0 ;
	for(int i = 0 ; i< N; ++i, j+=4){
		VTKfile <<    particles[j]    << "  "    <<   particles[j+1]    << "   "   <<   particles[j+2]      << "   "   ;
	}
	VTKfile <<std::endl<< "</DataArray> "<<std::endl
			<< "</Points> "<<std::endl
			<< "<PointData Scalars=\"PhysicalValue\" > "<<std::endl
			<< "<DataArray type=\"Float64\" Name=\"PhysicalValue\"  format=\"ascii\">"<<std::endl ;
	j = 0 ;
	for(int i = 0 ; i< N; ++i, j+=4){
		VTKfile <<    particles[j+3]    << " "   ;
	}
	VTKfile <<std::endl << "</DataArray>"<<std::endl
			<< "	</PointData>"<<std::endl
			<< "	<CellData>"<<" </CellData>"<<std::endl
			<< "	<Verts>"<<std::endl
			<< "	<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">"<<std::endl ;
	for(int i = 0 ; i< N; ++i){
		VTKfile <<   i   << " "   ;
	}
	VTKfile<<std::endl << "</DataArray>" <<std::endl
			<< "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">"<<std::endl ;
	for(int i = 1 ; i< N+1; ++i){
		VTKfile <<   i   << " "   ;
	}
	VTKfile<<std::endl  << "</DataArray>"<<std::endl
			<< "	</Verts>"<<std::endl
			<< "<Lines></Lines>"<<std::endl
			<< "<Strips></Strips>"<<std::endl
			<< "<Polys></Polys>"<<std::endl
			<< "</Piece>"<<std::endl
			<< "</PolyData>"<<std::endl
			<< "</VTKFile>"<<std::endl;
} ;
//

#endif