shm_2params

#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.2int 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;     }