PeTar
N-body code for collisional gravitational systems
galpy_interface.h
Go to the documentation of this file.
1 #include <iostream>
2 #include <sstream>
3 #include <iomanip>
4 #include <cassert>
5 #include <string>
6 #include <vector>
7 #include <getopt.h>
8 #include <cmath>
9 #include <algorithm>
10 #include "../src/io.hpp"
11 
12 #include <integrateFullOrbit.h>
13 
15 class IOParamsGalpy{
16 public:
17  IOParamsContainer input_par_store;
18  IOParams<std::string> type_args; // potential type and argument list
19  IOParams<std::string> pre_define_type; // pre defined type name
20  IOParams<std::string> config_filename; // configure file name
21  //IOParams<std::string> unit_set; // unscale, bovy
22  IOParams<double> rscale;
23  //IOParams<double> tscale;
24  IOParams<double> vscale;
25  //IOParams<double> fscale;
26  //IOParams<double> pscale;
27 
28  bool print_flag;
29 
30  IOParamsGalpy(): input_par_store(),
31  type_args (input_par_store, "__NONE__", "galpy-type-arg", "Add potential types and argumentsto the potential list in the center of the galactic reference frame"),
32  pre_define_type (input_par_store, "__NONE__", "galpy-set", "Add a pre-defined potential set to the potential list, options are: MWPotential2014"),
33  config_filename(input_par_store, "__NONE__", "galpy-conf-file", "A configure file of the time- and space-dependent potentials"),
34  //unit_set(input_par_store, "unscale", "galpy-units", "Units conversion set: 'unscale': no conversion; all scale factors are 1.0; 'bovy': radial unit: 8 kpc; velocity unit: 220 km/s"),
35  rscale(input_par_store, 1.0, "galpy-rscale", "Radius scale factor from unit of the input particle data (IN) to Galpy distance unit (1.0)"),
36  //tscale(input_par_store, 1.0, "galpy-tscale", "Time scale factor (rscale/vscale) from unit of the input particle data (IN) to Galpy time (1.0)"),
37  vscale(input_par_store, 1.0, "galpy-vscale", "Velocity scale factor from unit of the input particle data (IN) to Galpy velocity unit (1.0)"),
38  //fscale(input_par_store, 1.0, "galpy-fscale", "Acceleration scale factor (vscale^2/rscale) from unit of the input particle data (IN) to Galpy acceleration unit (1.0)"),
39  //pscale(input_par_store, 1.0, "galpy-pscale", "Potential scale factor (vscale^2) from unit of the input particle data (IN) to Galpy potential unit (1.0)"),
40  print_flag(false) {}
41 
43 
49  int read(int argc, char *argv[], const int opt_used_pre=0) {
50  static int galpy_flag=-1;
51  const struct option long_options[] = {
52  {type_args.key, required_argument, &galpy_flag, 0},
53  {config_filename.key, required_argument, &galpy_flag, 1},
54  {pre_define_type.key, required_argument, &galpy_flag, 2},
55  {rscale.key, required_argument, &galpy_flag, 3},
56  //{tscale.key, required_argument, &galpy_flag, 4},
57  {vscale.key, required_argument, &galpy_flag, 5},
58  //{fscale.key, required_argument, &galpy_flag, 6},
59  //{pscale.key, required_argument, &galpy_flag, 7},
60  {"help", no_argument, 0, 'h'},
61  {0,0,0,0}
62  };
63 
64  int opt_used=opt_used_pre;
65  int copt;
66  int option_index;
67  std::string fname_par;
68  optind = 0;
69  while ((copt = getopt_long(argc, argv, "-p:h", long_options, &option_index)) != -1)
70  switch (copt) {
71  case 0:
72  switch (galpy_flag) {
73  case 0:
74  type_args.value = optarg;
75  if (print_flag) type_args.print(std::cout);
76  opt_used+=2;
77  break;
78  case 1:
79  config_filename.value = optarg;
80  if(print_flag) config_filename.print(std::cout);
81  opt_used+=2;
82  break;
83  case 2:
84  pre_define_type.value = optarg;
85  if (print_flag) pre_define_type.print(std::cout);
86  opt_used+=2;
87  break;
88  case 3:
89  rscale.value = atof(optarg);
90  if(print_flag) rscale.print(std::cout);
91  opt_used+=2;
92  break;
93  //case 4:
94  // tscale.value = atof(optarg);
95  // if(print_flag) tscale.print(std::cout);
96  // opt_used+=2;
97  // break;
98  case 5:
99  vscale.value = atof(optarg);
100  if(print_flag) vscale.print(std::cout);
101  opt_used+=2;
102  break;
103  //case 6:
104  // fscale.value = atof(optarg);
105  // if(print_flag) fscale.print(std::cout);
106  // opt_used+=2;
107  // break;
108  //case 7:
109  // pscale.value = atof(optarg);
110  // if(print_flag) pscale.print(std::cout);
111  // opt_used+=2;
112  break;
113  default:
114  break;
115  }
116  break;
117  case 'p':
118  fname_par = optarg;
119  if(print_flag) {
120  std::string fgalpy_par = fname_par+".galpy";
121  FILE* fpar_in;
122  if( (fpar_in = fopen(fgalpy_par.c_str(),"r")) == NULL) {
123  fprintf(stderr,"Error: Cannot open file %s.\n", fgalpy_par.c_str());
124  abort();
125  }
126  input_par_store.readAscii(fpar_in);
127  fclose(fpar_in);
128  }
129  opt_used+=2;
130 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
131  input_par_store.mpi_broadcast();
132  PS::Comm::barrier();
133 #endif
134  break;
135  case 'h':
136  if(print_flag){
137  std::cout<<"Galpy options:"<<std::endl;
138  input_par_store.printHelp(std::cout, 2, 10, 23);
139  std::cout<<"***PS: use petar.galpy.help to check how to setup --galpy-type-arg and --galpy-conf-file."
140  <<std::endl;
141  }
142  return -1;
143  case '?':
144  opt_used +=2;
145  break;
146  default:
147  break;
148  }
149 
150  if(print_flag) std::cout<<"----- Finish reading input options of Galpy -----\n";
151 
152  return opt_used;
153  }
154 
156  void setBovyUnit(const bool print_flag=true) {
157  double kms_pcmyr=1.022712165045695; // pc = 30856775814913673 m, yr = 365.25 days
158  double vbase=220.0;
159  double rbase=8.0;
160  double vb_pcmyr = vbase*kms_pcmyr;
161  rscale.value = 0.001/rbase; // pc to solar position in Milkyway
162  vscale.value = 1.0/vb_pcmyr; // pc/Myr to solar velocity in Milkyway
163 
164  if (print_flag) {
165  double tscale = rscale.value/vscale.value; // myr to time unit in galpy
166  double fscale = vscale.value*vscale.value/rscale.value; // pc/Myr^2 to acc unit in galpy
167  double pscale = vscale.value*vscale.value; // pc^2/Myr^2 to pot unit in galpy
168  double GMscale = pscale*rscale.value; // pc^3/Myr^2
169  std::cout<<"----- Unit conversion for Galpy -----\n"
170  <<"rscale = "<<rscale.value<<" [8 kpc] / [pc]\n"
171  <<"vscale = "<<vscale.value<<" [220 km/s] / [pc/Myr]\n"
172  <<"tscale = "<<tscale <<" [Solar orbital period/2 pi] / Myr\n"
173  <<"fscale = "<<fscale <<" [galpy acceleration unit] / [pc/Myr^2]\n"
174  <<"pscale = "<<pscale <<" [galpy potential unit] / [pc^2/Myr^2]\n"
175  <<"GMscale = "<<GMscale<<" [galpy gravitational constant * mass unit] / [pc^3/Myr^2]\n";
176  }
177  }
178 };
179 
181 class FRWModel{
182 public:
183  double time; // current time (Myr)
184  double dt_max; // maximum time step to integrate a (Myr)
185  double a; //scale factor
186  double H0; // Hubble constant (km s^-1 Mpc^-1)
187  double H0_myr; // scaled Hubble constant (Myr^-1)
188  double omega_energy; // normalized dark energy density
189  double omega_radiation; // normalized radiation density
190  double omega_matter; // normalized matter density
191  double dt_scale; // if -1, evolve a backwards with -dt, if 1, evolve forwards with dt, but time integration is not multiplied by dt_scale
192 
193  inline double getMatterDensity(const double _a) const {
194  return omega_matter/_a;
195  }
196 
197  inline double getRadiationDensity(const double _a) const {
198  return omega_radiation/(_a*_a);
199  }
200 
201  inline double getEnergyDensity(const double _a) const {
202  return omega_energy*(_a*_a);
203  }
204 
205  inline double getCurvatureDensity() const {
206  return 1 - omega_energy - omega_radiation - omega_matter;
207  }
208 
209  inline double getDensity(const double _a) const {
210  return getMatterDensity(_a) + getRadiationDensity(_a) + getEnergyDensity(_a) + getCurvatureDensity();
211  }
212 
213  inline void calcH0Myr() {
214  double kms_pcmyr=1.022712165045695; // pc = 30856775814913673 m, yr = 365.25 days
215  H0_myr = H0*kms_pcmyr*1e-6;
216  }
217 
219 
223  inline double calcAdot(const double _a) const {
224  return H0_myr*sqrt(getDensity(_a));
225  }
226 
228 
231  void updateA(const double _dt) {
232  double adot = calcAdot(a);
233  double a_tmp = a + adot*dt_scale*_dt;
234 
235  for (int k=0; k<10; k++) {
236  double adot_new = calcAdot(a_tmp);
237  a_tmp = a + 0.5*(adot_new+adot)*dt_scale*_dt;
238  }
239  a = a_tmp;
240  time += _dt;
241  }
242 
244 
247  void updateAwithDtmin(const double _dt) {
248  assert(_dt>=0);
249  assert(dt_max>0);
250 
251  if (_dt==0) return;
252 
253  double dt = _dt;
254  int nstep = 1;
255  if (dt>dt_max) {
256  dt = dt_max;
257  nstep = int(_dt/dt_max);
258  }
259  for (int i=0; i<nstep; i++) updateA(dt);
260  if (dt*nstep<_dt) {
261  dt = _dt - dt*nstep;
262  updateA(dt);
263  }
264  }
265 
267  double getH() const {
268  return calcAdot(a)/a;
269  }
270 
271 };
272 
274 
277 class MWPotentialEvolve{
278 private:
279  template <class Tstream>
280  void labelCheck(Tstream& fconf, const char* match) {
281  std::string label;
282  fconf>>label;
283  if (label!=match) {
284  std::cerr<<"Galpy MWPotentialEvolve config: reading label error, should be "<<match<<" given "<<label<<std::endl;
285  abort();
286  }
287  }
288 
289  void eofCheck(std::ifstream& fconf, const char* message) {
290  if(fconf.eof()) {
291  std::cerr<<"Galpy MWPotentialEvolve config: reading "<<message<<" fails! File reaches EOF."<<std::endl;
292  abort();
293  }
294  }
295 
296 public:
297  FRWModel frw;
298  struct {
299  double m_vir_halo;
300  double c_halo; // concentration
301  double ac_halo; // formation epoch of halo (Wechsler 2002, Zhao 2003)
302  double m_disk;
303  double ra_disk;
304  double rb_disk;
305  double rho_bulge;
306  double alpha_bulge;
307  double rcut_bulge;
308  } init;
309  struct {
310  double m_vir_halo;
311  double r_vir_halo;
312  double rho_c_halo;
313  double gm_halo;
314  double rs_halo;
315  double gm_disk;
316  double ra_disk;
317  double rb_disk;
318  double grho_bulge;
319  double alpha_bulge;
320  double rcut_bulge;
321  } now;
322 
324 
343  void initialFromFile(const std::string& _filename, const double& _time) {
344 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
345  int my_rank = PS::Comm::getRank();
346  if (my_rank==0) {
347 #endif
348  std::ifstream fconf;
349  fconf.open(_filename.c_str(), std::ifstream::in);
350  if (!fconf.is_open()) {
351  std::cerr<<"Galpy MWPotentialEvolve config: configure file "<<_filename.c_str()<<" cannot be open!"<<std::endl;
352  abort();
353  }
354 
355  labelCheck(fconf, "Time");
356  fconf>>frw.time;
357  eofCheck(fconf, "current time (Myr)");
358  labelCheck(fconf, "a");
359  fconf>>frw.a;
360  eofCheck(fconf, "scale factor");
361  labelCheck(fconf, "H0");
362  fconf>>frw.H0;
363  eofCheck(fconf, "Hubble constant");
364  labelCheck(fconf, "Omega");
365  fconf>>frw.omega_energy>>frw.omega_radiation>>frw.omega_matter;
366  eofCheck(fconf, "cosmological densities");
367  labelCheck(fconf, "dt_scale");
368  fconf>>frw.dt_scale;
369  eofCheck(fconf, "time scale");
370  labelCheck(fconf, "dt_max");
371  fconf>>frw.dt_max;
372 
373  eofCheck(fconf, "maximum time step (Myr)");
374  labelCheck(fconf, "Halo");
375  fconf>>init.m_vir_halo>>init.c_halo>>init.ac_halo;
376  eofCheck(fconf, "Halo parameters");
377  labelCheck(fconf, "Disk");
378  fconf>>init.m_disk>>init.ra_disk>>init.rb_disk;
379  eofCheck(fconf, "Disk parameters");
380  labelCheck(fconf, "Bulge");
381  fconf>>init.rho_bulge>>init.alpha_bulge>>init.rcut_bulge;
382  eofCheck(fconf, "Bulge parameters");
383 
384  fconf.close();
385 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
386  }
387  PS::Comm::broadcast(&frw, 1, 0);
388  PS::Comm::broadcast(&init, 1, 0);
389 #endif
390 
391  if (frw.time!=_time) {
392  std::cerr<<"Error: the time recorded in the MWpotential parameter file, "<<frw.time<<", is not equal to the current system time, "<<_time<<"!"<<std::endl;
393  abort();
394  }
395 
396  frw.calcH0Myr();
397  }
398 
400 
408  void readDataFromFile(const std::string& _filename, const double& _time) {
409  const int npar=17;
410  double pars[npar];
411  #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
412  int my_rank = PS::Comm::getRank();
413  if (my_rank==0) {
414  #endif
415  std::ifstream fconf;
416  fconf.open(_filename.c_str(), std::ifstream::in);
417  if (!fconf.is_open()) {
418  std::cerr<<"Galpy MWPotentialEvolve config: configure file "<<_filename.c_str()<<" cannot be open!"<<std::endl;
419  abort();
420  }
421 
422  for (int i=0; i<npar; i++) {
423  fconf>>pars[i];
424  if (fconf.eof()) {
425  std::cerr<<"Galpy MWPotentialEvolve config: reading number of parameters is "<<npar<<" parameters, only "<<i<<" given!"<<std::endl;
426  abort();
427  }
428  }
429  fconf.close();
430 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
431  }
432  PS::Comm::broadcast(pars, npar, 0);
433 #endif
434 
435  frw.time = pars[0];
436  if (frw.time!=_time) {
437  std::cerr<<"Error: the time recorded in the MWpotential parameter file, "<<frw.time<<", is not equal to the current system time, "<<_time<<"!"<<std::endl;
438  abort();
439  }
440 
441  frw.a = pars[1];
442  frw.H0 = pars[2];
443 
444  frw.calcH0Myr();
445  frw.omega_energy = pars[3];
446  frw.omega_radiation = pars[4];
447  frw.omega_matter = pars[5];
448  frw.dt_scale = pars[6];
449  frw.dt_max = pars[7];
450 
451  init.m_vir_halo = pars[8];
452  init.c_halo = pars[9];
453  init.ac_halo = pars[10];
454  init.m_disk = pars[11];
455  init.ra_disk = pars[12];
456  init.rb_disk = pars[13];
457  init.rho_bulge = pars[14];
458  init.alpha_bulge = pars[15];
459  init.rcut_bulge = pars[16];
460  }
461 
463 
466  void writeDataToFile(const std::string& _filename) {
467  std::ofstream fout;
468  fout.open(_filename.c_str(), std::ifstream::out);
469  if (!fout.is_open()) {
470  std::cerr<<"Error: Galpy potential parameter file to write, "<<_filename<<", cannot be open!"<<std::endl;
471  abort();
472  }
473  fout<<std::setprecision(14)
474  <<frw.time<<" "
475  <<frw.a<<" "
476  <<frw.H0<<" "
477  <<frw.omega_energy<<" "
478  <<frw.omega_radiation<<" "
479  <<frw.omega_matter<<" "
480  <<frw.dt_scale<<" "
481  <<frw.dt_max<<" "
482  <<init.m_vir_halo<<" "
483  <<init.c_halo<<" "
484  <<init.ac_halo<<" "
485  <<init.m_disk<<" "
486  <<init.ra_disk<<" "
487  <<init.rb_disk<<" "
488  <<init.rho_bulge<<" "
489  <<init.alpha_bulge<<" "
490  <<init.rcut_bulge
491  <<std::endl;
492  fout.close();
493  }
494 
496 
498  void calcMWPotentialParameters(const double& _time) {
499 
500  double dt = _time-frw.time;
501  frw.updateAwithDtmin(dt);
502 
503  const double G_astro = 0.0044983099795944; // Msun, pc, Myr
504  const double pi = 3.141592653589793;
505 
506  double z = 1/frw.a-1; // redshift
507  double Ht = frw.getH(); // Myr^-1
508  double H0 = frw.H0_myr; // Myr^-1
509  double mass_density_m1 = frw.getMatterDensity(frw.a)/frw.getDensity(frw.a) - 1;
510  double mass_density_m1_z0 = frw.getMatterDensity(1.0) - 1;
511  double delta = 18*pi*pi + 82*mass_density_m1 - 39*mass_density_m1*mass_density_m1;
512  double delta0 = 18*pi*pi + 82*mass_density_m1_z0 - 39*mass_density_m1_z0*mass_density_m1_z0;
513  now.rho_c_halo = 3*Ht*Ht/(8*pi*G_astro); // astro unit
514  double rho_c_halo0 = 3*H0*H0/(8*pi*G_astro); // astro unit
515 
516  // halo, update Mvir, Rvir, c (Wechsler 2002, Zhao 2003) in astro unit
517  double mass_z_factor = std::exp(-2*init.ac_halo*z);
518  now.m_vir_halo = init.m_vir_halo*mass_z_factor; // Msun
519 
520  now.r_vir_halo = std::pow(3*now.m_vir_halo/(4*pi*delta*now.rho_c_halo),1.0/3.0);
521  double rvir_z_factor = std::pow(mass_z_factor*delta0*rho_c_halo0/(delta*now.rho_c_halo), 1.0/3.0);
522 
523  // Galpy parameter for halo in Galpy unit
524  double c0 = init.c_halo; // present-day concentration r_vir/r_s
525  double c = c0/(z+1); // concentration at reshift z
526  double c_factor = 1.0/(std::log(1+c) - c/(1+c));
527 
528  now.gm_halo = G_astro*now.m_vir_halo*c_factor;
529  now.rs_halo = now.r_vir_halo/c;
530 
531  // disk & budge, update mass and radial scale factors (Bullock & Johnston 2005)
532  now.gm_disk = G_astro*init.m_disk*mass_z_factor;
533  now.ra_disk = init.ra_disk*rvir_z_factor;
534  now.rb_disk = init.rb_disk*rvir_z_factor;
535 
536  now.grho_bulge = G_astro*init.rho_bulge*mass_z_factor/std::pow(rvir_z_factor,3-init.alpha_bulge);
537  now.alpha_bulge = init.alpha_bulge;
538  now.rcut_bulge = init.rcut_bulge*rvir_z_factor;
539 
540  }
541 };
542 
544 struct PotentialSetPar{
545  int mode; // mode of origin: 0: galactic frame; 1: local particle system frame
546  double gm; // G*mass of potential (used for computing acc from particle system)
547  double pos[3]; // position
548  double vel[3]; // velocity
549  double acc[3]; // acceleration
550 
551  PotentialSetPar(): mode(-1), gm(0), pos{0.0,0.0,0.0}, vel{0.0,0.0,0.0}, acc{0.0,0.0,0.0} {}
552 
554  void writeData(std::ostream& fout) {
555  fout<<mode<<" "
556  <<gm<<" ";
557  for (std::size_t i=0; i<3; i++) fout<<pos[i]<<" ";
558  for (std::size_t i=0; i<3; i++) fout<<vel[i]<<" ";
559  for (std::size_t i=0; i<3; i++) fout<<acc[i]<<" ";
560  }
561 
562  void readData(std::ifstream& fin) {
563  fin>>mode>>gm;
564  fin>>pos[0]>>pos[1]>>pos[2];
565  fin>>vel[0]>>vel[1]>>vel[2];
566  fin>>acc[0]>>acc[1]>>acc[2];
567  }
568 
570  void setOrigin(const int _mode, const double _gm=0, const double* _pos=NULL, const double* _vel=NULL, const double* _acc=NULL) {
571  assert(_mode>=0&&_mode<=2);
572  mode = _mode;
573  gm = _gm;
574  if (_pos!=NULL) {
575  pos[0] = _pos[0];
576  pos[1] = _pos[1];
577  pos[2] = _pos[2];
578  }
579  if (_vel!=NULL) {
580  vel[0] = _vel[0];
581  vel[1] = _vel[1];
582  vel[2] = _vel[2];
583  }
584  if (_acc!=NULL) {
585  acc[0] = _acc[0];
586  acc[1] = _acc[1];
587  acc[2] = _acc[2];
588  }
589  }
590 
591  void clear() {
592  mode = -1;
593  gm = 0.0;
594  pos[0] = pos[1] = pos[2] = 0.0;
595  vel[0] = vel[1] = vel[2] = 0.0;
596  acc[0] = acc[1] = acc[2] = 0.0;
597  }
598 
599 };
600 
602 struct PotentialSet{
603  int npot; // number of potential models
604  struct potentialArg* arguments; //potential arguments array for Gaply
605 
606  PotentialSet(): npot(0), arguments(NULL) {}
607 
609 
614  void generatePotentialArgs(const int _npot, int* _type, double* _arg) {
615  assert(_npot>0);
616  npot = _npot;
617  arguments = new struct potentialArg[npot];
618 #ifdef GALPY_VERSION_1_7_1
619  parse_leapFuncArgs_Full(npot, arguments, &_type, &_arg);
620 #else
621  parse_leapFuncArgs_Full(npot, arguments, &_type, &_arg, NULL);
622 #endif
623  }
624 
625  void clear() {
626  if (arguments!=NULL) {
627  assert(npot>0);
628  free_potentialArgs(npot, arguments);
629  free(arguments);
630  }
631  npot = 0;
632  arguments = NULL;
633  }
634 
635  ~PotentialSet() {
636  clear();
637  }
638 };
639 
641 struct ChangeArgument{
642  int index; // index of argument for change
643  int mode; // change mode (1: linear; 2: exponential)
644  double rate; // change rate: linear: cofficient a in a*t; exp: cofficient a in exp(a*t)
645 
647  void writeData(std::ostream& fout) {
648  fout<<index<<" "<<mode<<" "<<rate<<" ";
649  }
650 
652  void readData(std::ifstream& fin) {
653  fin>>index>>mode>>rate;
654  }
655 };
656 
658 class GalpyManager{
659 private:
660  template <class Tstream>
661  void labelCheck(Tstream& fconf, const char* match) {
662  std::string label;
663  fconf>>label;
664  if (label!=match) {
665  std::cerr<<"Galpy config: reading label error, should be "<<match<<" given "<<label<<std::endl;
666  abort();
667  }
668  }
669 
670  void eofCheck(std::ifstream& fconf, const char* message) {
671  if(fconf.eof()) {
672  std::cerr<<"Galpy config: reading "<<message<<" fails! File reaches EOF."<<std::endl;
673  abort();
674  }
675  }
676 
677  void gtZeroCheck(const int n, const char* message) {
678  if (n<=0) {
679  std::cerr<<"Galpy config: "<<message<<" <=0!"<<std::endl;
680  abort();
681  }
682  }
683 
684  void geZeroCheck(const int n, const char* message) {
685  if (n<0) {
686  std::cerr<<"Galpy config: "<<message<<" <0!"<<std::endl;
687  abort();
688  }
689  }
690 
691 
693 
699  template <class ttype>
700  void resizeArray(const int n_diff, const int index, std::vector<ttype>& array, std::vector<int>& array_offset) {
701  int offset = array_offset[index];
702  if (n_diff!=0) {
703  if (n_diff>0) {
704  std::vector<ttype> data(n_diff,ttype());
705  array.insert(array.begin()+offset, data.begin(), data.end());
706  }
707  else
708  array.erase(array.begin()+offset, array.begin()+offset-n_diff);
709  for (size_t i=index+1; i<array_offset.size(); i++)
710  array_offset[i] += n_diff;
711  }
712  }
713 
715 
720  template <class ttype>
721  void eraseArray(const int index, std::vector<ttype>& array, std::vector<int>& array_offset) {
722  int n = array_offset[index+1] - array_offset[index];
723  int offset = array_offset[index];
724  if (n>0) {
725  array.erase(array.begin()+offset, array.begin()+offset+n);
726  }
727  for (size_t i=index+1; i<array_offset.size(); i++)
728  array_offset[i] = array_offset[i+1]-n;
729  array_offset.pop_back();
730  assert(array_offset.back()==int(array.size()));
731  }
732 
733 
734 public:
735 
736  // for MPI communication, data IO and initialization
737  std::vector<int> pot_type_offset; // set offset
738  std::vector<int> pot_type; // types for each set
739  std::vector<int> pot_args_offset; // arguments of pot for each set
740  std::vector<double> pot_args; // set offset
741  double time; // current time, update in evolveChangingArguments
742  std::vector<ChangeArgument> change_args; // changing argument index to evolve
743  std::vector<int> change_args_offset; // set offset
744  // galpy potential arguments
745  std::vector<PotentialSetPar> pot_set_pars; // potential parameters for each set
746  std::vector<PotentialSet> pot_sets; // potential arguments for each set
747  double update_time;
748  // unit scaling
749  double rscale;
750  double vscale;
751  double tscale;
752  double fscale;
753  double pscale;
754  double gmscale;
755  // IO
756  std::ifstream fconf;
757  std::string set_name;
758  std::string set_parfile;
759  // evolving Milkyway potential
760  MWPotentialEvolve mw_evolve;
761 
762  GalpyManager(): pot_type_offset(), pot_type(),
763  pot_args_offset(), pot_args(), change_args(), change_args_offset(),
764  pot_set_pars(), pot_sets(), update_time(0.0), rscale(1.0), vscale(1.0), tscale(1.0), fscale(1.0), pscale(1.0), gmscale(1.0), fconf(), set_name(), set_parfile(), mw_evolve() {}
765 
767  void printData(std::ostream& fout) {
768  fout<<"Galpy parameters, time: "<<time;
769  fout<<" Next update time: "<<update_time;
770  fout<<std::endl;
771  int nset = pot_set_pars.size();
772  for (int k=0; k<nset; k++) {
773  auto& pot_set_par_k = pot_set_pars[k];
774  fout<<"Potential set "<<k+1<<" Mode: "<<pot_set_par_k.mode
775  <<" GM: "<<pot_set_par_k.gm
776  <<" Pos: "<<pot_set_par_k.pos[0]<<" "<<pot_set_par_k.pos[1]<<" "<<pot_set_par_k.pos[2]
777  <<" Vel: "<<pot_set_par_k.vel[0]<<" "<<pot_set_par_k.vel[1]<<" "<<pot_set_par_k.vel[2]
778  <<" Acc: "<<pot_set_par_k.acc[0]<<" "<<pot_set_par_k.acc[1]<<" "<<pot_set_par_k.acc[2]
779  <<"\nPotential type indice: ";
780  for (int i=pot_type_offset[k]; i<pot_type_offset[k+1]; i++)
781  fout<<pot_type[i]<<" ";
782  fout<<"\nPotential arguments: ";
783  for (int i=pot_args_offset[k]; i<pot_args_offset[k+1]; i++)
784  fout<<pot_args[i]<<" ";
785  fout<<"\nChange argument [index mode rate]:";
786  for (int i=change_args_offset[k]; i<change_args_offset[k+1]; i++) {
787  fout<<"["<<change_args[i].index
788  <<" "<<change_args[i].mode
789  <<" "<<change_args[i].rate<<"] ";
790  }
791  fout<<std::endl;
792 
793  if (set_name=="MWPotentialEvolve")
794  fout<<"Scale factor a: "<<mw_evolve.frw.a<<std::endl;
795  }
796  }
797 
798 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
799  void broadcastDataMPI() {
800 
801  int nset;
802  int my_rank = PS::Comm::getRank();
803  if (my_rank==0) nset = pot_set_pars.size();
804  PS::Comm::broadcast(&nset, 1, 0);
805 
806  // update potentials
807  if (nset>=0) {
808 
809  PS::Comm::broadcast(&update_time, 1, 0);
810  if (my_rank==0) {
811  assert((int)pot_type_offset.size()==nset+1);
812  assert((int)pot_args_offset.size()==nset+1);
813  assert((int)change_args_offset.size()==nset+1);
814  }
815  else {
816  pot_set_pars.resize(nset);
817  pot_type_offset.resize(nset+1);
818  pot_args_offset.resize(nset+1);
819  change_args_offset.resize(nset+1);
820  }
821  PS::Comm::broadcast(pot_set_pars.data(), pot_set_pars.size(), 0);
822  PS::Comm::broadcast(pot_type_offset.data(), pot_type_offset.size(), 0);
823  PS::Comm::broadcast(pot_args_offset.data(), pot_args_offset.size(), 0);
824  PS::Comm::broadcast(change_args_offset.data(), change_args_offset.size(), 0);
825 
826  if (my_rank==0) {
827  assert((int)pot_type.size()==pot_type_offset.back());
828  assert((int)pot_args.size()==pot_args_offset.back());
829  assert((int)change_args.size()==change_args_offset.back());
830  }
831  else {
832  pot_type.resize(pot_type_offset.back());
833  pot_args.resize(pot_args_offset.back());
834  change_args.resize(change_args_offset.back());
835  }
836  PS::Comm::broadcast(pot_type.data(), pot_type.size(), 0);
837  if (pot_args.size()>0) PS::Comm::broadcast(pot_args.data(), pot_args.size(), 0);
838  if (change_args.size()>0) PS::Comm::broadcast(change_args.data(), change_args.size(), 0);
839  }
840  }
841 
842 #endif
843 
845 
852  void initial(const IOParamsGalpy& _input, const double _time, const std::string _conf_name, const bool _restart_flag, const bool _print_flag=false) {
853  // unit scale
854  rscale = _input.rscale.value;
855  vscale = _input.vscale.value;
856  tscale = rscale/vscale;
857  fscale = vscale*vscale/rscale;
858  pscale = vscale*vscale;
859  gmscale = pscale*rscale;
860 
861  // initial offset
862  pot_type_offset.push_back(0);
863  pot_args_offset.push_back(0);
864  change_args_offset.push_back(0);
865 
866  set_name = _input.pre_define_type.value;
867  std::size_t ipar = set_name.find_first_of(":");
868  if (ipar!=std::string::npos) {
869  set_parfile = set_name.substr(ipar+1);
870  set_name = set_name.substr(0, ipar);
871  }
872 
873  // if restart, read the configure file from _conf_name instead of _input.pre_define_type.
874  if (_restart_flag) set_parfile = _conf_name;
875 
876  if (set_name=="MWPotential2014") {
877 
878  pot_set_pars.emplace_back();
879  auto& pot_set_par_k = pot_set_pars.back();
880  pot_set_par_k.setOrigin(0);
881 
882  pot_type.insert(pot_type.end(), {15,5,9});
883  pot_type_offset.push_back(pot_type.size());
884 
885  // Bovy unit
886  //double pot_args[8] = {0.0299946, 1.8, 0.2375,
887  // 0.7574802, 0.375, 0.035,
888  // 4.85223053, 2.0};
889 
890  // Astro unit
891  // PowerSphericalPotentialwithCutoff bulge
892  // rho(r) = (3-a) G*M_g / (4 pi R_g^(3-a)) * (r_1/r)^a * exp ((-r/rc)^2)
893  // args(3): (3-a) G*M_g / (4 pi R_g^(3-a)), a, rc (r_1 = 1)
894  // unit: GMscale/rscale^(3-a), 1, rscale
895  // MiyamotoNagai disk
896  // Phi(R,z) = -G*M_g / sqrt(R^2+(a+sqrt(z^2+b^2))^2)
897  // args(3): G*M_g, a, b
898  // unit: GMscale, rscale, rscale
899  // NFW halo
900  // rho(r) = G*M_g /(4 pi a^3) * 1/((r/a)*(1+r/a)^2)
901  // args(2): G*M_g, a
902  // unit: GMscale, rscale
903  pot_args.insert(pot_args.end(), {251.63858935563147, 1.8, 1899.9999999999998,
904  306770418.38588977, 3000.0, 280.0,
905  1965095308.192175, 16000.0});
906  pot_args_offset.push_back(pot_args.size());
907 
908  change_args_offset.push_back(0);
909 
910  }
911 
912  // MWpotentialEvolve should be the first set
913  if (set_name=="MWPotentialEvolve") {
914 
915  if (_restart_flag)
916  mw_evolve.readDataFromFile(set_parfile, _time);
917  else
918  mw_evolve.initialFromFile(set_parfile, _time);
919  updateMWPotentialEvolve(_time, _print_flag, true);
920  }
921 
922  // add pre-defined type-argu groups
923  std::string type_args = _input.type_args.value;
924  // Update types and arguments from type-args string
925  addPotentialFromString(type_args, true, _print_flag);
926 
927  // add type arguments from configure file if exist
928  if (_input.config_filename.value!="__NONE__") {
929  set_name="configure";
930 
931 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
932  int my_rank = PS::Comm::getRank();
933  if (my_rank==0) {
934 #endif
935  fconf.open(_input.config_filename.value.c_str(), std::ifstream::in);
936  if (!fconf.is_open()) {
937  std::cerr<<"Error: Galpy configure file "<<_input.config_filename.value.c_str()<<" cannot be open!"<<std::endl;
938  abort();
939  }
940  labelCheck(fconf, "Time");
941  if(fconf.eof()) fconf.close();
942  fconf>>update_time;
943 
944 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
945  }
946  PS::Comm::broadcast(&update_time, 1, 0);
947 #endif
948  // if restart, read configure file at the current time
949  if (_restart_flag) {
950  // first escape past time to avoid overwritting restart configure file
951  updatePotentialFromFile(_time, false);
952  readDataFromFile(set_parfile, _print_flag);
953  }
954  }
955 
956  updatePotentialSet();
957 
958  if(_print_flag) {
959  printData(std::cout);
960  std::cout<<"----- Finish initialize Galpy potential -----\n";
961  }
962 
963  }
964 
966 
969  void updatePotentialSet() {
970  freePotentialArgs();
971  int nset = pot_set_pars.size();
972  pot_sets.resize(nset);
973  for (int k=0; k<nset; k++) {
974  pot_sets[k].generatePotentialArgs(pot_type_offset[k+1]-pot_type_offset[k], &(pot_type[pot_type_offset[k]]), &(pot_args[pot_args_offset[k]]));
975  }
976  }
977 
979 
986  bool addPotentialFromString(const std::string& _type_args, const bool _reset_flag, const bool _print_flag) {
987  if (_type_args!="__NONE__") {
988 
989  pot_set_pars.emplace_back();
990  auto& pot_set_par_k = pot_set_pars.back();
991  pot_set_par_k.setOrigin(0);
992 
993  std::vector<std::string> type_args_pair;
994 
995  // split type-arg groups
996  std::size_t istart = 0;
997  std::size_t inext = _type_args.find_first_of("|");
998  while (inext!=std::string::npos) {
999  type_args_pair.push_back(_type_args.substr(istart,inext-istart));
1000  istart = inext+1;
1001  inext = _type_args.find_first_of("|",istart);
1002  }
1003  if (istart!=_type_args.size()) type_args_pair.push_back(_type_args.substr(istart));
1004 
1005  if (_print_flag) std::cout<<"Galpy Potential combination list:\n";
1006  // loop group
1007  for (std::size_t i=0; i<type_args_pair.size(); i++) {
1008  // get types
1009  istart = 0;
1010  inext = type_args_pair[i].find_first_of(":");
1011  if (inext==std::string::npos) {
1012  std::cerr<<"Error: potential type index delimiter ':' is not found in input string "<<type_args_pair[i]<<std::endl;
1013  abort();
1014  }
1015  std::string type_str = type_args_pair[i].substr(0,inext);
1016  std::string arg_str = type_args_pair[i].substr(inext+1);
1017 
1018  inext = type_str.find_first_of(",");
1019  if (_print_flag) std::cout<<"Type index: ";
1020  while (inext!=std::string::npos) {
1021  int type_i=std::stoi(type_str.substr(istart, inext-istart));
1022  pot_type.push_back(type_i);
1023  if (_print_flag) std::cout<<type_i<<" ";
1024  istart = inext+1;
1025  inext = type_str.find_first_of(",",istart);
1026  }
1027  if (istart!=type_str.size()) {
1028  int type_i=std::stoi(type_str.substr(istart));
1029  pot_type.push_back(type_i);
1030  if (_print_flag) std::cout<<type_i<<" ";
1031  }
1032  if (_print_flag) std::cout<<"args:";
1033 
1034  // get arguments
1035  istart = 0;
1036  inext = arg_str.find_first_of(",",istart);
1037  while (inext!=std::string::npos) {
1038  double arg_i = std::stod(arg_str.substr(istart, inext-istart));
1039  pot_args.push_back(arg_i);
1040  if (_print_flag) std::cout<<" "<<arg_i;
1041  istart = inext+1;
1042  inext = arg_str.find_first_of(",",istart);
1043  }
1044  if(istart!=arg_str.size()) {
1045  double arg_i = std::stod(arg_str.substr(istart));
1046  pot_args.push_back(arg_i);
1047  if (_print_flag) std::cout<<" "<<arg_i;
1048  }
1049  if (_print_flag) std::cout<<std::endl;
1050  }
1051  pot_type_offset.push_back(pot_type.size());
1052  pot_args_offset.push_back(pot_args.size());
1053  change_args_offset.push_back(0);
1054 
1055  return true;
1056  }
1057  return false;
1058  }
1059 
1061 
1067  bool updateMWPotentialEvolve(const double& _system_time, const bool _print_flag, const bool _initial_flag) {
1068 
1069  mw_evolve.calcMWPotentialParameters(_system_time);
1070  auto& mwpot = mw_evolve.now;
1071  double pot_arg_mw[8] = {mwpot.grho_bulge*gmscale, mwpot.alpha_bulge, mwpot.rcut_bulge*rscale,
1072  mwpot.gm_disk*gmscale, mwpot.ra_disk*rscale, mwpot.rb_disk*rscale,
1073  mwpot.gm_halo*gmscale, mwpot.rs_halo*rscale};
1074 
1075  if (_initial_flag) {
1076  pot_set_pars.emplace_back();
1077  auto& pot_set_par_k = pot_set_pars.back();
1078  pot_set_par_k.setOrigin(0);
1079 
1080  pot_type.insert(pot_type.end(), {15,5,9});
1081  pot_type_offset.push_back(pot_type.size());
1082 
1083  pot_args.insert(pot_args.end(), pot_arg_mw, pot_arg_mw+8);
1084  pot_args_offset.push_back(pot_args.size());
1085 
1086  change_args_offset.push_back(0);
1087  }
1088  else{
1089  assert(pot_type[0]==15);
1090  assert(pot_type[1]==5);
1091  assert(pot_type[2]==9);
1092  for (int i=0; i<8; i++) pot_args[i] = pot_arg_mw[i];
1093  }
1094 
1095  if (_print_flag) {
1096  std::cout<<std::setprecision(14)
1097  <<" Evolve MilkyWay potential: FRW a: "<<mw_evolve.frw.a
1098  <<" NFW Halo: M_vir[Msun]: "<<mwpot.m_vir_halo<<" R_vir[pc]: "<<mwpot.r_vir_halo<<std::endl;
1099  }
1100 
1101  return true;
1102  }
1103 
1105 
1111  bool updatePotentialFromFile(const double& _system_time, const bool _print_flag) {
1112  double time_scaled= _system_time*tscale;
1113  int nset=-1;
1114 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1115  int my_rank = PS::Comm::getRank();
1116  if (my_rank==0) {
1117 #endif
1118  if (fconf.is_open()) {
1119  if (time_scaled>=update_time) {
1120 
1121  double update_time_next;
1122 
1123  while(true) {
1124  labelCheck(fconf, "Task");
1125  std::string task;
1126  fconf>>task;
1127  eofCheck(fconf, "task"); int n_set_old = pot_set_pars.size();
1128  if (task=="add") {
1129  labelCheck(fconf, "Nset");
1130  fconf>>nset;
1131  eofCheck(fconf, "number of new potential sets");
1132  gtZeroCheck(nset, "number of new potential sets");
1133 
1134  for (int k=0; k<nset; k++) {
1135  pot_set_pars.push_back(PotentialSetPar());
1136  auto& pot_set_par_k = pot_set_pars.back();
1137  // set
1138  labelCheck(fconf, "Set");
1139  int set_index;
1140  fconf>>set_index;
1141  eofCheck(fconf, "set index");
1142  if (set_index!=k) {
1143  std::cerr<<"Galpy config: reading set index not match the expected one, should be "<<k<<" given "<<set_index<<std::endl;
1144  abort();
1145  }
1146  // Ntype
1147  labelCheck(fconf, "Ntype");
1148  int n_pot_k=-1;
1149  fconf>>n_pot_k;
1150  eofCheck(fconf, "number of potential types");
1151  gtZeroCheck(n_pot_k, "number of potential types");
1152 
1153  // mode
1154  labelCheck(fconf, "Mode");
1155  int mode_k;
1156  fconf>>mode_k;
1157  eofCheck(fconf, "mode");
1158  if (mode_k<0 || mode_k>2) {
1159  std::cerr<<"Galpy config: mode should be 0, 1, 2, given "<<mode_k<<std::endl;
1160  abort();
1161  }
1162  pot_set_par_k.mode = mode_k;
1163 
1164  // Mass
1165  labelCheck(fconf, "GM");
1166  fconf>>pot_set_par_k.gm;
1167  eofCheck(fconf, "g*m");
1168  // Pos
1169  labelCheck(fconf, "Pos");
1170  for (int j=0; j<3; j++) fconf>>pot_set_par_k.pos[j];
1171  eofCheck(fconf, "position");
1172  // Vel
1173  labelCheck(fconf, "Vel");
1174  for (int j=0; j<3; j++) fconf>>pot_set_par_k.vel[j];
1175  eofCheck(fconf, "velocity");
1176 
1177  // type, arg and change_args are not directly save to potential set because of MPI communication
1178  // types
1179  labelCheck(fconf, "Type");
1180  for (int i=0; i<n_pot_k; i++) {
1181  int type_i;
1182  fconf>>type_i;
1183  eofCheck(fconf, "potential type");
1184  pot_type.push_back(type_i);
1185  }
1186  pot_type_offset.push_back(pot_type.size());
1187 
1188  // arguments
1189  std::string line;
1190  std::getline(fconf, line); // skip 2nd line
1191  std::getline(fconf, line);
1192  std::istringstream fin(line);
1193  double arg_i;
1194  int n_arg=0;
1195  labelCheck(fin, "Arg");
1196  while (fin>>arg_i) {
1197  pot_args.push_back(arg_i);
1198  n_arg++;
1199  }
1200  pot_args_offset.push_back(pot_args.size());
1201 
1202  // time-dependent arguments
1203  labelCheck(fconf, "Nchange");
1204  int n_change;
1205  fconf>>n_change;
1206  eofCheck(fconf, "number of changing arguments");
1207  geZeroCheck(n_change, "number of changing arguments");
1208  int offset = change_args.size();
1209  if (n_change>0) {
1210  change_args.resize(offset+n_change);
1211  labelCheck(fconf, "Index");
1212  for (int i=0; i<n_change; i++) {
1213  int index_i;
1214  fconf>>index_i;
1215  eofCheck(fconf, "index of the changing argument");
1216  if (index_i<-1 || index_i>=n_arg) {
1217  std::cerr<<"Galpy config: index of the changing argument should be >=-1 and <"<<n_arg<<", given "<<index_i<<std::endl;
1218  abort();
1219  }
1220  change_args[offset+i].index = index_i;
1221  }
1222  labelCheck(fconf, "ChangeMode");
1223  for (int i=0; i<n_change; i++) {
1224  int change_mode_i;
1225  fconf>>change_mode_i;
1226  eofCheck(fconf, "changing mode");
1227  if (!(change_mode_i==1||change_mode_i==2)) {
1228  std::cerr<<"Galpy config: change mode must be 1 or 2, given "<<change_mode_i<<std::endl;
1229  abort();
1230  }
1231  change_args[offset+i].mode = change_mode_i;
1232  }
1233  labelCheck(fconf, "ChangeRate");
1234  for (int i=0; i<n_change; i++) {
1235  double rate_i;
1236  fconf>>rate_i;
1237  eofCheck(fconf, "changing rate");
1238  change_args[offset+i].rate = rate_i;
1239  }
1240  }
1241  change_args_offset.push_back(change_args.size());
1242  }
1243  nset += n_set_old;
1244  }
1245  else if (task=="update") {
1246  labelCheck(fconf, "Nset");
1247  fconf>>nset;
1248  eofCheck(fconf, "number of update potential sets");
1249  gtZeroCheck(nset, "number of update potential sets");
1250 
1251  for (int i=0; i<nset; i++) {
1252  // set
1253  labelCheck(fconf, "Set");
1254  int k;
1255  fconf>>k;
1256  eofCheck(fconf, "set index");
1257  if (k<0 || k>=n_set_old) {
1258  std::cerr<<"Galpy config: changing set index incorrect, should be from 0 to "<<n_set_old-1<<", given "<<k<<std::endl;
1259  abort();
1260  }
1261  auto& pot_set_par_k = pot_set_pars[k];
1262 
1263  // Ntype
1264  labelCheck(fconf, "Ntype");
1265  int n_pot_k=-1;
1266  fconf>>n_pot_k;
1267  eofCheck(fconf, "number of potential types");
1268  gtZeroCheck(n_pot_k, "number of potential types");
1269 
1270  // mode
1271  labelCheck(fconf, "Mode");
1272  int mode_k;
1273  fconf>>mode_k;
1274  eofCheck(fconf, "mode");
1275  if (mode_k<-2 || mode_k==-1 || mode_k>2) {
1276  std::cerr<<"Galpy config: mode should be 0, 1, 2 or -2, given "<<mode_k<<std::endl;
1277  abort();
1278  }
1279 
1280  // mass
1281  labelCheck(fconf, "GM");
1282  fconf>>pot_set_par_k.gm;
1283  eofCheck(fconf, "gm");
1284 
1285  // update position and velocity for mode 2
1286  if (mode_k!=-2) {
1287  labelCheck(fconf, "Pos");
1288  for (int j=0; j<3; j++) fconf>>pot_set_par_k.pos[j];
1289  eofCheck(fconf, "position");
1290  labelCheck(fconf, "Vel");
1291  for (int j=0; j<3; j++) fconf>>pot_set_par_k.vel[j];
1292  eofCheck(fconf, "velocity");
1293  }
1294  else {
1295  double dtmp;
1296  labelCheck(fconf, "Pos");
1297  for (int j=0; j<3; j++) fconf>>dtmp;
1298  eofCheck(fconf, "position");
1299  labelCheck(fconf, "Vel");
1300  for (int j=0; j<3; j++) fconf>>dtmp;
1301  eofCheck(fconf, "velocity");
1302  mode_k = 2;
1303  }
1304  pot_set_par_k.mode = mode_k;
1305 
1306  // types
1307  // modify array offset if number of pot in the set changes
1308  int n_pot_diff = n_pot_k - (pot_type_offset[k+1]-pot_type_offset[k]);
1309  resizeArray(n_pot_diff, k, pot_type, pot_type_offset);
1310  int offset = pot_type_offset[k];
1311 
1312  labelCheck(fconf, "Type");
1313  for (int j=0; j<n_pot_k; j++) {
1314  int type_j;
1315  fconf>>type_j;
1316  eofCheck(fconf, "potential type");
1317  pot_type[offset+j] = type_j;
1318  }
1319 
1320  // arguments
1321  std::string line;
1322  std::getline(fconf, line); // skip 2nd line
1323  std::getline(fconf, line);
1324  std::istringstream fin(line);
1325  labelCheck(fin, "Arg");
1326  std::vector<double> arg_k;
1327  double arg_tmp;
1328  while (fin>>arg_tmp) {
1329  arg_k.push_back(arg_tmp);
1330  }
1331  int n_arg_diff = int(arg_k.size()) - (pot_args_offset[k+1] - pot_args_offset[k]);
1332  resizeArray(n_arg_diff, k, pot_args, pot_args_offset);
1333  offset = pot_args_offset[k];
1334  for (size_t j=0; j<arg_k.size(); j++) {
1335  pot_args[offset+j] = arg_k[j];
1336  }
1337 
1338  // time-dependent arguments
1339  labelCheck(fconf, "Nchange");
1340  int n_change;
1341  fconf>>n_change;
1342  eofCheck(fconf, "number of changing arguments");
1343  geZeroCheck(n_change, "changing arguments");
1344  int n_change_diff = n_change - (change_args_offset[k+1] - change_args_offset[k]);
1345 
1346  resizeArray(n_change_diff, k, change_args, change_args_offset);
1347  if (n_change>0) {
1348  offset = change_args_offset[k];
1349  labelCheck(fconf, "Index");
1350  for (int j=0; j<n_change; j++) {
1351  int index_j;
1352  fconf>>index_j;
1353  eofCheck(fconf, "index of the changing argument");
1354  if (index_j<-1 || index_j>=int(arg_k.size())) {
1355  std::cerr<<"Galpy config: index of the changing argument should be >=-1 and <"<<(arg_k.size())<<", given "<<index_j<<std::endl;
1356  abort();
1357  }
1358  change_args[offset+j].index = index_j;
1359  }
1360  labelCheck(fconf, "ChangeMode");
1361  for (int j=0; j<n_change; j++) {
1362  int change_mode_j;
1363  fconf>>change_mode_j;
1364  eofCheck(fconf, "changing mode");
1365  if (!(change_mode_j==1||change_mode_j==2)) {
1366  std::cerr<<"Galpy config: change mode must be 1 or 2, given "<<change_mode_j<<std::endl;
1367  abort();
1368  }
1369  change_args[offset+j].mode = change_mode_j;
1370  }
1371  labelCheck(fconf, "ChangeRate");
1372  for (int j=0; j<n_change; j++) {
1373  double rate_j;
1374  fconf>>rate_j;
1375  eofCheck(fconf, "changing rate");
1376  change_args[offset+j].rate = rate_j;
1377  }
1378  }
1379  }
1380  nset = n_set_old;
1381  }
1382  else if (task=="remove") {
1383  labelCheck(fconf, "Nset");
1384  fconf>>nset;
1385  eofCheck(fconf, "number of removing potential sets");
1386  gtZeroCheck(nset, "number of removing potential sets");
1387 
1388  labelCheck(fconf, "Index");
1389  std::vector<int> rm_indices;
1390  for (int i=0; i<nset; i++) {
1391  int k;
1392  fconf>>k;
1393  eofCheck(fconf, "removing set index");
1394  rm_indices.push_back(k);
1395  }
1396  std::sort(rm_indices.begin(), rm_indices.end());
1397  for (int i=0; i<nset; i++) {
1398  int idel = rm_indices[i]-i; // each time of remove, the following indices reduce by one, total reducing number is i
1399  pot_set_pars.erase(pot_set_pars.begin()+idel);
1400  eraseArray(idel, pot_type, pot_type_offset);
1401  eraseArray(idel, pot_args, pot_args_offset);
1402  eraseArray(idel, change_args, change_args_offset);
1403  }
1404  nset = n_set_old - nset;
1405  geZeroCheck(nset, "number of potential sets after remove");
1406  }
1407  else {
1408  std::cerr<<"Galpy config: Task must be add, update or remove, given "<<task<<std::endl;
1409  abort();
1410  }
1411 
1412  std::string label;
1413  fconf>>label;
1414  if(fconf.eof()) {
1415  fconf.close();
1416  break;
1417  }
1418  if(label!="Time") {
1419  std::cerr<<"Galpy config: reading label error, should be Time given "<<label<<std::endl;
1420  abort();
1421  }
1422  fconf>>update_time_next;
1423  if (update_time_next>time_scaled) {
1424  update_time = update_time_next;
1425  break;
1426  }
1427  };
1428  }
1429  }
1430 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1431  }
1432  broadcastDataMPI();
1433  PS::Comm::broadcast(&nset, 1, 0);
1434 #endif
1435  return (nset>=0);
1436  }
1437 
1439 
1443  void updatePotential(const double& _system_time, const bool _print_flag) {
1444 
1445  int nset = pot_set_pars.size();
1446  assert((int)pot_type_offset.size()==nset+1);
1447  assert((int)pot_args_offset.size()==nset+1);
1448  assert((int)change_args_offset.size()==nset+1);
1449  assert((int)pot_type.size()==pot_type_offset.back());
1450  assert((int)pot_args.size()==pot_args_offset.back());
1451  assert((int)change_args.size()==change_args_offset.back());
1452 
1453  bool update_flag = false;
1454  if (set_name=="MWPotentialEvolve") {
1455  update_flag = updateMWPotentialEvolve(_system_time, false, false);
1456  }
1457  else if (set_name=="configure") {
1458  update_flag = updatePotentialFromFile(_system_time, _print_flag);
1459  // update time-dependent argument
1460  bool change_flag = evolveChangingArguments(_system_time*tscale, false);
1461  if (update_flag && _print_flag) printData(std::cout);
1462  update_flag = (update_flag || change_flag);
1463  }
1464 
1465  // generate galpy potential argument
1466  if (update_flag) updatePotentialSet();
1467  }
1468 
1470 
1474  void writePotentialPars(const std::string& _filename, const double& _system_time) {
1475  if (set_name=="MWPotentialEvolve") {
1476  assert(mw_evolve.frw.time==_system_time);
1477  mw_evolve.writeDataToFile(_filename);
1478  }
1479  else if (set_name=="configure") {
1480  assert(time==_system_time*tscale);
1481  writeDataToFile(_filename);
1482  }
1483  }
1484 
1486  static void printReference(std::ostream & fout, const int offset=4) {
1487  for (int i=0; i<offset; i++) fout<<" ";
1488  fout<<"Galpy: Bovy J., 2015, ApJS, 216, 29"<<std::endl;
1489  }
1490 
1492  void resetPotAcc() {
1493  for (size_t k=0; k<pot_set_pars.size(); k++) {
1494  double* acc_pot = pot_set_pars[k].acc;
1495  acc_pot[0] = acc_pot[1] = acc_pot[2] = 0.0;
1496  }
1497  }
1498 
1500 
1503  void kickMovePot(const double dt) {
1504  double dt_scale = dt*tscale;
1505  for (size_t k=0; k<pot_set_pars.size(); k++) {
1506  int mode_k = pot_set_pars[k].mode;
1507  assert(mode_k>=0||mode_k<=2);
1508  if (mode_k==2) {
1509  double* vel_k = pot_set_pars[k].vel;
1510  double* acc_k = pot_set_pars[k].acc;
1511  vel_k[0] += acc_k[0]*dt_scale;
1512  vel_k[1] += acc_k[1]*dt_scale;
1513  vel_k[2] += acc_k[2]*dt_scale;
1514  }
1515  }
1516  }
1517 
1519 
1522  void driftMovePot(const double dt) {
1523  double dt_scale = dt*tscale;
1524  for (size_t k=0; k<pot_set_pars.size(); k++) {
1525  int mode_k = pot_set_pars[k].mode;
1526  assert(mode_k>=0||mode_k<=2);
1527  if (mode_k==2) {
1528  double* vel_k = pot_set_pars[k].vel;
1529  double* pos_k = pot_set_pars[k].pos;
1530  pos_k[0] += vel_k[0]*dt_scale;
1531  pos_k[1] += vel_k[1]*dt_scale;
1532  pos_k[2] += vel_k[2]*dt_scale;
1533  }
1534  }
1535  }
1536 
1537 
1539 
1543  void calcMovePotAccFromPot(const double _time, const double* pos_pcm) {
1544  assert(pot_sets.size()==pot_set_pars.size());
1545  double t = _time*tscale;
1546 
1547  for (size_t k=0; k<pot_set_pars.size(); k++) {
1548  int mode_k = pot_set_pars[k].mode;
1549  assert(mode_k>=0||mode_k<=2);
1550  if (mode_k==2) {
1551  double* pos_k = pot_set_pars[k].pos; // galactic frame
1552  double* acc_pot = pot_set_pars[k].acc;
1553  for (size_t j=0; j<pot_set_pars.size(); j++) {
1554  if (k==j) continue;
1555  int mode_j = pot_set_pars[j].mode;
1556  double* pos_j = pot_set_pars[j].pos;
1557  double dx,dy,dz;
1558  if (mode_j==1) {
1559  // should be consistent in galactic frame
1560  dx = pos_k[0]-pos_j[0]-pos_pcm[0];
1561  dy = pos_k[1]-pos_j[1]-pos_pcm[1];
1562  dz = pos_k[2]-pos_j[2]-pos_pcm[2];
1563  }
1564  else {
1565  dx = pos_k[0]-pos_j[0];
1566  dy = pos_k[1]-pos_j[1];
1567  dz = pos_k[2]-pos_j[2];
1568  }
1569  double rxy= std::sqrt(dx*dx+dy*dy);
1570  //double phi= std::acos(dx/rxy);
1571  double phi= std::atan2(dy, dx);
1572  double sinphi = dy/rxy;
1573  double cosphi = dx/rxy;
1574  auto& pot_args = pot_sets[j].arguments;
1575  int npot = pot_sets[j].npot;
1576  double acc_rxy = calcRforce(rxy, dz, phi, t, npot, pot_args);
1577  double acc_z = calczforce(rxy, dz, phi, t, npot, pot_args);
1578 #if (defined GALPY_VERSION_1_7_9) || (defined GALPY_VERSION_1_7_1)
1579  double acc_phi = calcPhiforce(rxy, dz, phi, t, npot, pot_args);
1580 #else
1581  double acc_phi = calcphitorque(rxy, dz, phi, t, npot, pot_args);
1582 #endif
1583  if (rxy>0.0) {
1584  acc_pot[0] += (cosphi*acc_rxy - sinphi*acc_phi/rxy);
1585  acc_pot[1] += (sinphi*acc_rxy + cosphi*acc_phi/rxy);
1586  acc_pot[2] += acc_z;
1587  }
1588  }
1589  }
1590  }
1591  }
1592 
1593 
1595 
1603  void calcAccPot(double* acc, double& pot, const double _time, const double gm, const double* pos_g, const double* pos_l) {
1604  assert(pot_sets.size()==pot_set_pars.size());
1605  int nset = pot_sets.size();
1606  if (nset>0) {
1607  double t = _time*tscale;
1608 
1609  // galactic frame and rest frame of particle system
1610  double x[2] = {pos_g[0]*rscale, pos_l[0]*rscale};
1611  double y[2] = {pos_g[1]*rscale, pos_l[1]*rscale};
1612  double z[2] = {pos_g[2]*rscale, pos_l[2]*rscale};
1613 
1614  pot = 0;
1615  acc[0] = acc[1] = acc[2] = 0.0;
1616 
1617 
1618 
1619  for (int k=0; k<nset; k++) {
1620  int mode_k = pot_set_pars[k].mode;
1621  assert(mode_k>=0||mode_k<=2);
1622  int npot = pot_sets[k].npot;
1623  double* pos_k = pot_set_pars[k].pos;
1624  int i = (mode_k & 1); // get first bit to select frame (0: galactic; 1: rest)
1625  // frame is consistent
1626  double dx = x[i]-pos_k[0];
1627  double dy = y[i]-pos_k[1];
1628  double dz = z[i]-pos_k[2];
1629  double rxy= std::sqrt(dx*dx+dy*dy);
1630  double phi= std::atan2(dy, dx);
1631  double sinphi = dy/rxy;
1632  double cosphi = dx/rxy;
1633 
1634  auto& pot_args = pot_sets[k].arguments;
1635  double acc_rxy = calcRforce(rxy, dz, phi, t, npot, pot_args);
1636  double acc_z = calczforce(rxy, dz, phi, t, npot, pot_args);
1637 #if (defined GALPY_VERSION_1_7_9) || (defined GALPY_VERSION_1_7_1)
1638  double acc_phi = calcPhiforce(rxy, dz, phi, t, npot, pot_args);
1639 #else
1640  double acc_phi = calcphitorque(rxy, dz, phi, t, npot, pot_args);
1641 #endif
1642 // double pot_i = evaluatePotentials(rxy, dz, phi, t, npot, pot_args);
1643  double pot_i = evaluatePotentials(rxy, dz, npot, pot_args);
1644  double gm_pot = pot_set_pars[k].gm;
1645  if (rxy>0.0) {
1646  assert(!std::isinf(acc_rxy));
1647  assert(!std::isnan(acc_rxy));
1648  assert(!std::isinf(acc_phi));
1649  assert(!std::isnan(acc_phi));
1650  assert(!std::isinf(pot));
1651  assert(!std::isnan(pot));
1652  pot += pot_i;
1653  double acc_x = (cosphi*acc_rxy - sinphi*acc_phi/rxy);
1654  double acc_y = (sinphi*acc_rxy + cosphi*acc_phi/rxy);
1655  acc[0] += acc_x;
1656  acc[1] += acc_y;
1657  acc[2] += acc_z;
1658  if (mode_k==2) {
1659  double* acc_pot = pot_set_pars[k].acc;
1660  acc_pot[0] -= gm*acc_x/gm_pot; // anti-acceleration to potential set origin
1661  acc_pot[1] -= gm*acc_y/gm_pot;
1662  acc_pot[2] -= gm*acc_z/gm_pot;
1663  }
1664  }
1665  }
1666  pot /= pscale;
1667  acc[0] /= fscale;
1668  acc[1] /= fscale;
1669  acc[2] /= fscale;
1670 
1671  //pot = acc[0]*x+acc[1]*y+acc[2]*z;
1672  }
1673  else {
1674  acc[0] = acc[1] = acc[2] = pot = 0.0;
1675  }
1676  }
1677 
1679 
1693  void writeDataToFile(const std::string& _filename) {
1694  std::ofstream fout;
1695  fout.open(_filename.c_str(), std::ifstream::out);
1696  if (!fout.is_open()) {
1697  std::cerr<<"Error: Galpy potential parameter file to write, "<<_filename<<", cannot be open!"<<std::endl;
1698  abort();
1699  }
1700  fout<<std::setprecision(14);
1701 
1702  std::size_t n_set = pot_set_pars.size();
1703  std::size_t n_pot = pot_type.size();
1704  std::size_t n_pot_offset = pot_type_offset.size();
1705  std::size_t n_arg = pot_args.size();
1706  std::size_t n_arg_offset = pot_args_offset.size();
1707  std::size_t n_change = change_args.size();
1708  std::size_t n_change_offset = change_args_offset.size();
1709  fout<<time<<" "
1710  <<n_set<<" "
1711  <<std::endl;
1712  for (std::size_t i=0; i<n_set; i++) {
1713  pot_set_pars[i].writeData(fout);
1714  fout<<std::endl;
1715  }
1716  fout<<std::endl<<n_pot_offset<<" ";
1717  for (std::size_t i=0; i<n_pot_offset; i++) {
1718  fout<<pot_type_offset[i]<<" ";
1719  }
1720  fout<<std::endl<<n_pot<<" ";
1721  for (std::size_t i=0; i<n_pot; i++) {
1722  fout<<pot_type[i]<<" ";
1723  }
1724  fout<<std::endl<<n_arg_offset<<" ";
1725  for (std::size_t i=0; i<n_arg_offset; i++) {
1726  fout<<pot_args_offset[i]<<" ";
1727  }
1728  fout<<std::endl<<n_arg<<" ";
1729  for (std::size_t i=0; i<n_arg; i++) {
1730  fout<<pot_args[i]<<" ";
1731  }
1732  fout<<std::endl<<n_change_offset<<" ";
1733  for (std::size_t i=0; i<n_change_offset; i++) {
1734  fout<<change_args_offset[i]<<" ";
1735  }
1736  fout<<std::endl<<n_change<<" ";
1737  for (std::size_t i=0; i<n_change; i++) {
1738  change_args[i].writeData(fout);
1739  }
1740  fout<<std::endl;
1741 
1742  fout.close();
1743  }
1744 
1746  /*
1747  @param[in] _filename: file to save data
1748  @param[in] _print_flag: if true, print reading parameters
1749  */
1750  void readDataFromFile(const std::string& _filename, const bool _print_flag) {
1751 
1752  std::size_t n_set, n_pot, n_pot_offset, n_arg, n_arg_offset, n_change, n_change_offset;
1753 
1754 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1755  int my_rank = PS::Comm::getRank();
1756 
1757  if (my_rank==0) {
1758 #endif
1759  std::ifstream fin;
1760  fin.open(_filename.c_str(), std::ifstream::in);
1761  if (!fin.is_open()) {
1762  std::cerr<<"Error: Galpy configure file for restart "<<_filename.c_str()<<" cannot be open!"<<std::endl;
1763  abort();
1764  }
1765 
1766  fin>>time>>n_set;
1767  eofCheck(fin, "time and number of potential sets");
1768  geZeroCheck(n_set, "number of potential sets");
1769  pot_set_pars.resize(n_set);
1770  for (std::size_t i=0; i<n_set; i++) {
1771  pot_set_pars[i].readData(fin);
1772  eofCheck(fin, "potential set parameter");
1773  }
1774  fin>>n_pot_offset;
1775  assert(n_pot_offset==n_set+1);
1776  pot_type_offset.resize(n_pot_offset);
1777  for (std::size_t i=0; i<n_pot_offset; i++) {
1778  fin>>pot_type_offset[i];
1779  }
1780  eofCheck(fin, "potential type array offset");
1781 
1782  fin>>n_pot;
1783  geZeroCheck(n_pot, "number of potential types");
1784  pot_type.resize(n_pot);
1785  for (std::size_t i=0; i<n_pot; i++) {
1786  fin>>pot_type[i];
1787  }
1788  eofCheck(fin, "potential type");
1789 
1790  fin>>n_arg_offset;
1791  pot_args_offset.resize(n_arg_offset);
1792  assert(n_arg_offset==n_set+1);
1793  for (std::size_t i=0; i<n_arg_offset; i++) {
1794  fin>>pot_args_offset[i];
1795  }
1796  eofCheck(fin, "potential argument array offset");
1797 
1798  fin>>n_arg;
1799  geZeroCheck(n_arg, "number of potential arguments");
1800  pot_args.resize(n_arg);
1801  for (std::size_t i=0; i<n_arg; i++) {
1802  fin>>pot_args[i];
1803  }
1804  eofCheck(fin, "potential arguments");
1805 
1806  fin>>n_change_offset;
1807  change_args_offset.resize(n_change_offset);
1808  assert(n_change_offset==n_set+1);
1809  for (std::size_t i=0; i<n_change_offset; i++) {
1810  fin>>change_args_offset[i];
1811  }
1812  eofCheck(fin, "changing argument offset");
1813 
1814  fin>>n_change;
1815  geZeroCheck(n_change, "number of changing arguments");
1816  change_args.resize(n_change);
1817  for (std::size_t i=0; i<n_change; i++) {
1818  change_args[i].readData(fin);
1819  }
1820  eofCheck(fin, "changing arguments");
1821 
1822  fin.close();
1823 
1824 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1825  }
1826 
1827  broadcastDataMPI();
1828 #endif
1829 
1830  if (_print_flag) printData(std::cout);
1831  }
1832 
1834 
1843  bool evolveChangingArguments(const double& _time, const bool _print_flag) {
1844  if (change_args.size()>0) {
1845  double dt = _time - time;
1846  if (dt>0) {
1847  if (_print_flag) std::cout<<"Galpy update argument, time [Galpy unit]: "<<_time<<" dt: "<<dt<<std::endl;
1848  std::size_t n_set = change_args_offset.size()-1;
1849  for (std::size_t i=0; i<n_set; i++) {
1850  for (int k=change_args_offset[i]; k<change_args_offset[i+1]; k++) {
1851  int index_local = change_args[k].index;
1852  double* change_arg_ptr = NULL;
1853  if (index_local==-1) // change gm instead of arguments
1854  change_arg_ptr = &(pot_set_pars[i].gm);
1855  else {
1856  int index = index_local + pot_args_offset[i]; // change args index counting for individual set, need to obtain the global index of pot_args.
1857  change_arg_ptr = &(pot_args[index]);
1858  }
1859  if (change_args[k].mode==1) {
1860  *change_arg_ptr += change_args[k].rate*dt;
1861  if (_print_flag) std::cout<<"Set "<<i+1<<" Index: "<<change_args[k].index<<" "<<" new_argument(linear): "<<*change_arg_ptr;
1862  }
1863  else if(change_args[k].mode==2) {
1864  *change_arg_ptr *= std::exp(change_args[k].rate*dt);
1865  if (_print_flag) std::cout<<"Set "<<i+1<<" Index: "<<change_args[k].index<<" "<<" new_argument(exp): "<<*change_arg_ptr;
1866  }
1867  if (_print_flag) std::cout<<std::endl;
1868  }
1869  }
1870  time = _time;
1871  return true;
1872  }
1873  }
1874  time = _time;
1875  return false;
1876  }
1877 
1878  void freePotentialArgs() {
1879  if (!pot_sets.empty()) {
1880  for (size_t i=0; i<pot_sets.size(); i++) pot_sets[i].clear();
1881  pot_sets.resize(0);
1882  }
1883  }
1884 
1885  void clear() {
1886  freePotentialArgs();
1887  update_time = 0;
1888  if (fconf.is_open()) fconf.close();
1889  }
1890 
1891  ~GalpyManager() {
1892  clear();
1893  }
1894 };
1895 
1896 
PotentialSetPar::gm
double gm
Definition: galpy_interface.h:546
PotentialSet::npot
int npot
Definition: galpy_interface.h:603
Status::ParticleCM::vel
PS::F64vec vel
Definition: status.hpp:21
GalpyManager::pscale
double pscale
Definition: galpy_interface.h:753
FRWModel::a
double a
Definition: galpy_interface.h:185
GalpyManager::set_parfile
std::string set_parfile
Definition: galpy_interface.h:758
FileHeader::n_body
long long int n_body
Definition: io.hpp:213
GalpyManager::readDataFromFile
void readDataFromFile(const std::string &_filename, const bool _print_flag)
read configure data for restart
Definition: galpy_interface.h:1750
Status::pcm
struct Status::ParticleCM pcm
FileHeaderNoOffset::nfile
long long int nfile
Definition: format_transfer.cxx:47
galpy_help.n_pot
n_pot
Definition: galpy_help.py:306
FPSoft
Definition: soft_ptcl.hpp:26
galpy_help.savePotTypeArg
def savePotTypeArg(config_filename, n_pot, pot_type, pot_arg)
Definition: galpy_help.py:45
ChangeArgument::mode
int mode
Definition: galpy_interface.h:643
SystemSoft
PS::ParticleSystem< FPSoft > SystemSoft
Definition: format_transfer.cxx:8
GalpyManager::fscale
double fscale
Definition: galpy_interface.h:752
main
int main(int argc, char *argv[])
Definition: format_transfer.cxx:92
GalpyManager::printReference
static void printReference(std::ostream &fout, const int offset=4)
print reference to cite
Definition: galpy_interface.h:1486
GalpyManager::pot_set_pars
std::vector< PotentialSetPar > pot_set_pars
Definition: galpy_interface.h:745
IOParams::value
Type value
Definition: io.hpp:43
PotentialSetPar::mode
int mode
Definition: galpy_interface.h:545
GalpyManager::writePotentialPars
void writePotentialPars(const std::string &_filename, const double &_system_time)
write potential parameters for restart
Definition: galpy_interface.h:1474
FileHeaderNoOffset::writeAscii
void writeAscii(FILE *fp) const
Definition: format_transfer.cxx:70
FileHeaderNoOffset::FileHeaderNoOffset
FileHeaderNoOffset(const long long int ni, const long long int n, const double t)
Definition: format_transfer.cxx:54
MWPotentialEvolve::m_vir_halo
double m_vir_halo
Definition: galpy_interface.h:299
galpy_help.listPot
def listPot()
Definition: galpy_help.py:140
ChangeArgument
changing argument parameter
Definition: galpy_interface.h:641
GalpyManager::GalpyManager
GalpyManager()
Definition: galpy_interface.h:762
MWPotentialEvolve::gm_halo
double gm_halo
Definition: galpy_interface.h:313
galpy_pot_movie.fp
fp
Definition: galpy_pot_movie.py:131
soft_ptcl.hpp
MWPotentialEvolve::alpha_bulge
double alpha_bulge
Definition: galpy_interface.h:306
MWPotentialEvolve::initialFromFile
void initialFromFile(const std::string &_filename, const double &_time)
initial potential parameter from a file
Definition: galpy_interface.h:343
FileHeader::time
double time
Definition: io.hpp:214
GalpyManager::updatePotential
void updatePotential(const double &_system_time, const bool _print_flag)
Update Potential types and arguments based on time.
Definition: galpy_interface.h:1443
PotentialSet::~PotentialSet
~PotentialSet()
Definition: galpy_interface.h:635
IOParamsContainer::readAscii
void readAscii(FILE *_fin)
Definition: io.hpp:112
galpy_help.printCPot
def printCPot(pot_name)
Definition: galpy_help.py:217
FPSoftWriteArtificial
redefine readAscii to read group data format, write in I64 format
Definition: format_transfer.cxx:11
IOParamsGalpy::rscale
IOParams< double > rscale
Definition: galpy_interface.h:22
galpy_help.printPotTypeArg
def printPotTypeArg(pot_name, pot_module, pot_instance, print_front_offset=0, print_long_list=False)
Definition: galpy_help.py:70
galpy_help.getPotInstance
def getPotInstance(pot_name)
Definition: galpy_help.py:7
IOParamsContainer::printHelp
void printHelp(std::ostream &os, const int _offset_short_key=2, const int _offset_long_key=1, const int _width_key=23) const
Definition: io.hpp:198
FRWModel::dt_max
double dt_max
Definition: galpy_interface.h:184
FRWModel::getDensity
double getDensity(const double _a) const
Definition: galpy_interface.h:209
galpy_pot_movie.dt
dt
Definition: galpy_pot_movie.py:134
GalpyManager::pot_type_offset
std::vector< int > pot_type_offset
Definition: galpy_interface.h:737
FRWModel::omega_radiation
double omega_radiation
Definition: galpy_interface.h:189
FRWModel
FRW cosmological model for calculating scale factor (redshift)
Definition: galpy_interface.h:181
GalpyManager::initial
void initial(const IOParamsGalpy &_input, const double _time, const std::string _conf_name, const bool _restart_flag, const bool _print_flag=false)
initialization function
Definition: galpy_interface.h:852
FRWModel::updateAwithDtmin
void updateAwithDtmin(const double _dt)
update scale factor with maximum time step dt_max (first order Euler method)
Definition: galpy_interface.h:247
FRWModel::getEnergyDensity
double getEnergyDensity(const double _a) const
Definition: galpy_interface.h:201
FileHeaderNoOffset::n_body
long long int n_body
Definition: format_transfer.cxx:48
MWPotentialEvolve::c_halo
double c_halo
Definition: galpy_interface.h:300
MWPotentialEvolve::init
struct MWPotentialEvolve::@6 init
galpy_help.printPotTitle
def printPotTitle()
Definition: galpy_help.py:133
FRWModel::getMatterDensity
double getMatterDensity(const double _a) const
Definition: galpy_interface.h:193
GalpyManager::updateMWPotentialEvolve
bool updateMWPotentialEvolve(const double &_system_time, const bool _print_flag, const bool _initial_flag)
Update MWpotential evolve (Gomez et al. 2010)
Definition: galpy_interface.h:1067
MWPotentialEvolve::rho_c_halo
double rho_c_halo
Definition: galpy_interface.h:312
MWPotentialEvolve::frw
FRWModel frw
Definition: galpy_interface.h:297
PotentialSet::arguments
struct potentialArg * arguments
Definition: galpy_interface.h:604
PotentialSetPar::writeData
void writeData(std::ostream &fout)
print data in one line
Definition: galpy_interface.h:554
GalpyManager::calcAccPot
void calcAccPot(double *acc, double &pot, const double _time, const double gm, const double *pos_g, const double *pos_l)
calculate acceleration and potential at give position
Definition: galpy_interface.h:1603
galpy_help.listCPot
def listCPot()
Definition: galpy_help.py:239
IOParamsGalpy::print_flag
bool print_flag
Definition: galpy_interface.h:28
PotentialSetPar::pos
double pos[3]
Definition: galpy_interface.h:547
IOParams::print
void print(std::ostream &os) const
Definition: io.hpp:54
FileHeaderNoOffset
old version without c.m. pos and vel offset
Definition: format_transfer.cxx:45
MWPotentialEvolve::ac_halo
double ac_halo
Definition: galpy_interface.h:301
GalpyManager::kickMovePot
void kickMovePot(const double dt)
kick velocity of moving potential (mode 2)
Definition: galpy_interface.h:1503
BSEManager::getBSEName
static std::string getBSEName()
Definition: bse_interface.h:1198
MWPotentialEvolve::now
struct MWPotentialEvolve::@7 now
galpy_help.usage
def usage()
Definition: galpy_help.py:271
MWPotentialEvolve::ra_disk
double ra_disk
Definition: galpy_interface.h:303
GalpyManager
A class to manager the API to Galpy.
Definition: galpy_interface.h:658
MWPotentialEvolve::calcMWPotentialParameters
void calcMWPotentialParameters(const double &_time)
calculate MW potential based on the input time
Definition: galpy_interface.h:498
FileHeader::readBinary
int readBinary(FILE *fp)
Definition: io.hpp:269
Status
class for measure the status of the system
Definition: status.hpp:7
IOParamsGalpy::input_par_store
IOParamsContainer input_par_store
Definition: galpy_interface.h:17
PIKG::S64
int64_t S64
Definition: pikg_vector.hpp:23
PotentialSetPar::vel
double vel[3]
Definition: galpy_interface.h:548
MWPotentialEvolve::gm_disk
double gm_disk
Definition: galpy_interface.h:315
IOParams::key
const char * key
Definition: io.hpp:44
FileHeader::writeAscii
void writeAscii(FILE *fp) const
Definition: io.hpp:265
GalpyManager::time
double time
Definition: galpy_interface.h:741
FPSoft::acc
PS::F64vec acc
Definition: soft_ptcl.hpp:28
io.hpp
GalpyManager::mw_evolve
MWPotentialEvolve mw_evolve
Definition: galpy_interface.h:760
GalpyManager::update_time
double update_time
Definition: galpy_interface.h:747
GalpyManager::writeDataToFile
void writeDataToFile(const std::string &_filename)
write data for restart
Definition: galpy_interface.h:1693
FRWModel::omega_matter
double omega_matter
Definition: galpy_interface.h:190
GalpyManager::freePotentialArgs
void freePotentialArgs()
Definition: galpy_interface.h:1878
PotentialSet::generatePotentialArgs
void generatePotentialArgs(const int _npot, int *_type, double *_arg)
generate Galpy potential arguments
Definition: galpy_interface.h:614
MWPotentialEvolve
data for MWPotential evolve
Definition: galpy_interface.h:277
PotentialSet::PotentialSet
PotentialSet()
Definition: galpy_interface.h:606
PotentialSetPar::clear
void clear()
Definition: galpy_interface.h:591
GalpyManager::printData
void printData(std::ostream &fout)
print current potential data
Definition: galpy_interface.h:767
galpy_help.printSpliter
def printSpliter()
Definition: galpy_help.py:137
FileHeaderNoOffset::writeBinary
void writeBinary(FILE *fp) const
Definition: format_transfer.cxx:84
FRWModel::dt_scale
double dt_scale
Definition: galpy_interface.h:191
MWPotentialEvolve::writeDataToFile
void writeDataToFile(const std::string &_filename)
write data for restart
Definition: galpy_interface.h:466
IOParamsGalpy::setBovyUnit
void setBovyUnit(const bool print_flag=true)
set Bovy unit scale factor (solar distance and velocityr referring to the Galactic center) converge f...
Definition: galpy_interface.h:156
MWPotentialEvolve::rcut_bulge
double rcut_bulge
Definition: galpy_interface.h:307
FileHeaderNoOffset::FileHeaderNoOffset
FileHeaderNoOffset()
Definition: format_transfer.cxx:50
FRWModel::omega_energy
double omega_energy
Definition: galpy_interface.h:188
GalpyManager::clear
void clear()
Definition: galpy_interface.h:1885
MWPotentialEvolve::rho_bulge
double rho_bulge
Definition: galpy_interface.h:305
GalpyManager::pot_type
std::vector< int > pot_type
Definition: galpy_interface.h:738
PotentialSetPar::acc
double acc[3]
Definition: galpy_interface.h:549
Ptcl::r_search
PS::F64 r_search
Definition: ptcl.hpp:38
GalpyManager::rscale
double rscale
Definition: galpy_interface.h:749
MWPotentialEvolve::m_disk
double m_disk
Definition: galpy_interface.h:302
Ptcl::changeover
ChangeOver changeover
Definition: ptcl.hpp:41
IOParamsGalpy::config_filename
IOParams< std::string > config_filename
Definition: galpy_interface.h:20
GalpyManager::change_args
std::vector< ChangeArgument > change_args
Definition: galpy_interface.h:742
status.hpp
MWPotentialEvolve::rs_halo
double rs_halo
Definition: galpy_interface.h:314
FRWModel::calcH0Myr
void calcH0Myr()
Definition: galpy_interface.h:213
MWPotentialEvolve::grho_bulge
double grho_bulge
Definition: galpy_interface.h:318
Status::ParticleCM::pos
PS::F64vec pos
Definition: status.hpp:20
PotentialSet::clear
void clear()
Definition: galpy_interface.h:625
ChangeArgument::rate
double rate
Definition: galpy_interface.h:644
FileHeaderNoOffset::time
double time
Definition: format_transfer.cxx:49
GalpyManager::~GalpyManager
~GalpyManager()
Definition: galpy_interface.h:1891
FileHeaderNoOffset::readAscii
int readAscii(FILE *fp)
Definition: format_transfer.cxx:60
ChangeArgument::index
int index
Definition: galpy_interface.h:642
IOParams< std::string >
GalpyManager::pot_sets
std::vector< PotentialSet > pot_sets
Definition: galpy_interface.h:746
Status::calcAndShiftCenterOfMass
void calcAndShiftCenterOfMass(Tsoft *_tsys, const PS::S64 _n, const int _mode=3, const bool initial_flag=false)
calculate the center of system and shift particle systems to center frame
Definition: status.hpp:228
FRWModel::H0_myr
double H0_myr
Definition: galpy_interface.h:187
FPSoftWriteArtificial::readAscii
void readAscii(FILE *_fin)
Definition: format_transfer.cxx:24
IOParamsGalpy::read
int read(int argc, char *argv[], const int opt_used_pre=0)
reading parameters from GNU option API
Definition: galpy_interface.h:49
IOParamsGalpy
IO parameters for Galpy manager.
Definition: galpy_interface.h:15
IOParamsGalpy::type_args
IOParams< std::string > type_args
Definition: galpy_interface.h:18
ChangeArgument::readData
void readData(std::ifstream &fin)
print data in one line
Definition: galpy_interface.h:652
FileHeader::nfile
long long int nfile
Definition: io.hpp:212
ArtificialParticleInformation::readAscii
void readAscii(FILE *_fin)
read class data with ASCII format
Definition: artificial_particles.hpp:203
galpy_help.pot_module
pot_module
Definition: galpy_help.py:303
IOParamsGalpy::pre_define_type
IOParams< std::string > pre_define_type
Definition: galpy_interface.h:19
FRWModel::getRadiationDensity
double getRadiationDensity(const double _a) const
Definition: galpy_interface.h:197
GroupDataDeliver::artificial
ArtificialParticleInformation artificial
Definition: ptcl.hpp:18
IOParamsGalpy::IOParamsGalpy
IOParamsGalpy()
Definition: galpy_interface.h:30
GalpyManager::updatePotentialSet
void updatePotentialSet()
generate potentail args
Definition: galpy_interface.h:969
FRWModel::updateA
void updateA(const double _dt)
update scale factor by time step _dt (first order Euler method)
Definition: galpy_interface.h:231
MWPotentialEvolve::r_vir_halo
double r_vir_halo
Definition: galpy_interface.h:311
GalpyManager::fconf
std::ifstream fconf
Definition: galpy_interface.h:756
GalpyManager::set_name
std::string set_name
Definition: galpy_interface.h:757
GalpyManager::pot_args_offset
std::vector< int > pot_args_offset
Definition: galpy_interface.h:739
SystemSoftWriteArtificial
PS::ParticleSystem< FPSoftWriteArtificial > SystemSoftWriteArtificial
Definition: format_transfer.cxx:90
GalpyManager::gmscale
double gmscale
Definition: galpy_interface.h:754
galpy_pot_movie.nstep
nstep
Definition: galpy_pot_movie.py:134
ChangeOver::readAscii
void readAscii(FILE *_fin)
read class data to file with binary format
Definition: changeover.hpp:168
FRWModel::calcAdot
double calcAdot(const double _a) const
calculate da/dt
Definition: galpy_interface.h:223
FileHeaderNoOffset::readBinary
int readBinary(FILE *fp)
Definition: format_transfer.cxx:74
PotentialSetPar::setOrigin
void setOrigin(const int _mode, const double _gm=0, const double *_pos=NULL, const double *_vel=NULL, const double *_acc=NULL)
set position, velocity and acceleration
Definition: galpy_interface.h:570
GalpyManager::driftMovePot
void driftMovePot(const double dt)
drift position of moving potential (mode 2)
Definition: galpy_interface.h:1522
FRWModel::getH
double getH() const
get Hubble constant at current a
Definition: galpy_interface.h:267
PotentialSet
set of Galpy potentials sharing the same position and velocity of the origin
Definition: galpy_interface.h:602
PotentialSetPar
mode, position, velocity and acceleration for set of Galpy potentials
Definition: galpy_interface.h:544
GalpyManager::resetPotAcc
void resetPotAcc()
reset acceleraltion of potential
Definition: galpy_interface.h:1492
FRWModel::getCurvatureDensity
double getCurvatureDensity() const
Definition: galpy_interface.h:205
FRWModel::time
double time
Definition: galpy_interface.h:183
GalpyManager::updatePotentialFromFile
bool updatePotentialFromFile(const double &_system_time, const bool _print_flag)
Update types and arguments from configure file if update_time <= particle system time.
Definition: galpy_interface.h:1111
IOParamsGalpy::vscale
IOParams< double > vscale
Definition: galpy_interface.h:24
FileHeader
Definition: io.hpp:210
GalpyManager::vscale
double vscale
Definition: galpy_interface.h:750
GalpyManager::change_args_offset
std::vector< int > change_args_offset
Definition: galpy_interface.h:743
GalpyManager::addPotentialFromString
bool addPotentialFromString(const std::string &_type_args, const bool _reset_flag, const bool _print_flag)
Update types and arguments from type-args string.
Definition: galpy_interface.h:986
IOParamsContainer
IO Params container.
Definition: io.hpp:82
Ptcl::group_data
GroupDataDeliver group_data
Definition: ptcl.hpp:40
GalpyManager::pot_args
std::vector< double > pot_args
Definition: galpy_interface.h:740
MWPotentialEvolve::readDataFromFile
void readDataFromFile(const std::string &_filename, const double &_time)
read potential parameter from a file
Definition: galpy_interface.h:408
PotentialSetPar::PotentialSetPar
PotentialSetPar()
Definition: galpy_interface.h:551
ChangeArgument::writeData
void writeData(std::ostream &fout)
print data in one line
Definition: galpy_interface.h:647
GalpyManager::calcMovePotAccFromPot
void calcMovePotAccFromPot(const double _time, const double *pos_pcm)
calculate acceleration for moving potentials from other potentials
Definition: galpy_interface.h:1543
ParticleBase::readAscii
void readAscii(FILE *fp)
read class data with ASCII format
Definition: particle_base.hpp:161
GalpyManager::evolveChangingArguments
bool evolveChangingArguments(const double &_time, const bool _print_flag)
calculate new arguments for changing potentials
Definition: galpy_interface.h:1843
PotentialSetPar::readData
void readData(std::ifstream &fin)
Definition: galpy_interface.h:562
FRWModel::H0
double H0
Definition: galpy_interface.h:186
MWPotentialEvolve::rb_disk
double rb_disk
Definition: galpy_interface.h:304
GalpyManager::tscale
double tscale
Definition: galpy_interface.h:751