FGenerateDistribution.hpp 9.31 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// ===================================================================================
// Copyright ScalFmm 2014 INRIA
//
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// 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 and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
16
17
18
#ifndef FGENERATEDISTRIBUTION_HPP
#define FGENERATEDISTRIBUTION_HPP

COULAUD Olivier's avatar
COULAUD Olivier committed
19
// @author O. Coulaud
COULAUD Olivier's avatar
COULAUD Olivier committed
20

21
22
23
24
25
26
27
28
29
30
#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 */
31

32
void initRandom() {
33
	srand48( static_cast<long int>(time(nullptr))) ;
34
} ;
35
template <class FReal>
36
FReal getRandom() {
COULAUD Olivier's avatar
COULAUD Olivier committed
37
38
	return static_cast<FReal>(drand48());
	//return static_cast<FReal>(rand()/FReal(RAND_MAX));
39
} ;
COULAUD Olivier's avatar
COULAUD Olivier committed
40
41
42
43
44
45
46
47
//!  \fn   unifRandonPointsOnUnitCube(const int N , FReal * points)

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

//!
//! \param N the number of points uniformly randomly sample on the unit cube
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//! \example  generateDistributions.cpp
48
template <class FReal>
49
void unifRandonPointsOnUnitCube(const FSize N , FReal * points) {
COULAUD Olivier's avatar
COULAUD Olivier committed
50
51
52
	//
	initRandom() ;
	int j = 0;
53
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
COULAUD Olivier's avatar
COULAUD Olivier committed
54
		//
55
56
57
        points[j]	  =	getRandom<FReal>()  ;
        points[j+1] =	getRandom<FReal>()  ;
        points[j+2] =	getRandom<FReal>()  ;
COULAUD Olivier's avatar
COULAUD Olivier committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
		//
	}
};
//!  \fn   unifRandonPointsOnCube(const int N , FReal * points)

//! \brief Generate N points uniformly distributed on the cube of length R

//!
//! \param N the number of points uniformly randomly sample on the unit cube
//! \param Lx the the X-length of the  cube
//! \param Ly the the Y-length of the  cube
//! \param Lz the the Z-length of the  cube
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//! \example  generateDistributions.cpp
72
template <class FReal>
73
void unifRandonPointsOnCube(const FSize N , const FReal& Lx,  const FReal &Ly,  const FReal& Lz, FReal * points) {
COULAUD Olivier's avatar
COULAUD Olivier committed
74
75
	//
	unifRandonPointsOnUnitCube(N , points) ;
76
77
    FSize j =0 ;
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
COULAUD Olivier's avatar
COULAUD Olivier committed
78
79
80
81
82
		points[j]	   *= Lx ;
		points[j+1]  *= Ly ;
		points[j+2]  *= Lz ;
	}
};
83
84
85
86
87
88
89
90
//!  \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
91
template <class FReal>
92
void unifRandonPointsOnUnitSphere(const FSize N , FReal * points) {
93
94
95
	FReal u, v, theta, phi, sinPhi ;
	//
	initRandom() ;
96
97
    FSize j = 0 ;
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
98
		//
99
100
        u = getRandom<FReal>() ;  v = getRandom<FReal>() ;
        theta  = FMath::FTwoPi<FReal>()*u ;
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
		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....
//!
121
template <class FReal>
122
void nonunifRandonPointsOnElipsoid(const FSize N , const FReal &a, const FReal &b, const FReal &c, FReal * points) {
123
124
	//
	FReal u, v , cosu ;
125
126
    FSize j =0 ;
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
127
128
        u = getRandom<FReal>() ;  v = getRandom<FReal>() ;
        u  = FMath::FPi<FReal>()*u - FMath::FPiDiv2<FReal>();   v   = FMath::FTwoPi<FReal>()*v - FMath::FPi<FReal>();
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
		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....
//!
145
template <class FReal>
146
void unifRandonPointsOnProlate(const FSize N , const FReal &a, const FReal &c, FReal * points){
147
148
149
150
	//
	FReal u, w,v ,ksi ;
	FReal e = (a*a*a*a)/(c*c*c*c) ;
	bool isgood = false;
151
    FSize j =0 , cpt =0 ;
152
	//
153
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
154
155
156
		// Select a random point on the prolate
		do {
			cpt++	;
157
158
            u = getRandom<FReal>() ;  v = getRandom<FReal>() ;
            u  = 2.0*u - 1.0;   v   = FMath::FTwoPi<FReal>()*v;
159
160
161
162
163
			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 ?
164
            ksi = a*getRandom<FReal>()  ;
165
166
167
168
169
			//			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);
170
    std::cout << "Total tested points: "<< cpt << " % of rejected points: "<<100*static_cast<FReal>(cpt-N)/static_cast<FReal>(cpt) << " %" <<std::endl;
171
172
173

} ;

174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
//!  \fn  unifRandonPointsOnHyperPara(const int N , const FReal &a, const FReal &b, const FReal &c, FReal * points)

//! \brief  Generate N points uniformly distributed on the hyperbolic paraboloid 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....
//!
template <class FReal>
void unifRandonPointsOnHyperPara(const FSize N , const FReal &a, const FReal &b, const FReal &c, FReal * points) {
    //
    FReal u, v ;
    FSize j =0 ;
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
        u = 2.0*getRandom<FReal>() - 1.0 ;  v = 2.0*getRandom<FReal>() - 1.0 ;
        points[j]    = a*u ;
        points[j+1]  = b*v ;
        points[j+2]  = c*(u*u - v*v)  ;
    }
};


199
200
201
202
203
204
205
206
207
//!  \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....
//!
208
template <class FReal>
209
void unifRandonPointsOnSphere(const FSize N , const FReal R, FReal * points) {
210
211
	//
	unifRandonPointsOnUnitSphere(N , points) ;
212
213
    FSize j =0 ;
    for (FSize i = 0 ; i< N ; ++i, j+=4)  {
214
215
216
217
218
		points[j]	   *= R ;
		points[j+1]  *= R ;
		points[j+2]  *= R ;
	}
};
COULAUD Olivier's avatar
COULAUD Olivier committed
219
//!  \fn void plummerDist(int & cpt, const FReal &R)
220
221
222
223
224
225
226
227

//! \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
//!
228
template <class FReal>
229
FReal  plummerDist(FSize cpt, const FReal &R) {
230
231
232
	//
	FReal radius ,u ;
	do  {
COULAUD Olivier's avatar
COULAUD Olivier committed
233
		//
234
        u        = FMath::pow (getRandom<FReal>() , 2.0/3.0) ;
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
		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....
253
template <class FReal>
254
void unifRandonPlummer(const FSize N , const FReal R, const FReal M, FReal * points) {
255
256
257
258
	//
	unifRandonPointsOnUnitSphere(N , points) ;
	//
	FReal r , rm= 0.0;
259
    //	FReal Coeff =  3.0*M/(4.0*FMath::FPi<FReal>()*R*R*R) ;
260
	//am1 = 0 ;//1/FMath::pow(1+R*R,2.5);
261
262
    FSize cpt = 0 ;
    for (FSize i = 0,j=0 ; i< N ; ++i, j+=4)  {
263
264
265
266
267
268
269
270
271
		// 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: "
272
            <<100*static_cast<FReal>(cpt-N)/static_cast<FReal>(cpt) << " %" <<std::endl;
273
274

} ;
COULAUD Olivier's avatar
COULAUD Olivier committed
275
//
276
#endif