PeTar
N-body code for collisional gravitational systems
particle_distribution_generator.hpp
Go to the documentation of this file.
1 #pragma once
2 #include "Common/binary_tree.h"
3 
5 class ParticleDistributionGenerator {
6  // f(x) = x/sigma^2*exp(-x^2/(2sigma^2))
7  // F(x) = 1.0 - exp(-x^2/(2sigma^2))
8  // x = sqrt(-2*sigma^2*ln(1-F))
9  // 2sigma^2 = <e^2> or <i^2>
10  // <> means R.M.S.
11  static double RayleighDistribution(const double sigma){
12  static PS::MTTS mt;
13  static bool first = true;
14  if(first){
15  mt.init_genrand( PS::Comm::getRank() );
16  first = false;
17  }
18  double F = mt.genrand_res53();
19  /*
20  double ret = 0.0;
21  do{
22  ret = sqrt( -2.0*sigma*sigma*log(1.0-F));
23  }while(ret >= 1.0);
24  return ret;
25  */
26  return sqrt( -2.0*sigma*sigma*log(1.0-F));
27  }
28 
29  // p(a)da = C*a^-0.5*da
30  // P(a) = 2*C*(a^0.5 - a_in^0.5)
31  // a = (P/(2C) + a_in^0.5)^2
32  static double HayashiDistribution(const double a_in, const double a_out){
33  static PS::MTTS mt;
34  static bool first = true;
35  if(first){
36  mt.init_genrand( PS::Comm::getRank() );
37  first = false;
38  }
39  const double C = 0.5 / (sqrt(a_out) - sqrt(a_in));
40  double P = mt.genrand_res53();
41  double ret = P/(2*C) + sqrt(a_in);
42  ret *= ret;
43  return ret;
44  }
45 
46 
47  // p(a)da = C*a^(p+1)*da
48  // P(a) = 2*C*(a^0.5 - a_in^0.5)
49  static double HayashiDistributionWithIceLine(const double a_in, const double a_out,
50  const double a_ice, const double f_ice=1.0, const double p_sigma=-1.5){
51  //std::cout<<"a_in="<<a_in<<std::endl;
52  //const PS::F64 p_sigma = -1.5;
53  const PS::F64 p_mass = 2.0 + p_sigma;
54  //std::cout<<"p_mass="<<p_mass<<std::endl;
55  static PS::MTTS mt;
56  //std::cout<<"check 0"<<std::endl;
57  static bool first = true;
58  //std::cout<<"check 1"<<std::endl;
59  //std::cout<<"PS::Comm::getRank()="<<PS::Comm::getRank()<<std::endl;
60  if(first){
61  mt.init_genrand( PS::Comm::getRank() );
62  first = false;
63  }
64  //std::cout<<"first="<<first<<std::endl;
65  double P = mt.genrand_res53();
66  double ret = 0.0;
67  double C = 0.0;
68  double P_ice = 0.0;
69  if(a_ice <= a_in || a_ice >= a_out){
70  //std::cout<<"check a"<<std::endl;
71  C = p_mass / (pow(a_out, p_mass) - pow(a_in, p_mass));
72  ret = pow(P*p_mass/C + pow(a_in, p_mass), 1.0/p_mass);
73  }
74  else{
75  //std::cout<<"check b"<<std::endl;
76  C = p_mass / ( f_ice*(pow(a_out, p_mass)-pow(a_ice, p_mass)) + pow(a_ice, p_mass) - pow(a_in, p_mass) );
77  P_ice = C/p_mass*(pow(a_ice, p_mass) - pow(a_in, p_mass));
78  //std::cout<<"C="<<C<<std::endl;
79  //std::cout<<"P_ice="<<P_ice<<std::endl;
80  if(P < P_ice){
81  ret = pow(P*p_mass/C + pow(a_in, p_mass), 1.0/p_mass);
82  }
83  else{
84  ret = pow( ((P*p_mass/C + pow(a_in, p_mass) - pow(a_ice, p_mass))/f_ice + pow(a_ice, p_mass)), 1.0/p_mass); // beyond the ice line
85  }
86  }
87  return ret;
88  }
89 
90 
91  static double AU2CM(const double x){
92  return x * 1.49597871e13;
93  }
94 
95  static double CM2AU(const double x){
96  return x / 1.49597871e13;
97  }
98 
99 public:
100  // [M] = 1 [g]
101  // [L] = 1 [cm]
102  // [T] = 3.871e3[sec]
103  // [V] = 2.5833118e-4 [cm/sec] = 2.5833118e-9 [km/sec]
104  // 1yr = 8.146732e3[T]
105  // or
106  // G = 1 [L]^3/([M][T]^2)
107  // [L] = 1[AU] = 1.495978707e8[km]
108  // [M] = 1[Msun] = 1.989e30[kg]
109  // [T] = 5.02198050479e6 [sec] = 0.15924595715 [yr]
110  // [V] = 29.7886203575 [km/sec]
111  static void makeKeplerDisk(PS::F64 & mass_planet_glb,
112  PS::F64 *& mass,
113  PS::F64vec *& pos,
114  PS::F64vec *& vel,
115  const long long int n_glb,
116  const long long int n_loc,
117  const double a_in, // [AU]
118  const double a_out, // [AU]
119  const double e_rms, // normalized
120  const double i_rms, // normalized
121  const double dens = 10.0, // [g/cm^2]
122  const double mass_sun = 1.0, //[m_sun]
123  const double a_ice = 0.0,
124  const double f_ice = 1.0,
125  const double power = -1.5,
126  const int seed = 0,
127  const double gravitational_constant = 1.0
128  ){
129  static const double mass_sun_gram = 1.989e33; //[g]
130  PS::MTTS mt;
131  //mt.init_genrand( PS::Comm::getRank() );
132  mt.init_genrand( PS::Comm::getRank()+seed*PS::Comm::getNumberOfProc() );
133  static const double PI = atan(1.0) * 4.0;
134  const double AU = AU2CM(1.0);
135  //mass_planet_glb = 4.0 * PI * dens * ( sqrt(a_out) - sqrt(a_in) ) * AU * AU / mass_sun_gram; // [Msun]
136  mass_planet_glb = 2.0 * PI * dens / (2.0+power) * ( pow(a_out, 2.0+power) - pow(a_in, 2.0+power) ) * AU * AU / mass_sun_gram; // [Msun]
137  const double m_planet = mass_planet_glb / n_glb;
138  mass = new double[n_loc];
139  pos = new PS::F64vec[n_loc];
140  vel = new PS::F64vec[n_loc];
141  const double h = pow(2.0*m_planet / (3.0*mass_sun), 1.0/3.0);
142  PS::F64 e_ave = 0.0;
143  const PS::F64 e_sigma = sqrt(0.5*e_rms*e_rms); // this is right procedure
144  const PS::F64 i_sigma = sqrt(0.5*i_rms*i_rms);
145  //const PS::F64 e_sigma = e_rms;
146  //const PS::F64 i_sigma = i_rms;
147  for(long long int i=0; i<n_loc; i++){
148  mass[i] = m_planet;
149  COMM::Binary bin;
150  bin.semi = HayashiDistributionWithIceLine(a_in, a_out, a_ice, f_ice, power);
151  bin.ecc = RayleighDistribution(e_sigma) * h;
152  bin.incline = RayleighDistribution(i_sigma) * h;
153  bin.rot_horizon = 2.0 * PI * mt.genrand_res53();
154  bin.rot_self = 2.0 * PI * mt.genrand_res53();
155  double l = 2.0 * PI * mt.genrand_res53(); // mean anomayl
156  bin.ecca = bin.calcEccAnomaly(l, bin.ecc); // eccentric anomayl
157  bin.m1 = mass_sun;
158  bin.m2 = m_planet;
159  ParticleBase sun,planet;
160  bin.calcParticles(sun, planet, gravitational_constant);
161  pos[i]=planet.pos;
162  vel[i]=planet.vel;
163  e_ave += bin.ecc*bin.ecc;
164  }
165  PS::F64 e_ave_glb = sqrt(PS::Comm::getSum(e_ave)/n_glb);
166  if(PS::Comm::getRank() == 0){
167  std::cerr<<"e_ave_glb="<<e_ave_glb<<std::endl;
168  std::cerr<<"e_ave_hill="<<e_ave_glb/h<<std::endl;
169  std::cerr<<"e_rms="<<e_rms<<std::endl;
170  }
171  }
172 
173  static void makePlummerModel(const double mass_glb,
174  const long long int n_glb,
175  const long long int n_loc,
176  double *& mass,
177  PS::F64vec *& pos,
178  PS::F64vec *& vel,
179  const double eng = -0.25,
180  const int seed = 0){
181 
182  assert(eng < 0.0);
183  static const double PI = atan(1.0) * 4.0;
184  const double r_cutoff = 22.8 / (-3.0 * PI * mass_glb * mass_glb / (64.0 * -0.25)); // 22.8 is cutoff in Nbody units
185  //const double r_cutoff = 22.8 * 0.25;
186  mass = new double[n_loc];
187  pos = new PS::F64vec[n_loc];
188  vel = new PS::F64vec[n_loc];
189 
190  PS::MTTS mt;
191  mt.init_genrand( PS::Comm::getRank() );
192  for(int i=0; i<n_loc; i++){
193  mass[i] = mass_glb / n_glb;
194  double r_tmp = 9999.9;
195  while(r_tmp > r_cutoff){
196  double m_tmp = mt.genrand_res53();
197  r_tmp = 1.0 / sqrt( pow(m_tmp, (-2.0/3.0)) - 1.0);
198  }
199  double phi = 2.0 * PI * mt.genrand_res53();
200  double cth = 2.0 * (mt.genrand_real2() - 0.5);
201  double sth = sqrt(1.0 - cth*cth);
202  pos[i][0] = r_tmp * sth * cos(phi);
203  pos[i][1] = r_tmp * sth * sin(phi);
204  pos[i][2] = r_tmp * cth;
205  while(1){
206  const double v_max = 0.1;
207  const double v_try = mt.genrand_res53();
208  const double v_crit = v_max * mt.genrand_res53();
209  if(v_crit < v_try * v_try * pow( (1.0 - v_try * v_try), 3.5) ){
210  const double ve = sqrt(2.0) * pow( (r_tmp*r_tmp + 1.0), -0.25 );
211  phi = 2.0 * PI * mt.genrand_res53();
212  cth = 2.0 * (mt.genrand_res53() - 0.5);
213  sth = sqrt(1.0 - cth*cth);
214  vel[i][0] = ve * v_try * sth * cos(phi);
215  vel[i][1] = ve * v_try * sth * sin(phi);
216  vel[i][2] = ve * v_try * cth;
217  break;
218  }
219  }
220  }
221 
222  PS::F64vec cm_pos = 0.0;
223  PS::F64vec cm_vel = 0.0;
224  double cm_mass = 0.0;
225  for(PS::S32 i=0; i<n_loc; i++){
226  cm_pos += mass[i] * pos[i];
227  cm_vel += mass[i] * vel[i];
228  cm_mass += mass[i];
229  }
230  cm_pos /= cm_mass;
231  cm_vel /= cm_mass;
232  for(PS::S32 i=0; i<n_loc; i++){
233  pos[i] -= cm_pos;
234  vel[i] -= cm_vel;
235  }
236 
237  const double r_scale = -3.0 * PI * mass_glb * mass_glb / (64.0 * eng);
238  const double coef = 1.0 / sqrt(r_scale);
239  for(PS::S32 i=0; i<n_loc; i++){
240  pos[i] *= r_scale;
241  vel[i] *= coef;
242  }
243 
244  double r_max_sq = -1.0;
245  for(int i=0; i<n_loc; i++){
246  if(r_max_sq < pos[i] * pos[i]){
247  r_max_sq = pos[i] * pos[i];
248  }
249  }
250  }
251 
252 };
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
ParticleDistributionGenerator
Particle distribution generator.
Definition: particle_distribution_generator.hpp:5
ParticleDistributionGenerator::makePlummerModel
static void makePlummerModel(const double mass_glb, const long long int n_glb, const long long int n_loc, double *&mass, PS::F64vec *&pos, PS::F64vec *&vel, const double eng=-0.25, const int seed=0)
Definition: particle_distribution_generator.hpp:173
ParticleBase
Basic particle class.
Definition: particle_base.hpp:20
PIKG::F64vec
Vector3< F64 > F64vec
Definition: pikg_vector.hpp:167
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
ParticleBase::pos
PS::F64vec pos
Definition: particle_base.hpp:24
ParticleBase::vel
PS::F64vec vel
Definition: particle_base.hpp:25
ParticleDistributionGenerator::makeKeplerDisk
static void makeKeplerDisk(PS::F64 &mass_planet_glb, PS::F64 *&mass, PS::F64vec *&pos, PS::F64vec *&vel, const long long int n_glb, const long long int n_loc, const double a_in, const double a_out, const double e_rms, const double i_rms, const double dens=10.0, const double mass_sun=1.0, const double a_ice=0.0, const double f_ice=1.0, const double power=-1.5, const int seed=0, const double gravitational_constant=1.0)
Definition: particle_distribution_generator.hpp:111