/*
 * cpldosc_fastvf_main.c
 *
 */


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

#include "cpldosc_fast_proto.h"


int main (int argc, char *argv[])
    {
    double y[2],p[5];
    double fastvf[2];

    if (argc != 9)
        {
        fprintf(stderr,"use: %s v1 v2 q1 q2 sigma1 sigma2 omega time_duration\n",argv[0]);
        exit(-1);
        }

    double v1     = atof(argv[1]);
    double v2     = atof(argv[2]);
    double q1     = atof(argv[3]);
    double q2     = atof(argv[4]);
    double sigma1 = atof(argv[5]);
    double sigma2 = atof(argv[6]);
    double omega  = atof(argv[7]);
    double t1     = atof(argv[8]);

    const gsl_odeiv_step_type * T  = gsl_odeiv_step_rk8pd;
    gsl_odeiv_step * step    = gsl_odeiv_step_alloc(T, 2);
    gsl_odeiv_control * control = gsl_odeiv_control_y_new(1e-12, 0.0);
    gsl_odeiv_evolve * evolve  = gsl_odeiv_evolve_alloc(2);

    p[0]  =  sigma1;
    p[1]  =  sigma2;
    p[2]  =  omega;
    p[3]  =  q1;
    p[4]  =  q2;

    y[0] = v1;
    y[1] = v2;

    gsl_odeiv_system sys = {cpldosc_fast_vf, cpldosc_fast_jac, 2, &(p[0])};

    double t  = 0.0;

    double h = 1e-6;
    double ztol = 1.0e-10;
    while (t < t1)
        {
        int status = gsl_odeiv_evolve_apply(evolve, control, step, &sys, &t, t1, &h, y);
        if (status != GSL_SUCCESS)
            {
            fprintf(stderr,"status=%d\n",status);
            break;
            }
        printf("%16.12e %16.12e %16.12e %16.12e %16.12e\n", t, y[0], y[1],q1,q2);
        cpldosc_fast_vf(0.0,y,fastvf,p);
        double m = sqrt(fastvf[0]*fastvf[0]+fastvf[1]*fastvf[1]);
        if (m < ztol)
            break;
        }

    gsl_odeiv_evolve_free(evolve);
    gsl_odeiv_control_free(control);
    gsl_odeiv_step_free(step);
    return 0;
    }