Skip to content

Commit

Permalink
Merge branch 'master' into fix361
Browse files Browse the repository at this point in the history
  • Loading branch information
jobovy committed Nov 28, 2018
2 parents 59eec24 + f3b6dcc commit 99981c5
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 20 deletions.
3 changes: 3 additions & 0 deletions HISTORY.txt
Expand Up @@ -20,6 +20,9 @@ v1.5 (????-??-??)

- DehnenSmoothWrapperPotential can now decay rather than grow a
potential by setting ``decay=True``.

- Allow orbit integrations to be KeyboardInterrupted on Windows as well
(#362 by Henry Leung)

v1.4 (2018-09-09)
=================
Expand Down
26 changes: 15 additions & 11 deletions galpy/util/bovy_rk.c
Expand Up @@ -34,9 +34,7 @@ POSSIBILITY OF SUCH DAMAGE.
#include <math.h>
#include <bovy_symplecticode.h>
#include <bovy_rk.h>
#ifndef _WIN32
#include "signal.h"
#endif
#define _MAX_STEPCHANGE_POWERTWO 3.
#define _MIN_STEPCHANGE_POWERTWO -3.
#define _MAX_STEPREDUCE 10000.
Expand Down Expand Up @@ -99,6 +97,8 @@ void bovy_rk4(void (*func)(double t, double *q, double *a,
memset(&action, 0, sizeof(struct sigaction));
action.sa_handler= handle_sigint;
sigaction(SIGINT,&action,NULL);
#else
if (SetConsoleCtrlHandler(CtrlHandler, TRUE)) {}
#endif
for (ii=0; ii < (nt-1); ii++){
if ( interrupted ) {
Expand Down Expand Up @@ -201,6 +201,8 @@ void bovy_rk6(void (*func)(double t, double *q, double *a,
memset(&action, 0, sizeof(struct sigaction));
action.sa_handler= handle_sigint;
sigaction(SIGINT,&action,NULL);
#else
if (SetConsoleCtrlHandler(CtrlHandler, TRUE)) {}
#endif
for (ii=0; ii < (nt-1); ii++){
if ( interrupted ) {
Expand Down Expand Up @@ -277,14 +279,14 @@ void bovy_rk6_onestep(void (*func)(double t, double *q, double *a,
func(tn+2.*dt/3.,ynk,a,nargs,potentialArgs);
for (ii=0; ii < dim; ii++) *(yn1+ii) += 81. * dt * *(a+ii) / 120.;
for (ii=0; ii < dim; ii++) *(k3+ii)= dt * *(a+ii);
for (ii=0; ii < dim; ii++) *(ynk+ii)= *(yn+ii) + ( *(k1+ii)
for (ii=0; ii < dim; ii++) *(ynk+ii)= *(yn+ii) + ( *(k1+ii)
+ 4. * *(k2+ii)
- *(k3+ii))/12.;
//calculate k4
func(tn+dt/3.,ynk,a,nargs,potentialArgs);
for (ii=0; ii < dim; ii++) *(yn1+ii) += 81.* dt * *(a+ii) / 120.;
for (ii=0; ii < dim; ii++) *(k4+ii)= dt * *(a+ii);
for (ii=0; ii < dim; ii++) *(ynk+ii)= *(yn+ii) + ( -*(k1+ii)
for (ii=0; ii < dim; ii++) *(ynk+ii)= *(yn+ii) + ( -*(k1+ii)
+ 18. * *(k2+ii)
- 3. * *(k3+ii)
-6.* *(k4+ii))/16.;
Expand All @@ -307,7 +309,7 @@ void bovy_rk6_onestep(void (*func)(double t, double *q, double *a,
-64. * *(k5+ii))/44.;
//calculate k7
func(tn+dt,ynk,a,nargs,potentialArgs);
for (ii=0; ii < dim; ii++) *(yn1+ii) += 11.* dt * *(a+ii) / 120.;
for (ii=0; ii < dim; ii++) *(yn1+ii) += 11.* dt * *(a+ii) / 120.;
//yn1 is new value
}

Expand Down Expand Up @@ -360,10 +362,10 @@ double rk4_estimate_step(void (*func)(double t, double *y, double *a,int nargs,
err+= exp(2.*log(fabs(*(y1+ii)-*(y2+ii)))-2.* *(scale+ii));
}
err= sqrt(err/dim);
if ( ceil(pow(err,1./5.)) > 1.
if ( ceil(pow(err,1./5.)) > 1.
&& init_dt / dt * ceil(pow(err,1./5.)) < _MAX_DT_REDUCE)
dt/= ceil(pow(err,1./5.));
else
else
break;
}
//free what we allocated
Expand All @@ -378,7 +380,7 @@ double rk4_estimate_step(void (*func)(double t, double *y, double *a,int nargs,
//printf("%f\n",dt);
//fflush(stdout);
return dt;
}
}
double rk6_estimate_step(void (*func)(double t, double *y, double *a,int nargs, struct potentialArg *),
int dim, double *yo,
double dt, double *t,
Expand Down Expand Up @@ -436,10 +438,10 @@ double rk6_estimate_step(void (*func)(double t, double *y, double *a,int nargs,
err+= exp(2.*log(fabs(*(y1+ii)-*(y2+ii)))-2.* *(scale+ii));
}
err= sqrt(err/dim);
if ( ceil(pow(err,1./7.)) > 1.
if ( ceil(pow(err,1./7.)) > 1.
&& init_dt / dt * ceil(pow(err,1./7.)) < _MAX_DT_REDUCE)
dt/= ceil(pow(err,1./7.));
else
else
break;
}
//free what we allocated
Expand All @@ -459,7 +461,7 @@ double rk6_estimate_step(void (*func)(double t, double *y, double *a,int nargs,
//printf("%f\n",dt);
//fflush(stdout);
return dt;
}
}
/*
Runge-Kutta Dormand-Prince 5/4 integrator
Usage:
Expand Down Expand Up @@ -525,6 +527,8 @@ void bovy_dopr54(void (*func)(double t, double *q, double *a,
memset(&action, 0, sizeof(struct sigaction));
action.sa_handler= handle_sigint;
sigaction(SIGINT,&action,NULL);
#else
if (SetConsoleCtrlHandler(CtrlHandler, TRUE)) {}
#endif
for (ii=0; ii < (nt-1); ii++){
if ( interrupted ) {
Expand Down
27 changes: 24 additions & 3 deletions galpy/util/bovy_symplecticode.c
Expand Up @@ -34,16 +34,31 @@ POSSIBILITY OF SUCH DAMAGE.
#include <math.h>
#include <bovy_symplecticode.h>
#define _MAX_DT_REDUCE 10000.
//CTRL-C only works on UNIX systems with signal library
#ifndef _WIN32
#include "signal.h"
volatile sig_atomic_t interrupted= 0;

// handle CTRL-C differently on UNIX systems and Windows
#ifndef _WIN32
void handle_sigint(int signum)
{
interrupted= 1;
}
#else
int interrupted= 0;
#include <windows.h>
BOOL WINAPI CtrlHandler(DWORD fdwCtrlType)
{
switch (fdwCtrlType)
{
// Handle the CTRL-C signal.
case CTRL_C_EVENT:
interrupted= 1;
// needed to avoid other control handlers like python from being called before us
return TRUE;
default:
return FALSE;
}
}

#endif
static inline void leapfrog_leapq(int dim, double *q,double *p,double dt,
double *qn){
Expand Down Expand Up @@ -124,6 +139,8 @@ void leapfrog(void (*func)(double t, double *q, double *a,
memset(&action, 0, sizeof(struct sigaction));
action.sa_handler= handle_sigint;
sigaction(SIGINT,&action,NULL);
#else
if (SetConsoleCtrlHandler(CtrlHandler, TRUE)){}
#endif
for (ii=0; ii < (nt-1); ii++){
if ( interrupted ) {
Expand Down Expand Up @@ -242,6 +259,8 @@ void symplec4(void (*func)(double t, double *q, double *a,
memset(&action, 0, sizeof(struct sigaction));
action.sa_handler= handle_sigint;
sigaction(SIGINT,&action,NULL);
#else
if (SetConsoleCtrlHandler(CtrlHandler, TRUE)) {}
#endif
for (ii=0; ii < (nt-1); ii++){
if ( interrupted ) {
Expand Down Expand Up @@ -395,6 +414,8 @@ void symplec6(void (*func)(double t, double *q, double *a,
memset(&action, 0, sizeof(struct sigaction));
action.sa_handler= handle_sigint;
sigaction(SIGINT,&action,NULL);
#else
if (SetConsoleCtrlHandler(CtrlHandler, TRUE)) {}
#endif
for (ii=0; ii < (nt-1); ii++){
if ( interrupted ) {
Expand Down
9 changes: 3 additions & 6 deletions galpy/util/bovy_symplecticode.h
Expand Up @@ -34,23 +34,20 @@ POSSIBILITY OF SUCH DAMAGE.
#ifdef __cplusplus
extern "C" {
#endif
#ifndef _WIN32
#include "signal.h"
#endif
#include <galpy_potentials.h>
/*
Global variables
*/
#ifndef _WIN32
extern volatile sig_atomic_t interrupted;
#else
extern int interrupted;
#endif
/*
Function declarations
*/
#ifndef _WIN32
void handle_sigint(int);
#else
#include "windows.h"
BOOL WINAPI CtrlHandler(DWORD fdwCtrlType);
#endif
void leapfrog(void (*func)(double, double *, double *,
int, struct potentialArg *),
Expand Down

0 comments on commit 99981c5

Please sign in to comment.