swirling_strength.h

 
#ifndef VA_SWIRLING_STRENGTH_H
#define VA_SWIRLING_STRENGTH_H
 
#include "common_defs.h"
#include "linalg3.h"
 
/***************************************************************************/
VA_DEVICE_FUN void valib_swirling_strength(VA_REAL *A,
                                           VA_REAL *lambda_ci,
                                           VA_REAL *lambda_cr)
{
    VA_REAL Q, R, delta;
 
    // determine the discriminant of the cubic characteristic equation
    valib_delta(A, &delta);
 
    // find the second invariant
    valib_second_invariant3(A, &Q);
 
    // find the third invariant (determinant)
    valib_determinant3(A, &R);
 
    if (delta <= 0) {
        // only real eigenvalues exist
        *lambda_ci = 0.;
        *lambda_cr = 0.;
    }
    else {
        // one real eigenvalue and a pair of complex conjugate eigenvalues
        VA_REAL aux = sqrt(delta);
 
        VA_REAL aux2;
        aux2 = 0.5*R + aux;
        VA_REAL g1 = valib_sign(1., aux2) * pow(fabs( aux2 ), 1./3.);
 
        aux2 = 0.5*R - aux;
        VA_REAL g2 = valib_sign(1., aux2) * pow(fabs( aux2 ), 1./3.);
 
        *lambda_ci = (g1 - g2) * sqrt(3.) / 2.;
        *lambda_cr = -(g1 + g2) / 2.;
    }
}
 
#endif // VA_SWIRLING_STRENGTH_H