Skip to content

Commit

Permalink
Perform non-symplectic integration in C in cylindrical coordinates, s…
Browse files Browse the repository at this point in the history 
…olving Hamilton's equations
  • Loading branch information
@jobovy
jobovy committed Dec 24, 2018
1 parent ff3b82c commit 444aa14
Show file tree
Hide file tree
Showing 7 changed files with 49 additions and 77 deletions.
34 changes: 3 additions & 31 deletions galpy/orbit/FullOrbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -612,41 +612,13 @@ def metric(x,y):
scaling=scaling,metric=metric)
elif ext_loaded and \
(method.lower() == 'leapfrog_c' or method.lower() == 'symplec4_c' \
or method.lower() == 'symplec6_c'):
or method.lower() == 'symplec6_c' or method.lower() == 'rk4_c' \
or method.lower() == 'rk6_c' or method.lower() == 'dopr54_c' \
or method.lower() == 'dop853_c'):
warnings.warn("Using C implementation to integrate orbits",
galpyWarningVerbose)
#integrate
out, msg= integrateFullOrbit_c(pot,nu.copy(vxvv),t,method,dt=dt)
elif ext_loaded and \
(method.lower() == 'rk4_c' \
or method.lower() == 'rk6_c' \
or method.lower() == 'dopr54_c' \
or method.lower() == 'dop853_c'):
warnings.warn("Using C implementation to integrate orbits",
galpyWarningVerbose)
#go to the rectangular frame
this_vxvv= nu.array([vxvv[0]*nu.cos(vxvv[5]),
vxvv[0]*nu.sin(vxvv[5]),
vxvv[3],
vxvv[1]*nu.cos(vxvv[5])-vxvv[2]*nu.sin(vxvv[5]),
vxvv[2]*nu.cos(vxvv[5])+vxvv[1]*nu.sin(vxvv[5]),
vxvv[4]])
#integrate
tmp_out, msg= integrateFullOrbit_c(pot,this_vxvv,
t,method,dt=dt)
#go back to the cylindrical frame
R= nu.sqrt(tmp_out[:,0]**2.+tmp_out[:,1]**2.)
phi= nu.arccos(tmp_out[:,0]/R)
phi[(tmp_out[:,1] < 0.)]= 2.*nu.pi-phi[(tmp_out[:,1] < 0.)]
vR= tmp_out[:,3]*nu.cos(phi)+tmp_out[:,4]*nu.sin(phi)
vT= tmp_out[:,4]*nu.cos(phi)-tmp_out[:,3]*nu.sin(phi)
out= nu.zeros((len(t),6))
out[:,0]= R
out[:,1]= vR
out[:,2]= vT
out[:,5]= phi
out[:,3]= tmp_out[:,2]
out[:,4]= tmp_out[:,5]
elif method.lower() == 'odeint' or method.lower() == 'dop853' or not ext_loaded:
init= [vxvv[0],vxvv[1],vxvv[0]*vxvv[2],vxvv[3],vxvv[4],vxvv[5]]
if method == 'dop853':
Expand Down
13 changes: 3 additions & 10 deletions galpy/orbit/linearOrbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,20 +172,13 @@ def _integrateLinearOrbit(vxvv,pot,t,method,dt):
return dop853(func=_linearEOM, x=vxvv, t=t, args=(pot,))
elif ext_loaded and \
(method.lower() == 'leapfrog_c' or method.lower() == 'symplec4_c'
or method.lower() == 'symplec6_c'):
or method.lower() == 'symplec6_c' or method.lower() == 'rk4_c' \
or method.lower() == 'rk6_c' or method.lower() == 'dopr54_c' \
or method.lower() == 'dop853_c'):
warnings.warn("Using C implementation to integrate orbits",
galpyWarningVerbose)
out, msg= integrateLinearOrbit_c(pot,nu.copy(vxvv),t,method,dt=dt)
return out
elif ext_loaded and \
(method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \
or method.lower() == 'rk6_c' \
or method.lower() == 'dopr54_c' \
or method.lower() == 'dop853_c'):
warnings.warn("Using C implementation to integrate orbits",
galpyWarningVerbose)
out, msg= integrateLinearOrbit_c(pot,nu.array(vxvv),t,method,dt=dt)
return out
elif method.lower() == 'odeint' or not ext_loaded:
return integrate.odeint(_linearEOM,vxvv,t,args=(pot,),rtol=10.**-8.)

Expand Down
23 changes: 19 additions & 4 deletions galpy/orbit/orbit_c_ext/integrateFullOrbit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@ void evalRectForce(double, double *, double *,
int, struct potentialArg *);
void evalRectDeriv(double, double *, double *,
int, struct potentialArg *);
void evalCylDeriv(double, double *, double *,
int, struct potentialArg *);
void evalRectDeriv_dxdv(double,double *, double *,
int, struct potentialArg *);
void symplec_drift_3d(double, double *);
Expand Down Expand Up @@ -388,6 +390,7 @@ EXPORT void integrateFullOrbit(double *yo,
int * err,
int odeint_type){
//Set up the forces, first count
int ii;
int dim= 6;
struct potentialArg * potentialArgs= (struct potentialArg *) malloc ( npot * sizeof (struct potentialArg) );
parse_leapFuncArgs_Full(npot,potentialArgs,&pot_type,&pot_args);
Expand Down Expand Up @@ -442,26 +445,28 @@ EXPORT void integrateFullOrbit(double *yo,
double *,int *);
void (*odeint_deriv_func)(double, double *, double *,
int,struct potentialArg *);
odeint_deriv_func= &evalCylDeriv;
// vT --> Lz
*(yo+2)*= *yo;
switch ( odeint_type ) {
case 1: //RK4
odeint_func= &bovy_rk4;
odeint_deriv_func= &evalRectDeriv;
break;
case 2: //RK6
odeint_func= &bovy_rk6;
odeint_deriv_func= &evalRectDeriv;
break;
case 5: //DOPR54
odeint_func= &bovy_dopr54;
odeint_deriv_func= &evalRectDeriv;
break;
case 6: //DOP853
odeint_func= &dop853;
odeint_deriv_func= &evalRectDeriv;
break;
}
odeint_func(odeint_deriv_func,dim,yo,nt,dt,t,npot,potentialArgs,rtol,atol,
result,err);
// Lz --> vT
for (ii=0; ii < nt; ii++)
*(result+2+ii*6)/= *(result+ii*6);
}
//Free allocated memory
free_potentialArgs(npot,potentialArgs);
Expand Down Expand Up @@ -599,6 +604,16 @@ void evalRectDeriv(double t, double *q, double *a,
*a++= sinphi*Rforce+1./R*cosphi*phiforce;
*a= zforce;
}
void evalCylDeriv(double t, double *y, double *a,
int nargs, struct potentialArg * potentialArgs){
*a = *(y+1);
*(a+1)= *(y+2) * *(y+2) / *y / *y / *y + \
calcRforce(*y,*(y+3),*(y+5),t,nargs,potentialArgs);
*(a+2)= calcPhiforce(*y,*(y+3),*(y+5),t,nargs,potentialArgs);
*(a+3)= *(y+4);
*(a+4)= calczforce(*y,*(y+3),*(y+5),t,nargs,potentialArgs);
*(a+5)= *(y+2) / *y / *y;
}

// LCOV_EXCL_START
void evalRectDeriv_dxdv(double t, double *q, double *a,
Expand Down
21 changes: 17 additions & 4 deletions galpy/orbit/orbit_c_ext/integratePlanarOrbit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@ void evalPlanarRectForce(double, double *, double *,
int, struct potentialArg *);
void evalPlanarRectDeriv(double, double *, double *,
int, struct potentialArg *);
void evalPlanarCylDeriv(double, double *, double *,
int, struct potentialArg *);
void evalPlanarRectDeriv_dxdv(double, double *, double *,
int, struct potentialArg *);
void symplec_drift_2d(double, double *);
Expand Down Expand Up @@ -380,6 +382,7 @@ EXPORT void integratePlanarOrbit(double *yo,
int * err,
int odeint_type){
//Set up the forces, first count
int ii;
int dim= 4;
struct potentialArg * potentialArgs= (struct potentialArg *) malloc ( npot * sizeof (struct potentialArg) );
parse_leapFuncArgs(npot,potentialArgs,&pot_type,&pot_args);
Expand Down Expand Up @@ -434,26 +437,28 @@ EXPORT void integratePlanarOrbit(double *yo,
double *,int *);
void (*odeint_deriv_func)(double, double *, double *,
int,struct potentialArg *);
odeint_deriv_func= &evalPlanarCylDeriv;
// vT --> Lz
*(yo+2)*= *yo;
switch ( odeint_type ) {
case 1: //RK4
odeint_func= &bovy_rk4;
odeint_deriv_func= &evalPlanarRectDeriv;
break;
case 2: //RK6
odeint_func= &bovy_rk6;
odeint_deriv_func= &evalPlanarRectDeriv;
break;
case 5: //DOPR54
odeint_func= &bovy_dopr54;
odeint_deriv_func= &evalPlanarRectDeriv;
break;
case 6: //DOP853
odeint_func= &dop853;
odeint_deriv_func= &evalPlanarRectDeriv;
break;
}
odeint_func(odeint_deriv_func,dim,yo,nt,dt,t,npot,potentialArgs,rtol,atol,
result,err);
// Lz --> vT
for (ii=0; ii < nt; ii++)
*(result+2+ii*4)/= *(result+ii*4);
}
//Free allocated memory
free_potentialArgs(npot,potentialArgs);
Expand Down Expand Up @@ -580,6 +585,14 @@ void evalPlanarRectDeriv(double t, double *q, double *a,
*a++= cosphi*Rforce-1./R*sinphi*phiforce;
*a= sinphi*Rforce+1./R*cosphi*phiforce;
}
void evalPlanarCylDeriv(double t, double *y, double *a,
int nargs, struct potentialArg * potentialArgs){
*a = *(y+1);
*(a+1)= *(y+2) * *(y+2) / *y / *y / *y + \
calcPlanarRforce(*y,*(y+3),t,nargs,potentialArgs);
*(a+2)= calcPlanarphiforce(*y,*(y+3),t,nargs,potentialArgs);
*(a+3)= *(y+2) / *y / *y;
}

void evalPlanarRectDeriv_dxdv(double t, double *q, double *a,
int nargs, struct potentialArg * potentialArgs){
Expand Down
29 changes: 3 additions & 26 deletions galpy/orbit/planarOrbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -625,35 +625,12 @@ def metric(x,y):
msg= 0
elif ext_loaded and \
(method.lower() == 'leapfrog_c' or method.lower() == 'symplec4_c' \
or method.lower() == 'symplec6_c'):
or method.lower() == 'symplec6_c' or method.lower() == 'rk4_c' \
or method.lower() == 'rk6_c' or method.lower() == 'dopr54_c' \
or method.lower() == 'dop853_c'):
warnings.warn("Using C implementation to integrate orbits",galpyWarningVerbose)
#integrate
out, msg= integratePlanarOrbit_c(pot,nu.copy(vxvv),t,method,dt=dt)
elif ext_loaded and \
(method.lower() == 'leapfrog_c' or method.lower() == 'rk4_c' \
or method.lower() == 'rk6_c' \
or method.lower() == 'dopr54_c' \
or method.lower() == 'dop853_c'):
warnings.warn("Using C implementation to integrate orbits",galpyWarningVerbose)
#go to the rectangular frame
this_vxvv= nu.array([vxvv[0]*nu.cos(vxvv[3]),
vxvv[0]*nu.sin(vxvv[3]),
vxvv[1]*nu.cos(vxvv[3])-vxvv[2]*nu.sin(vxvv[3]),
vxvv[2]*nu.cos(vxvv[3])+vxvv[1]*nu.sin(vxvv[3])])
#integrate
tmp_out, msg= integratePlanarOrbit_c(pot,this_vxvv,
t,method,dt=dt)
#go back to the cylindrical frame
R= nu.sqrt(tmp_out[:,0]**2.+tmp_out[:,1]**2.)
phi= nu.arccos(tmp_out[:,0]/R)
phi[(tmp_out[:,1] < 0.)]= 2.*nu.pi-phi[(tmp_out[:,1] < 0.)]
vR= tmp_out[:,2]*nu.cos(phi)+tmp_out[:,3]*nu.sin(phi)
vT= tmp_out[:,3]*nu.cos(phi)-tmp_out[:,2]*nu.sin(phi)
out= nu.zeros((len(t),4))
out[:,0]= R
out[:,1]= vR
out[:,2]= vT
out[:,3]= phi
elif method.lower() == 'odeint' or method.lower() == 'dop853' or not ext_loaded:
init= [vxvv[0],vxvv[1],vxvv[0]*vxvv[2],vxvv[3]]
if method == 'dop853':
Expand Down
4 changes: 2 additions & 2 deletions galpy/util/bovy_rk.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
C implementations of Runge-Kutta integrators
*/
/*
Copyright (c) 2011, Jo Bovy
Copyright (c) 2011, 2018 Jo Bovy
All rights reserved.
Redistribution and use in source and binary forms, with or without
Expand Down Expand Up @@ -60,7 +60,7 @@ Runge-Kutta 4 integrator
double *args: see above
double rtol, double atol: relative and absolute tolerance levels desired
Output:
double *result: result (nt blocks of size 2dim)
double *result: result (nt blocks of size dim)
int *err: error: -10 if interrupted by CTRL-C (SIGINT)
*/
void bovy_rk4(void (*func)(double t, double *q, double *a,
Expand Down
2 changes: 2 additions & 0 deletions tests/test_orbit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,13 +155,15 @@ def test_energy_jacobi_conservation():
tol= {}
tol['default']= -10.
tol['DoubleExponentialDiskPotential']= -6. #these are more difficult
tol['fullyRotatedTriaxialNFWPotential']= -9. #these are more difficult
jactol= {}
jactol['default']= -10.
jactol['RazorThinExponentialDiskPotential']= -9. #these are more difficult
jactol['DoubleExponentialDiskPotential']= -6. #these are more difficult
jactol['mockFlatDehnenBarPotential']= -8. #these are more difficult
jactol['mockFlatDehnenSmoothBarPotential']= -8. #these are more difficult
jactol['mockMovingObjectLongIntPotential']= -8. #these are more difficult
jactol['fullyRotatedTriaxialNFWPotential']= -9. #these are more difficult
jactol['mockSlowFlatEllipticalDiskPotential']= -6. #these are more difficult (and also not quite conserved)
jactol['mockSlowFlatSteadyLogSpiralPotential']= -8. #these are more difficult (and also not quite conserved)
jactol['mockSlowFlatDehnenSmoothBarPotential']= -8. #these are more difficult (and also not quite conserved)
Expand Down

0 comments on commit 444aa14

Please sign in to comment.