/* --------------------------------------

Duffing eqn solver

gcc $CFLAGS -Wall -Wextra -o duffing2x2 duffing2x2.c -lm -lgsl -lgslcblas

dx_i/dt = f_i

Jacobian: df/dx, df/dt

integrate from t0 to tmax in steps of tinc

Duffing eqn: x'' + a*x' + b*x^3 + c*x = f*cos(w*t)

./duffing a b c w f x0 x'0 t0 tmax tstep solver > file.dat
./duffing 0.5 1.0 -1.0 1.0 0.39 0.0 1.0 0.0 160.0 0.05 2 > file.dat

gnuplot (plot phase space)
set term wxt size 1200,1200 background '#c0c0c0' linewidth 1
set size nosquare
plot [-2:2] [-2:2] "duffing.dat" using 2:3 with lines

----------------------------------------*/
#define _GNU_SOURCE
#define HAVE_INLINE
#define GSL_RANGE_CHECK_OFF
#define GSL_DISABLE_DEPRECATED

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

// global structure for ODE system parameters
// alter as needed
typedef struct parameters {
double par1;
double par2;
double par3;
double par4;
double par5;
} parameters;

// -----------------------------
// GSL function for ODE system
// -----------------------------
int ode_eqn (double t, const double x[], double f[], void *params)
{ 
// cast void pointer to type parameters
double p1 = ((parameters *)params)->par1;
double p2 = ((parameters *)params)->par2;
double p3 = ((parameters *)params)->par3;
double p4 = ((parameters *)params)->par4;
double p5 = ((parameters *)params)->par5;

// define GSL functions, i.e. the 2D non-autonoumous Duffing system
f[0]=x[1];
f[1]=-p3*x[0]-p2*x[0]*x[0]*x[0]-p1*x[1]+p5*cos(p4*t);

return GSL_SUCCESS; 
}
// ----------------------------

//-----------------------------
// GSL function for Jacobian
// ----------------------------
int ode_jac (double t, const double x[], double *dfdx, double dfdt[], void *params)
{ 
// cast void pointer to type parameters
double p1 = ((parameters *)params)->par1;
double p2 = ((parameters *)params)->par2;
double p3 = ((parameters *)params)->par3;
double p4 = ((parameters *)params)->par4;
double p5 = ((parameters *)params)->par5;

*(dfdx)=0.0;
*(dfdx+1)=-p3-3.0*p2*x[0]*x[0];

*(dfdx+2)=1.0;
*(dfdx+3)=-p1;

dfdt[0]=0.0;
dfdt[1]=-p4*p5*sin(p4*t);
 
return GSL_SUCCESS; 
}
// -----------------------------


// -----------------------------
// main()
// -----------------------------

int main (int argc, char **argv)
{ 

// initial values for system in t0, x[]
double t, tinc, ti, t0, tmax, x[2];
int status, solver;
// global structure param_list contains ODE system parameters, coefficients
parameters param_list;
// vars/parameters for GSL ODE solver
double hstart, eps_abs, eps_rel, a_y, a_dydt, scale_abs[2]; 
// gsl_odeiv2_system sys is defined below
gsl_odeiv2_driver *drv;

// read initial conditions and time steps

param_list.par1=atof(argv[1]);
param_list.par2=atof(argv[2]);
param_list.par3=atof(argv[3]);
param_list.par4=atof(argv[4]);
param_list.par5=atof(argv[5]);
x[0]=atof(argv[6]);
x[1]=atof(argv[7]);
t0=atof(argv[8]);
tmax=atof(argv[9]);
tinc=atof(argv[10]);
solver=atoi(argv[11]);

// set up ode solvers

gsl_odeiv2_system sys = { ode_eqn, ode_jac, 2, &param_list };

hstart=1e-6;
eps_abs=1e-10;
eps_rel=1e-10;
a_y=1.0;
a_dydt=0.0;

// these weight values allow fine tuning of the convergence but
// we leave them at 1.0 here
scale_abs[0]=1.0;
scale_abs[1]=1.0;

// two choices: 1=rk8pd, 2=rk4imp
// but we always choose 2
if(solver==1) {
drv = gsl_odeiv2_driver_alloc_scaled_new(&sys, gsl_odeiv2_step_rk8pd, hstart, eps_abs, eps_rel, a_y, a_dydt, scale_abs);
}
else if(solver==2) {
drv = gsl_odeiv2_driver_alloc_scaled_new(&sys, gsl_odeiv2_step_rk4imp, hstart, eps_abs, eps_rel, a_y, a_dydt, scale_abs);
}
else {
printf("%s\n", "Error: No solver selected");
return -1;
} 

// set up times and open output file

// output initial values
fprintf(stdout, "%.5e %.5e %.5e\n", t0, x[0], x[1]);

t=t0; ti=t0+tinc;
// simulate next steps from t to ti
do {
status = gsl_odeiv2_driver_apply(drv, &t, ti, x);
if (status != GSL_SUCCESS)
{
printf("Error, return value=%d %s\n", status, gsl_strerror(status));
break;
} 
// output data for time step
fprintf(stdout, "%.5e %.5e %.5e\n", t, x[0], x[1]);
ti=t+tinc;
}
while(ti<=tmax);

gsl_odeiv2_driver_free(drv);
return 0;
}