#include <stdio.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

#define MASS 1.0
#define OMEGA 10.0
#define XINIT 1.0
#define PINIT 0.0

#define STEPSIZE 1E-8
#define BIGSTEP 1E-3

#define ABSERR 1e-6    
#define RELERR 1e-8
#define YERR 0.5
#define YPERR 0.2

int func(double t, const double y[], double dydt[], void *params)
       {
         double mass,w;
        
        mass=*(double *)params;
        params++;
        w=*(double *)params;
         printf("\n %lf %lf",mass,w);

        dydt[0]=y[1];
        dydt[1]=-mass*w*w*y[0];
        return GSL_SUCCESS;
       }

int jac(double t, const double y[], double *dfdy, double *dfdt, void *params)
       {
        double mass,w;
        gsl_matrix_view dfdy_mat=gsl_matrix_view_array(dfdy,2,2);
        gsl_matrix *m=&dfdy_mat.matrix;
       
        mass=*(double *)params;
        params++;
        w=*(double *)params;

        gsl_matrix_set(m,0,0,0.0);
        gsl_matrix_set(m,0,1,1/mass);
        gsl_matrix_set(m,1,0,-mass*w*w);
        gsl_matrix_set(m,1,1,0.0);
        dfdt[0]=0.0;
        dfdt[1]=0.0;
        return GSL_SUCCESS;
       }

 int integrate(double initial, double final, void *parameters, double *y)
              {
                const gsl_odeiv_step_type *T = gsl_odeiv_step_rk8pd;
                gsl_odeiv_step *s=gsl_odeiv_step_alloc(T,2);
                gsl_odeiv_control *c=gsl_odeiv_control_standard_new(ABSERR,RELERR,YERR,YPERR);
                gsl_odeiv_evolve *e=gsl_odeiv_evolve_alloc(2);
                double params[2];
                gsl_odeiv_system sys = {func,jac,2,params};
                double t=initial, t1=final;
                double stepsize=STEPSIZE;

                params[0]=*(double *)parameters;
                parameters++;
                params[1]=*(double *)parameters;

                while(t<t1)
                     {
                       int status=gsl_odeiv_evolve_apply(e,c,s,&sys,&t,t1,&stepsize,y);
                       if(status!=GSL_SUCCESS)
                          break;
                     }
                gsl_odeiv_evolve_free(e);
                gsl_odeiv_control_free(c);
                gsl_odeiv_step_free(s);
                return GSL_SUCCESS;
               }


int main(void)
     {
       double t, period=2.0*M_PI/OMEGA;
       double params[2]={MASS,OMEGA};
       double y[2]={XINIT,PINIT};
       double step=BIGSTEP;
       FILE *out;
       out=fopen("out.dat","w");
       int status;
       for(t=0.0;t<=period;t=t+step)
          {
           status=integrate(t,t+step,params,y);
           if(status!=GSL_SUCCESS)
             {
              printf("GSL execution of rk failed at t = %lf",t);
              break;
             }
           fprintf(out,"\n%lf %lf",y[0],y[1]);
         }
      return 0;
     }