/* function loaded by the rungekutte~ object to describe a forced
oscillator.  See the README in this directory. */

    /* parameters and state variables are given symbolic names here to make
    the code below more readable */
#define PARAM_K (params[0])             /* resonant angular freq */
#define PARAM_FORCING (params[1])       /* forcing to unit length */
#define PARAM_COUPLING (params[2])      /* coupling input */
#define STATE_X (state[0])              /* 2 state variables x, y */
#define STATE_Y (state[1])
#define DERIVATIVE_X (result[0])        /* 2 derivatives x-dot, y-dot */
#define DERIVATIVE_Y (result[1])

int rungekutte_nparam = 3;
int rungekutte_nstate = 2;

void rungekutte_derivative(float *result, float *state, float *params)
{
    float correction = PARAM_FORCING * (1 - STATE_X*STATE_X - STATE_Y*STATE_Y);
    if (correction < -PARAM_FORCING)
        correction = -PARAM_FORCING;
    DERIVATIVE_X = - PARAM_K * STATE_Y + correction * STATE_X + PARAM_COUPLING;
    DERIVATIVE_Y = PARAM_K * STATE_X + correction * STATE_Y;
}