Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement caching of the mass for the PowerSphericalPotentialwCutoff potential #720

Merged
merged 1 commit into from
Mar 12, 2025
+31 −12
Merged

Implement caching of the mass for the PowerSphericalPotentialwCutoff potential #720

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions galpy/orbit/integrateFullOrbit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ def _parse_pot(pot, potforactions=False, potfortorus=False):
elif isinstance(p, potential.PowerSphericalPotentialwCutoff):
pot_type.append(15)
pot_args.extend([p._amp, p.alpha, p.rc])
pot_args.extend([0.0, 0.0]) # for caching
elif isinstance(p, potential.MN3ExponentialDiskPotential):
# Three Miyamoto-Nagai disks
npot += 2
Expand Down
4 changes: 2 additions & 2 deletions galpy/orbit/orbit_c_ext/integrateFullOrbit.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@ void parse_leapFuncArgs_Full(int npot,
potentialArgs->ntfuncs= 0;
potentialArgs->requiresVelocity= false;
break;
case 15: //PowerSphericalwCutoffPotential, 3 arguments
case 15: //PowerSphericalwCutoffPotential, 5 arguments
potentialArgs->potentialEval= &PowerSphericalPotentialwCutoffEval;
potentialArgs->Rforce= &PowerSphericalPotentialwCutoffRforce;
potentialArgs->zforce= &PowerSphericalPotentialwCutoffzforce;
Expand All @@ -243,7 +243,7 @@ void parse_leapFuncArgs_Full(int npot,
//potentialArgs->R2deriv= &PowerSphericalPotentialR2deriv;
//potentialArgs->planarphi2deriv= &ZeroForce;
//potentialArgs->planarRphideriv= &ZeroForce;
potentialArgs->nargs= 3;
potentialArgs->nargs= 5;
potentialArgs->ntfuncs= 0;
potentialArgs->requiresVelocity= false;
break;
Expand Down
38 changes: 28 additions & 10 deletions galpy/potential/potential_c_ext/PowerSphericalPotentialwCutoff.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
#define M_PI 3.14159265358979323846
#endif
//PowerSphericalPotentialwCutoff
//3 arguments: amp, alpha, rc
//5 arguments: amp, alpha, rc, cached_r2, cached_force
double mass(double r2,double alpha, double rc){
return 2. * M_PI * pow ( rc , 3. - alpha ) * ( gsl_sf_gamma ( 1.5 - 0.5 * alpha ) * gsl_sf_gamma_inc_P ( 1.5 - 0.5 * alpha , r2 / rc / rc ) );
}
Expand All @@ -27,13 +27,22 @@
struct potentialArg * potentialArgs){
double * args= potentialArgs->args;
//Get args
double amp= *args++;
double alpha= *args++;
double rc= *args;
double amp= *args;
double alpha= *(args + 1);
double rc= *(args + 2);
//Radius
double r2= R*R+Z*Z;
double force;
if ( r2 == *(args + 3) ) {
force= *(args + 4);
} else {
force= -amp * mass(r2, alpha, rc) / pow(r2,1.5);
//Setup caching
*(args + 3)= r2;
*(args + 4)= force;
}
//Calculate Rforce
return - amp * mass (r2,alpha,rc) * R / pow(r2,1.5);
return force * R;
}
double PowerSphericalPotentialwCutoffPlanarRforce(double R,double phi,
double t,
Expand All @@ -53,13 +62,22 @@
struct potentialArg * potentialArgs){
double * args= potentialArgs->args;
//Get args
double amp= *args++;
double alpha= *args++;
double rc= *args;
double amp= *args;
double alpha= *(args + 1);
double rc= *(args + 2);
//Radius
double r2= R*R+Z*Z;
//Calculate Rforce
return - amp * mass (r2,alpha,rc) * Z / pow(r2,1.5);
double force;
if ( r2 == *(args + 3) ) {
force= *(args + 4);
} else {
force= -amp * mass(r2, alpha, rc) / pow(r2,1.5);

Check warning on line 74 in galpy/potential/potential_c_ext/PowerSphericalPotentialwCutoff.c

View check run for this annotation

Codecov / codecov/patch

galpy/potential/potential_c_ext/PowerSphericalPotentialwCutoff.c#L74 

Added line #L74 was not covered by tests
//Setup caching
*(args + 3)= r2;
*(args + 4)= force;

Check warning on line 77 in galpy/potential/potential_c_ext/PowerSphericalPotentialwCutoff.c

View check run for this annotation

Codecov / codecov/patch

galpy/potential/potential_c_ext/PowerSphericalPotentialwCutoff.c#L76-L77 

Added lines #L76 - L77 were not covered by tests
}
//Calculate zforce
return force * Z;
}
double PowerSphericalPotentialwCutoffPlanarR2deriv(double R,double phi,
double t,
Expand Down
Loading