Skip to content

Commit

Permalink
Also assign potential evaluation in integratePlanarOrbit; always prop…
Browse files Browse the repository at this point in the history 
…agate q for log halo such that we can use the same potentialEval function for planar and full
  • Loading branch information
@jobovy
jobovy committed Jan 18, 2018
1 parent f454db4 commit 5efef89
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 15 deletions.
5 changes: 3 additions & 2 deletions galpy/orbit_src/integratePlanarOrbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,9 +66,10 @@ def _parse_pot(pot):
and isinstance(p._Pot,potential.LogarithmicHaloPotential):
pot_type.append(0)
if p._Pot.isNonAxi:
pot_args.extend([p._Pot._amp,p._Pot._core2,p._Pot._1m1overb2])
pot_args.extend([p._Pot._amp,p._Pot._q,
p._Pot._core2,p._Pot._1m1overb2])
else:
pot_args.extend([p._Pot._amp,p._Pot._core2,2.]) # 1m1overb2 > 1: axi
pot_args.extend([p._Pot._amp,p._Pot._q,p._Pot._core2,2.]) # 1m1overb2 > 1: axi
elif isinstance(p,potential_src.planarPotential.planarPotentialFromFullPotential) \
and isinstance(p._Pot,potential.DehnenBarPotential):
pot_type.append(1)
Expand Down
26 changes: 24 additions & 2 deletions galpy/orbit_src/orbit_c_ext/integratePlanarOrbit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,14 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
init_potentialArgs(npot,potentialArgs);
for (ii=0; ii < npot; ii++){
switch ( *(*pot_type)++ ) {
case 0: //LogarithmicHaloPotential, 3 arguments
case 0: //LogarithmicHaloPotential, 4 arguments
potentialArgs->potentialEval= &LogarithmicHaloPotentialEval;
potentialArgs->planarRforce= &LogarithmicHaloPotentialPlanarRforce;
potentialArgs->planarphiforce= &LogarithmicHaloPotentialPlanarphiforce;
potentialArgs->planarR2deriv= &LogarithmicHaloPotentialPlanarR2deriv;
potentialArgs->planarphi2deriv= &LogarithmicHaloPotentialPlanarphi2deriv;
potentialArgs->planarRphideriv= &LogarithmicHaloPotentialPlanarRphideriv;
potentialArgs->nargs= 3;
potentialArgs->nargs= 4;
break;
case 1: //DehnenBarPotential, 6 arguments
potentialArgs->planarRforce= &DehnenBarPotentialPlanarRforce;
Expand Down Expand Up @@ -66,6 +67,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 6;
break;
case 5: //MiyamotoNagaiPotential, 3 arguments
potentialArgs->potentialEval= &MiyamotoNagaiPotentialEval;
potentialArgs->planarRforce= &MiyamotoNagaiPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &MiyamotoNagaiPotentialPlanarR2deriv;
Expand All @@ -82,6 +84,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 4;
break;
case 7: //PowerSphericalPotential, 2 arguments
potentialArgs->potentialEval= &PowerSphericalPotentialEval;
potentialArgs->planarRforce= &PowerSphericalPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &PowerSphericalPotentialPlanarR2deriv;
Expand All @@ -90,6 +93,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 8: //HernquistPotential, 2 arguments
potentialArgs->potentialEval= &HernquistPotentialEval;
potentialArgs->planarRforce= &HernquistPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &HernquistPotentialPlanarR2deriv;
Expand All @@ -98,6 +102,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 9: //NFWPotential, 2 arguments
potentialArgs->potentialEval= &NFWPotentialEval;
potentialArgs->planarRforce= &NFWPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &NFWPotentialPlanarR2deriv;
Expand All @@ -106,6 +111,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 10: //JaffePotential, 2 arguments
potentialArgs->potentialEval= &JaffePotentialEval;
potentialArgs->planarRforce= &JaffePotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &JaffePotentialPlanarR2deriv;
Expand All @@ -114,6 +120,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 11: //DoubleExponentialDiskPotential, XX arguments
potentialArgs->potentialEval= &DoubleExponentialDiskPotentialEval;
potentialArgs->planarRforce= &DoubleExponentialDiskPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
//potentialArgs->planarR2deriv= &DoubleExponentialDiskPotentialPlanarR2deriv;
Expand All @@ -123,6 +130,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= (int) (8 + 2 * *(*pot_args+5) + 4 * ( *(*pot_args+4) + 1 ));
break;
case 12: //FlattenedPowerPotential, 4 arguments
potentialArgs->potentialEval= &FlattenedPowerPotentialEval;
potentialArgs->planarRforce= &FlattenedPowerPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &FlattenedPowerPotentialPlanarR2deriv;
Expand All @@ -131,6 +139,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 3;
break;
case 14: //IsochronePotential, 2 arguments
potentialArgs->potentialEval= &IsochronePotentialEval;
potentialArgs->planarRforce= &IsochronePotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &IsochronePotentialPlanarR2deriv;
Expand All @@ -139,6 +148,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 15: //PowerSphericalPotentialwCutoff, 3 arguments
potentialArgs->potentialEval= &PowerSphericalPotentialwCutoffEval;
potentialArgs->planarRforce= &PowerSphericalPotentialwCutoffPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &PowerSphericalPotentialwCutoffPlanarR2deriv;
Expand All @@ -147,6 +157,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 3;
break;
case 16: //KuzminKutuzovStaeckelPotential, 3 arguments
potentialArgs->potentialEval= &KuzminKutuzovStaeckelPotentialEval;
potentialArgs->planarRforce= &KuzminKutuzovStaeckelPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &KuzminKutuzovStaeckelPotentialPlanarR2deriv;
Expand All @@ -155,6 +166,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 3;
break;
case 17: //PlummerPotential, 2 arguments
potentialArgs->potentialEval= &PlummerPotentialEval;
potentialArgs->planarRforce= &PlummerPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &PlummerPotentialPlanarR2deriv;
Expand All @@ -163,6 +175,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 18: //PseudoIsothermalPotential, 2 arguments
potentialArgs->potentialEval= &PseudoIsothermalPotentialEval;
potentialArgs->planarRforce= &PseudoIsothermalPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &PseudoIsothermalPotentialPlanarR2deriv;
Expand All @@ -171,6 +184,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 19: //KuzminDiskPotential, 2 arguments
potentialArgs->potentialEval= &KuzminDiskPotentialEval;
potentialArgs->planarRforce= &KuzminDiskPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &KuzminDiskPotentialPlanarR2deriv;
Expand All @@ -179,6 +193,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 20: //BurkertPotential, 2 arguments
potentialArgs->potentialEval= &BurkertPotentialEval;
potentialArgs->planarRforce= &BurkertPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->planarR2deriv= &BurkertPotentialPlanarR2deriv;
Expand All @@ -187,21 +202,25 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= 2;
break;
case 21: //TriaxialHernquistPotential, lots of arguments
potentialArgs->potentialEval= &TriaxialHernquistPotentialEval;
potentialArgs->planarRforce= &TriaxialHernquistPotentialPlanarRforce;
potentialArgs->planarphiforce= &TriaxialHernquistPotentialPlanarphiforce;
potentialArgs->nargs= (int) (21 + 2 * *(*pot_args+14));
break;
case 22: //TriaxialNFWPotential, lots of arguments
potentialArgs->potentialEval= &TriaxialNFWPotentialEval;
potentialArgs->planarRforce= &TriaxialNFWPotentialPlanarRforce;
potentialArgs->planarphiforce= &TriaxialNFWPotentialPlanarphiforce;
potentialArgs->nargs= (int) (21 + 2 * *(*pot_args+14));
break;
case 23: //TriaxialJaffePotential, lots of arguments
potentialArgs->potentialEval= &TriaxialJaffePotentialEval;
potentialArgs->planarRforce= &TriaxialJaffePotentialPlanarRforce;
potentialArgs->planarphiforce= &TriaxialJaffePotentialPlanarphiforce;
potentialArgs->nargs= (int) (21 + 2 * *(*pot_args+14));
break;
case 24: //SCFPotential, many arguments
potentialArgs->potentialEval= &SCFPotentialEval;
potentialArgs->planarRforce= &SCFPotentialPlanarRforce;
potentialArgs->planarphiforce= &SCFPotentialPlanarphiforce;
potentialArgs->planarR2deriv= &SCFPotentialPlanarR2deriv;
Expand All @@ -210,11 +229,13 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
potentialArgs->nargs= (int) (5 + (1 + *(*pot_args + 1)) * *(*pot_args+2) * *(*pot_args+3)* *(*pot_args+4) + 7);
break;
case 25: //SoftenedNeedleBarPotential, 13 arguments
potentialArgs->potentialEval= &SoftenedNeedleBarPotentialEval;
potentialArgs->planarRforce= &SoftenedNeedleBarPotentialPlanarRforce;
potentialArgs->planarphiforce= &SoftenedNeedleBarPotentialPlanarphiforce;
potentialArgs->nargs= (int) 13;
break;
case 26: //DiskSCFPotential, nsigma+3 arguments
potentialArgs->potentialEval= &DiskSCFPotentialEval;
potentialArgs->planarRforce= &DiskSCFPotentialPlanarRforce;
potentialArgs->planarphiforce= &ZeroPlanarForce;
potentialArgs->nargs= (int) **pot_args + 3;
Expand Down Expand Up @@ -245,6 +266,7 @@ void parse_leapFuncArgs(int npot,struct potentialArg * potentialArgs,
break;
//////////////////////////////// WRAPPERS /////////////////////////////////////
case -1: //DehnenSmoothWrapperPotential
potentialArgs->potentialEval= &DehnenSmoothWrapperPotentialEval;
potentialArgs->planarRforce= &DehnenSmoothWrapperPotentialPlanarRforce;
potentialArgs->planarphiforce= &DehnenSmoothWrapperPotentialPlanarphiforce;
potentialArgs->planarR2deriv= &DehnenSmoothWrapperPotentialPlanarR2deriv;
Expand Down
22 changes: 11 additions & 11 deletions galpy/potential_src/potential_c_ext/LogarithmicHaloPotential.c
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#include <math.h>
#include <galpy_potentials.h>
//LogarithmicHaloPotential
//3 (2) arguments: amp, c2, (and q), now also 1-1/b^2 for triaxial!
//4 (3) arguments: amp, c2, (and q), now also 1-1/b^2 for triaxial!
double LogarithmicHaloPotentialEval(double R,double Z, double phi,
double t,
struct potentialArg * potentialArgs){
Expand Down Expand Up @@ -42,8 +42,8 @@ double LogarithmicHaloPotentialPlanarRforce(double R,double phi,
double * args= potentialArgs->args;
//Get args
double amp= *args;
double c= *(args+1);
double onem1overb2= *(args+2);
double c= *(args+2); // skip q
double onem1overb2= *(args+3);
//Calculate Rforce
double Rt2;
if ( onem1overb2 < 1 ) {
Expand Down Expand Up @@ -95,8 +95,8 @@ double LogarithmicHaloPotentialPlanarphiforce(double R,double phi,
double * args= potentialArgs->args;
//Get args
double amp= *args;
double c= *(args+1);
double onem1overb2= *(args+2);
double c= *(args+2); // skip q
double onem1overb2= *(args+3);
//Calculate phiforce
double Rt2;
if ( onem1overb2 < 1 ) {
Expand All @@ -111,8 +111,8 @@ double LogarithmicHaloPotentialPlanarR2deriv(double R,double phi,
double * args= potentialArgs->args;
//Get args
double amp= *args;
double c= *(args+1);
double onem1overb2= *(args+2);
double c= *(args+2); // skip q
double onem1overb2= *(args+3);
//Calculate Rforce
double Rt2;
if ( onem1overb2 < 1 ) {
Expand All @@ -127,8 +127,8 @@ double LogarithmicHaloPotentialPlanarphi2deriv(double R,double phi,
double * args= potentialArgs->args;
//Get args
double amp= *args;
double c= *(args+1);
double onem1overb2= *(args+2);
double c= *(args+2); // skip q
double onem1overb2= *(args+3);
//Calculate Rforce
double Rt2;
if ( onem1overb2 < 1 ) {
Expand All @@ -144,8 +144,8 @@ double LogarithmicHaloPotentialPlanarRphideriv(double R,double phi,
double * args= potentialArgs->args;
//Get args
double amp= *args;
double c= *(args+1);
double onem1overb2= *(args+2);
double c= *(args+2); // skip q
double onem1overb2= *(args+3);
//Calculate Rforce
double Rt2;
if ( onem1overb2 < 1 ) {
Expand Down

0 comments on commit 5efef89

Please sign in to comment.