/*
 *  cpldosc_fast_eigvec0.c
 *
 *  Find an eigenvector associated with a zero eigenvalue of the
 *  Jacobian of the fast subsystem.  This code *assumes* that zero
 *  is an eigenvalue.
 *
 *  The vector p contains the parameters and the slow variables:
 *    p[0] = sigma1
 *    p[1] = sigma2
 *    p[2] = omega
 *    p[3] = q1
 *    p[4] = q2
 */


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

#include "cpldosc_fast_proto.h"

int cpldosc_fast_eigvec0(const double v[], const double p[], double ev[])
    {
    double dfdt[2];
    double J[2][2];

    /*
     *  Compute the Jacobian matrix.
     */
    cpldosc_fast_jac(0.0, v, (double *) J, dfdt, (void *) p);

    /*
     *  If the Jacobian is
     *      [ a  b ]
     *      [ c  d ]
     *  then either [b,-a] or [d,-c] could be used as an eigenvector
     *  (if the vector is not the zero vector). The next few lines
     *  of code choose the vector with the larger magnitude.
     */
    double n1 = sqrt(J[0][0]*J[0][0] + J[0][1]*J[0][1]);
    double n2 = sqrt(J[1][0]*J[1][0] + J[1][1]*J[1][1]);
    if (n1 > n2)
        {
        ev[0] =  J[0][1];
        ev[1] = -J[0][0];
        }
    else
        {
        ev[0] =  J[1][1];
        ev[1] = -J[1][0];
        }
    /*
     *  Make the vector a unit vector.
     */
    double n = sqrt(ev[0]*ev[0] + ev[1]*ev[1]);
    ev[0] = ev[0]/n;
    ev[1] = ev[1]/n;
    }