jobovy
commented
4 years ago Thanks for raising this issue, this is very interesting! (we discussed your code in one of my group meetings a few weeks ago!). This is definitely something I'd like to support and a goal I had started working towards recently, so I'd be happy to help out with any changes to galpy that might be necessary. The good news is that I think this is relatively easy to do.
The current situation is as follows:
-
A few months ago, I changed the C code so it compiles to a single shared library
libgalpy.so; this now makes it possible to link to it with-lgalpy, so that's all good (note that I've never tried it out in practice though...) -
The
galpy_potentials.hfile includes all of the potentials and functions to evaluate forces for arbitrary combinations of potentials, in particularcalcRforce,calczforce, andcalcphiforcewhich give the forces in cylindrical coordinates (no other coordinate system is defined, so you'd have to write conversions yourself [or we could add them togalpy]) as a function ofdouble R, double Z, double phi, double t, int nargs, struct potentialArg * potentialArgs, which are the cylindrical coordinates, time, and thennargsandpotentialArgswhich define the potential. So including thegalpy_potentials.hfile will give you access to the force evaluation functions you need -
Now, you want to get
nargsandpotentialArgsto define the potential. The waygalpypotentials are transferred to C is simply as two lists:- a list of integers that represent the type of each potential (e.g., '0' is a
LogarithmicPotential, '5' aMiyamotoNagaiPotential):int * pot_type - a list of doubles that gives the parameters of all of the potentials:
double * pot_args - I will explain below how you can get these two lists through the Python interface
nargsis just the number of potentials, so something you know...
- a list of integers that represent the type of each potential (e.g., '0' is a
-
To transform these lists to
potentialArgsyou needparse_leapFuncArgs_Full, which is defined ingalpy/orbit/orbit_c_ext/integrateFullOrbit.h(to make linking togalpyeasier, we might want to move that declaration togalpy_potentials.h). An example of how to use this is ingalpy/orbit/orbit_c_ext/integrateFullOrbit.c, here's a simplified version (the code there is a little more complicated, because it sets of N copies of the potential to use in multi-threading)struct potentialArg * potentialArgs= (struct potentialArg *) malloc ( npot * sizeof (struct potentialArg) ); parse_leapFuncArgs_Full(npot,potentialArgs,&pot_type,&pot_args);after running that code,
nargs=npotand you havepotentialArgsto use in the force evaluation functions. At the end of your code, you want to free the pointers withfree_potentialArgs(npot,potentialArgs); free(potentialArgs);because
potentialArgscontains a bunch of allocated memory that needs to be carefully freed. I think every code that wants to link togalpyhas to call these functions themselves when they initialize their simulation. -
To get the
pot_typeandpot_argslists, you can usefrom galpy.orbit.integrateFullOrbit import _parse_pot npot, pot_type, pot_args= _parse_pot(pot)For example,
from galpy.potential import MWPotential2014 npot, pot_type, pot_args= _parse_pot(MWPotential2014) print(npot,pot_type,pot_args) # 3 [15 5 9] [0.0299946 1.8 0.2375 0.7574802 0.375 0.035 4.85223053 2. ]Again, we might want to make the
_parse_potfunction more visible to make it easier to link togalpy.
Currently, galpy does not get installed to an obvious place to link to (e.g., /usr/local, I think the library is just put among the Python code in your .../site-packages/... directory and the header files aren't copied at all I think), so I think the way for you to start is to:
- Copy the header files by hand to, e.g.,
/usr/local/include, copy the library to/usr/local/lib(you can get the library by just compiling withpython setup.py build_ext) - Write some interface between your code and
galpy's C version that takes thenpot,pot_type,pot_argsdefinition of the potential and sets uppotentialArgs. - Then use the
calcRforce, etc., functions to evaluate the external forces in your code. I would probably add thegalpybindings in your code using preprocessor directives to only include these bindings when wanted.
I'm happy to help out and if this seems to be working, I can also work on getting galpy to install the header and library files to an installation directory and to restructure some of the code to make it easier to link to (e.g., put all functions in galpy_potentials.h, maybe rename to galpy.h).
Final note: the calcRforce etc. functions optionally take the velocity in cylindrical coordinates for velocity dependent forces, like dynamical friction. If you want to support those, also provide the velocity when you call calcRforce (velocities are optional using some preprocessor magic...).
Also cc-ing @webbjj, because he may be interested in this conversation.
lwang-astro
webbjj
Recently I have published a new N-body code, PeTar, aiming for simulating star clusters where many close encounters and binaries exist. The code is written in c++. I'd like to implement the external potentials as a tidal field for the clusters. I think galpy is a nice choice since it supports many kinds of potentials and also have c/c++ interface. I have tried to search the document and check the galpy_potential.h file, but not yet fully understand how to correctly implement the galpy as a c/c++ library. The major part of the documents seem to be for python interface. Is any test example writting in c/c++? That will be very helpful for me to understand the code.
What I'd like to have is a c/c++ API function to provide the 3D forces (x, y, z) for a given position. Besides, users can combine a group of potential functions like in the python interface.