Questions about the c/c++ API

[lwang-astro]

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.

[jobovy]

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.h file includes all of the potentials and functions to evaluate forces for arbitrary combinations of potentials, in particular calcRforce, calczforce, and calcphiforce which 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 to galpy]) as a function of double R, double Z, double phi, double t, int nargs, struct potentialArg * potentialArgs, which are the cylindrical coordinates, time, and then nargs and potentialArgs which define the potential. So including the galpy_potentials.h file will give you access to the force evaluation functions you need

• Now, you want to get nargs and potentialArgs to define the potential. The way galpy potentials 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' a MiyamotoNagaiPotential): 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
  • nargs is just the number of potentials, so something you know...

• To transform these lists to potentialArgs you need parse_leapFuncArgs_Full, which is defined in galpy/orbit/orbit_c_ext/integrateFullOrbit.h (to make linking to galpy easier, we might want to move that declaration to galpy_potentials.h). An example of how to use this is in galpy/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=npot and you have potentialArgs to use in the force evaluation functions. At the end of your code, you want to free the pointers with

  free_potentialArgs(npot,potentialArgs);
  free(potentialArgs);
  

  because potentialArgs contains a bunch of allocated memory that needs to be carefully freed. I think every code that wants to link to galpy has to call these functions themselves when they initialize their simulation.

• To get the pot_type and pot_args lists, you can use

  from 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_pot function more visible to make it easier to link to galpy.

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 with python setup.py build_ext)
• Write some interface between your code and galpy's C version that takes the npot,pot_type,pot_args definition of the potential and sets up potentialArgs.
• Then use the calcRforce, etc., functions to evaluate the external forces in your code. I would probably add the galpy bindings 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]

Thank you for the very detailed description! I am happy that you are interesting to help me to implement the galpy API for petar. I have managed to write a simple interface to use force calculation function. I tested the MWPotential2014 and the result is consistent with the python output. I still have a few questions.

1. The unit.
   In python interface, there is an unit converter. I am not sure how is the case in the c code. Basically, I prefer to get the unit for acceleration as pc/Myr^2 and potential as pc^2/Myr^2 in default and make possible options for users to choose different units.

   Besides, for the argument to calcRforce and calczforce, I guess the time unit is Myr and distance unit is 8kpc? Is it possible to make all units self-consistent for the input and output, i.e., all units are limited to be a combination of (Myr, pc, Msun), such as distance: pc; time: Myr; velocity: pc/Myr; acceleration pc/Myr^2; potential: pc^2/Myr^2? I think this would reduce the confusion and make the petar interface simpler (Myr, pc and Msun are one of the basic unit sets in Petar).

   A related question is that if the input distance unit is 8pc and output acceleration unit is pc/Myr^2, it means that there is an internal conversion of distance unit in the force calculation function. Then for the arguments of the potential models (pot_args), what should be the proper units?

2. The potential calculation
   In integrateFullOrbit.h, the force calculation is available, how about the potential value for a given position?

3. The list of pot_type and pot_args
   Right now I can manage to generate the pot_type and pot_args in the c++ interface. But without checking the python commander, is it possible to get the help information that show the corresponding type and definition of arguments for a potential? For exmple, some help information in a c function or a text file like:
   Name | Type | argument
   Miyamoto | 5 | A: def. B: def. ...
   ...

   For the PeTar users who want to directly set the type and argments of potentials without checking python galpy, this help information (will be printed in the petar commander with option "-h") can be very convenient. Right now I can get the type from the integrateFullOrbit.c, but the definitions of arguments are not shown there.

4. The library file
   From the pip installed galpy, I did not find the libgalpy.so. But I manage to compile the library myself by including all .c files in galpy directory and it seems working correctly. I think pip probably only provide the source but not the compiled library, since the latter can depends on the machine configuration. I also prefer to compile the lib manually, since in some supercomputer, the pip may not work and the galpy source code has to be downloaded and installed. Also there is no root permission to install the lib in /user/local.

5. OpenMP
   When I compile the galpy library, there is warning of OpenMP, how the multi threads are actually used in galpy? My idea is to add openmp parallel loop outside force calculation function, like
   #pragma omp parallel for
for (int i=0; i<n_particle; i++) {
calcAccPot(acc[i], pot[i], time, pos[i]);
}
   Where calcAccPot is the function to use calcRforce and calczforce to obtain the acceleration. In such case, will be some confliction existing for using the multi threads?

[lwang-astro]

By the way, my current implementation of the interface can be found in Git:Petar.

[jobovy]

Thank you for the very detailed description! I am happy that you are interesting to help me to implement the galpy API for petar. I have managed to write a simple interface to use force calculation function. I tested the MWPotential2014 and the result is consistent with the python output. I still have a few questions.

1. The unit.
   In python interface, there is an unit converter. I am not sure how is the case in the c code. Basically, I prefer to get the unit for acceleration as pc/Myr^2 and potential as pc^2/Myr^2 in default and make possible options for users to choose different units.
   Besides, for the argument to calcRforce and calczforce, I guess the time unit is Myr and distance unit is 8kpc? Is it possible to make all units self-consistent for the input and output, i.e., all units are limited to be a combination of (Myr, pc, Msun), such as distance: pc; time: Myr; velocity: pc/Myr; acceleration pc/Myr^2; potential: pc^2/Myr^2? I think this would reduce the confusion and make the petar interface simpler (Myr, pc and Msun are one of the basic unit sets in Petar).
   A related question is that if the input distance unit is 8pc and output acceleration unit is pc/Myr^2, it means that there is an internal conversion of distance unit in the force calculation function. Then for the arguments of the potential models (pot_args), what should be the proper units?

Glad to hear that you got a version to work!

The units of inputs and outputs in the C code are the internal units described here. Converting to and from internal units is done in the Python part of galpy and there is currently no way to do this in the C version.

Basically, you'd have to:

• Choose values for the unit conversion parameters ro and vo, which are in kpc and km/s (standard ro=8., vo=220.).
• Divide input velocities by ro to remove the kpc or pc if you want
• Convert the output to force units using an implementation of galpy.util.bovy_conversion.force_in_pcMyr2, which as you can see is a pretty simple function.

pot_args is also in internal units; I think the only consistent and long-termstable way to generate pot_args is through the Python interface.

2. The potential calculation
   In integrateFullOrbit.h, the force calculation is available, how about the potential value for a given position?

You can evaluate the potential with evaluatePotentials in galpy_potentials.h (the forces are also in that file, not in integrateFullOrbit.h).

3. The list of pot_type and pot_args
   Right now I can manage to generate the pot_type and pot_args in the c++ interface. But without checking the python commander, is it possible to get the help information that show the corresponding type and definition of arguments for a potential? For exmple, some help information in a c function or a text file like:
   Name | Type | argument
   Miyamoto | 5 | A: def. B: def. ...
   ...
   For the PeTar users who want to directly set the type and argments of potentials without checking python galpy, this help information (will be printed in the petar commander with option "-h") can be very convenient. Right now I can get the type from the integrateFullOrbit.c, but the definitions of arguments are not shown there.

The potentials are not documented in C; I would either direct people to the Python documentation (available online) or write an automated piece of code to transform the Python docstrings to C docstrings (because transferring them by hand will make it hard to keep them consistent). As I said above, I think the only good way for people to get pot_args is through the Python interface.

4. The library file
   From the pip installed galpy, I did not find the libgalpy.so. But I manage to compile the library myself by including all .c files in galpy directory and it seems working correctly. I think pip probably only provide the source but not the compiled library, since the latter can depends on the machine configuration. I also prefer to compile the lib manually, since in some supercomputer, the pip may not work and the galpy source code has to be downloaded and installed. Also there is no root permission to install the lib in /user/local.

I think pip puts the libgalpy.so file in the ../site-packages/... somewhere, but for now, I think the easiest way to get it is to download/clone the galpy source and run python setup.py build, which creates the libgalpy.so file in the current directory (note that it has a Python version specific extension). I'm not sure how one gets pip to install in something like /usr/local.

5. OpenMP
   When I compile the galpy library, there is warning of OpenMP, how the multi threads are actually used in galpy? My idea is to add openmp parallel loop outside force calculation function, like
   #pragma omp parallel for
for (int i=0; i<n_particle; i++) {
calcAccPot(acc[i], pot[i], time, pos[i]);
}
   Where calcAccPot is the function to use calcRforce and calczforce to obtain the acceleration. In such case, will be some confliction existing for using the multi threads?

Currently, OpenMP is only used to parallelize orbit integrations or the evaluation of action-angle coordinates for many orbits. So none of the code that you are using to evaluate forces uses OpenMP internally, so you should be good wrapping it in OpenMP pragmas.

[lwang-astro]

Thank you for the explanation! I have tried the _parse_pot (python3.6 pip installed galpy) and iterate all potentials listed in galpy.potential.
Firstly, I generate an instance for each type of a potential by:
pot_list=[p for p in dir(galpy.potential) if 'Potential' in p]
for pot_name in pot_list:
pot_module = galpy.potential.__getattribute__(pot_name)
pot_instance = pot_module()

This loop seems includes several function name that are not potential class. If I removed all functions, I still got several errors:

<class 'galpy.potential.GaussianAmplitudeWrapperPotential.GaussianAmplitudeWrapperPotential'> WrapperPotentials are only supported in 3D and 2D
<class 'galpy.potential.SnapshotRZPotential.InterpSnapshotRZPotential'> init() missing 1 required positional argument: 's'
<class 'galpy.potential.MovingObjectPotential.MovingObjectPotential'> init() missing 1 required positional argument: 'orbit'
<class 'galpy.potential.NumericalPotentialDerivativesMixin.NumericalPotentialDerivativesMixin'> init() missing 1 required positional argument: 'kwargs'
<class 'galpy.potential.Potential.PotentialError'> init() missing 1 required positional argument: 'value'
<class 'galpy.potential.SnapshotRZPotential.SnapshotRZPotential'> init() missing 1 required positional argument: 's'
<class 'galpy.potential.SolidBodyRotationWrapperPotential.SolidBodyRotationWrapperPotential'> WrapperPotentials are only supported in 3D and 2D
<class 'galpy.potential.interpRZPotential.interpRZPotential'> 'NoneType' object has no attribute '_ro'

It seems some potentials cannot be initialized with no arguments, or defaulted arguments are missing.
The error of "kwargs" in NumericalPotentialDerivativesMixin seems to indicate a bug that "**" is missed before kwargs.

Besides, some potential modules seems not import in init.py, thus cannot be used to initialize instances.
verticalPotential
WrapperPotential
EllipsoidalPotential

Some potentials give empty pot_type and pot_args from _parse_pot, is this means that these potentials are not implemented in c code?
CosmphiDiskPotential
EllipticalDiskPotential
FerrersPotential
HenonHeilesPotential
IsothermalDiskPotential
KGPotential
LopsidedDiskPotential
RazorThinExponentialDiskPotential
RingPotential
SphericalShellPotential
SteadyLogSpiralPotential
TransientLogSpiralPotential
TwoPowerSphericalPotential
TwoPowerTriaxialPotential
linearPotential
planarAxiPotential
planarPotential

Right now all these listed potentials I probably cannot used in the c interface. Besides, to generate a help function to show all pot_type and pot_args, It is necessary that the potential classes can support no arguments initalization.

On the other hand, some potentials have very long list of arguments obtained from _parse_pot. Is it possible to initialize this potentials without manually inputting all arguments (using the default values instead) in the c interface?

[jobovy]

Thank you for the explanation! I have tried the _parse_pot (python3.6 pip installed galpy) and iterate all potentials listed in galpy.potential.
Firstly, I generate an instance for each type of a potential by:
pot_list=[p for p in dir(galpy.potential) if 'Potential' in p]
for pot_name in pot_list:
pot_module = galpy.potential.__getattribute__(pot_name)
pot_instance = pot_module()

This loop seems includes several function name that are not potential class. If I removed all functions, I still got several errors:

<class 'galpy.potential.GaussianAmplitudeWrapperPotential.GaussianAmplitudeWrapperPotential'> WrapperPotentials are only supported in 3D and 2D
<class 'galpy.potential.SnapshotRZPotential.InterpSnapshotRZPotential'> init() missing 1 required positional argument: 's'
<class 'galpy.potential.MovingObjectPotential.MovingObjectPotential'> init() missing 1 required positional argument: 'orbit'
<class 'galpy.potential.NumericalPotentialDerivativesMixin.NumericalPotentialDerivativesMixin'> init() missing 1 required positional argument: 'kwargs'
<class 'galpy.potential.Potential.PotentialError'> init() missing 1 required positional argument: 'value'
<class 'galpy.potential.SnapshotRZPotential.SnapshotRZPotential'> init() missing 1 required positional argument: 's'
<class 'galpy.potential.SolidBodyRotationWrapperPotential.SolidBodyRotationWrapperPotential'> WrapperPotentials are only supported in 3D and 2D
<class 'galpy.potential.interpRZPotential.interpRZPotential'> 'NoneType' object has no attribute '_ro'

It seems some potentials cannot be initialized with no arguments, or defaulted arguments are missing.
The error of "kwargs" in NumericalPotentialDerivativesMixin seems to indicate a bug that "**" is missed before kwargs.

Some of these are objects that cannot be used on their own (wrappers need to wrap a potential), and some just don't have default arguments, because there isn't a useful default (interpolated potentials need to interpolate something, a moving object potential needs a moving object to initialize it). PotentialError is an exception class, and NumericalPotentialDerivativesMixin is a Mixin for use in Python only currently (to add numerical derivatives, there's no bug in it).

Besides, some potential modules seems not import in init.py, thus cannot be used to initialize instances.
verticalPotential
WrapperPotential
EllipsoidalPotential

Yes, these are not standard potentials, but just help with implementing other potentials.

Some potentials give empty pot_type and pot_args from _parse_pot, is this means that these potentials are not implemented in c code?
CosmphiDiskPotential
EllipticalDiskPotential
FerrersPotential
HenonHeilesPotential
IsothermalDiskPotential
KGPotential
LopsidedDiskPotential
RazorThinExponentialDiskPotential
RingPotential
SphericalShellPotential
SteadyLogSpiralPotential
TransientLogSpiralPotential
TwoPowerSphericalPotential
TwoPowerTriaxialPotential
linearPotential
planarAxiPotential
planarPotential

These are either again not regular potentials (linearPotential and such) or they are potentials without a C implementation (like SphericalShellPotential) or potentials that are only defined in 2D or 1D. So I don't think you want to support any of these.

There is a function _check_c in galpy.potential.Potential that you can use to check whether a potential has a C implementation. There's also _dim there to check whether it's a 2D potential.

Right now all these listed potentials I probably cannot used in the c interface. Besides, to generate a help function to show all pot_type and pot_args, It is necessary that the potential classes can support no arguments initalization.

On the other hand, some potentials have very long list of arguments obtained from _parse_pot. Is it possible to initialize this potentials without manually inputting all arguments (using the default values instead) in the c interface?

These are the potentials using interpolation or the triaxial ones that have to do a numerical integral using Gaussian quadrature. There's no way around this.

[jobovy]

I should perhaps add that I would imagine users of your code who want to use an external force from galpy to do the following:

1. Run galpy in Python to obtain the pot_type and pot_args parameters
2. Input these into your code

I get the sense that you would want to avoid 1), but I wouldn't be supportive of that, because I think users should know they are using galpy.

[lwang-astro]

Thank your for the explanation. Now I understand how to do next.

My propose is not exactly escaping step 1. I prefer to write a short python script as a PeTar tool that makes the step 1 much simpler. This is why I have the questions above. The script can print the pot_type and pot_args when users give a potential name and list all available potentials that can be used in the c interface. This can be done by one commander instead of opening python, importing libs and using functions to get the result step by step.

Besides, the input style in the galpy c interface for PeTar is not exactly like the Python one. Thus the script help users to transfer output of _parse_pot to the Petar input style. The style of PeTar input is based on the GNU getopt (options with the format of -[char] [arguments] or --[string] [arguments]), thus I prefer to make the input of the pot_type and pot_args be consistent with the PeTar input style.

As you can see, if a new user of Galpy try to touch the c interface, it is not straight forwards for them to get the information about _parse_pot. Even they finally successfully to get the idea to use _parse_pot, they will suffer the same problems like me and have questions about why some potentials are not working. Unless they check the Galpy documents in details or ask you directly, it is not easy to finger out what to do next. Thus, such kind of script can help users to avoid these techincal details.

[lwang-astro]

I have implemented the help script and it seems working well and I will continue more tests.

I have a few more questions:

1. There seems no _check_c function available in galpy.potential.Potential. I got the error:
   AttributeError: type object 'Potential' has no attribute '_check_c'
   Anyway, the _dim is available. Right now without _check_c, I check whether pot_type returned from _parse_pot is empty to identify whether the potential is supported by c. Not sure whether this is the correct way.

These are the potentials using interpolation or the triaxial ones that have to do a numerical integral using Gaussian quadrature. There's no way around this.

Do you mean these arguments can evolve in a simulation? From the _parse_pot, I can get a group of arguments for these type of potentials,
For example, for DiskSCFPotential, _parse_pot return:
(2,
array([24, 26], dtype=int32),
array([ 1.00000000e+00, 0.00000000e+00, 1.00000000e+01, 1.00000000e+01,
1.00000000e+00, 4.41180299e-01, 1.74634964e-18, 7.92714987e-01,
-1.12454613e-17, -1.36997135e+00, -2.19276441e-17, 6.15085603e+00,
8.26255685e-17, -3.93641232e+01, -3.94730929e-15, -1.00231468e-01,
8.16850680e-18, -2.37057429e-01, 6.24591799e-18, 3.66521735e-01,
-3.29776984e-17, -1.43977729e+00, 2.80369022e-16, 8.05804688e+00,
-3.35018932e-15, -1.05782007e-02, 7.37000516e-19, 1.68127165e-02,
3.64035330e-18, -2.98658913e-02, 9.39285386e-18, 1.15327085e-01,
-2.51734254e-17, -6.19559823e-01, -1.56253329e-17, 1.79125695e-02,
1.04838826e-18, 4.97267108e-03, 7.99638175e-19, 2.87868690e-03,
-1.20222817e-18, -2.45394591e-02, -7.57611599e-18, 1.60212615e-01,
3.33701572e-17, -6.72983557e-04, -4.25813745e-19, 9.70344791e-04,
2.45361654e-19, -2.10559873e-03, -8.93306929e-19, 7.30872323e-03,
-4.83163661e-19, -3.28932107e-02, 1.79756083e-17, -3.19236841e-03,
2.53258627e-20, -5.88086527e-04, 4.54835504e-20, 3.55540841e-04,
2.05116852e-21, -1.19120763e-03, 1.52641780e-18, 6.12148111e-03,
3.33872426e-18, -1.55389962e-04, 7.92624032e-21, -1.22492105e-04,
-5.52187420e-21, -4.59613396e-05, -7.57094313e-21, 3.72277075e-04,
-8.34602347e-20, -1.85288454e-03, -1.96989490e-18, 6.09188568e-04,
2.20950542e-20, -2.27545926e-06, 2.52728464e-20, 3.84412226e-05,
3.55311844e-20, -1.15246697e-04, 2.13226520e-19, 4.36217106e-04,
-2.09593779e-18, 2.12732704e-04, 1.60302830e-20, 2.48088780e-05,
3.47347292e-20, -6.35536263e-06, 3.59536675e-20, 2.64502059e-05,
-1.88636964e-19, -1.21508050e-04, -7.02771155e-20, -8.75766524e-05,
-2.82406227e-20, 6.65527544e-06, -1.47363792e-20, 9.29565774e-07,
8.01042101e-21, -1.05411611e-05, -5.52153961e-20, 3.94887507e-05,
2.21552941e-19, -1.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
3.00000000e+00, 0.00000000e+00, 1.25663706e+01, 3.33333333e-01,
0.00000000e+00, 3.70370370e-02]))

what are the definitaion of these values? Like the values at time zero?

3. For the DiskSCFPotential, when I create an instance
   pot=galpy.potential.DiskSCFPotential()
   I got the warning:
   DiskSCFPotential.py:478: RuntimeWarning: divide by zero encountered in double_scalars
out-= a*(s(r)*h(z)+d2s(r)*H(z)+2./r*ds(r)*(H(z)+z*dH(z)))
DiskSCFPotential.py:478: RuntimeWarning: invalid value encountered in double_scalars
out-= a*(s(r)*h(z)+d2s(r)*H(z)+2./r*ds(r)*(H(z)+z*dH(z)))
   Is this the correct way to create the instance for this potential? I did not get a similar warning for other potentials.

[jobovy]

Okay, so it sounds like you would still have your users install galpy and your script is a wrapper around it to get the arguments in the correct form (and provide help)? That seems fine to me, as long as it's clear to people that they are using galpy in the back. Another option would be add functions to galpy to produce the input to your code, like the to_nemo... functions and to_amuse functions available now.

I would say that I don't think most users should worry about the C interface, so I don't consider the difficulty in using it a big deal. I don't imagine that people using your code with galpy would have to really know anything about the galpy C interface and that just reading the galpy documentation should not be too much of a burden for them (many of them will probably already be familiar with galpy).

I'm able to access _check_c, e.g.,

from galpy.potential.Potential import _check_c                          
from galpy.potential import MWPotential2014                             
_check_c(MWPotential2014)  
# True

Does that not work for you?

The long parameter lists don't evolve, they are just static numbers that define the potential

(MovingObjectPotential is a special case, because it contains the orbit of a moving object whose force is calculated; for that you have to make sure that whatever time you run the simulation for overlaps with the time for which the orbit given to MovingObjectPotential is integrated; I imagine that your users probably wouldn't want to use MovingObjectPotential).

For the SCFPotential and DiskSCFPotential the long parameter list is mainly the expansion coefficients in the SCF expansion of the potential. I think that warning you get with the default DiskSCFPotential is fine.

If it were me, I would probably read the parameters of the potential from a file that contains them and then your interface to galpy would just have to write that file. That way you don't have to put ridiculously long argument lists when using potentials with many parameters.

[lwang-astro]

In petar -h and running time, users will see that they are using the galpy library. Also, the corresponding reference paper is shown. To use Galpy, users also need to install it first. So I think they will know that Galpy is used.

The help script is a small tool, when users run it in terminal, it printed the pot_type, pot_arg and the __doc__ of a given potential name, or list all supported potentials. It also provide a function to generate petar style input or a configure file from an instance of potential. If you would like to check details, it is shown in https://github.com/lwang-astro/PeTar/tree/master/galpy-interface.

I think many potential users of PeTar are interesting in the internal stellar dynamics of star clusters while they require an external potential to drive the tidal dissolution of the clusters. In most cases, they will probably only use the default Milkway potential or import existing potentials with default arguments. This help script is mostly useful for them to quickly set up the potential. Of course for users who want to design speficial potential models, they probably need to use the Galpy python interface to generate the instance first and then get type and argments.

For the potential with a long argument list, I am also thinking to add an additional configure file to read and use help script to generate the configure file. Thanks for suggestion.

For the _check_c, I think I have missed the
from galpy.potential.Potential import _check_c
When I add this, it works now.

In the filtered potential list from my help script, the MovingObjectPotential seems not in the list, so I think it is fine.

[jobovy]

Great, sounds like everything is working fine now then!

For future reference, I think this is a list of changes to make to galpy to make the interface easier:

• Install the shared-object library somewhere easily accessible as part of the installation (both python setup.py and pip install)
• Add all external declarations to a single header file (probably add all to galpy_potentials.h) and install it as galpy.h.
• Install the header file(s) somewhere easily accessible as well
• Make _parse_pot a public function somehow and document it better
• Make _check_c a public function to easily check whether a Potential or combination of potentials has a C implementation.
• Write some development documentation explaining the C code better

I think that's all I took away from this thread, but please add if there is anything I missed.

[lwang-astro]

Today I have tried to simulate a star cluster by using MWPotential2014. The result looks reasonable. I will do more test and compare with the result from Galpy python. Thanks for all the help!

I think your list of future reference is good. I have one more point to add. For the c interface, right now I include orbit/orbit_c_ext/integrateFullOrbit.h to get the parse_leapFuncArgs_Full function. I think you have mentioned that this declaration can be moved to galpy_potentials.h. Then, if only the force/potential calculation functions are needed in the interface, probably only the galpy_potentials.h is necessary to be included. This may simplify the interface.

[jobovy]

Great to hear that and thanks for the additional future point, I've added it to the list above.