/* suv.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

doublereal ccuij_(real *uu, real *vv)
{
    /* System generated locals */
    real ret_val;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static real ui[9]	/* was [3][3] */, vi[9]	/* was [3][3] */, wi[9]	/* 
	    was [3][3] */, ww[9]	/* was [3][3] */, dui, dvi, dwi;
    extern /* Subroutine */ int m3add_(real *, real *, real *);
    extern doublereal v3inv_(real *, real *);

/*       =============================================================== */
/* 	Correlation coefficient between two (3x3) anisotropic */
/* 	displacement matrices Uij and Vij. */
/* 	See Merritt (1999) Acta Crystallographica D55, 1997-2004. */
/*       =============================================================== */
    /* Parameter adjustments */
    vv -= 4;
    uu -= 4;

    /* Function Body */
    dui = 1.f / v3inv_(ui, &uu[4]);
    dvi = 1.f / v3inv_(vi, &vv[4]);
    if (dui <= 0.f || dvi <= 0.f) {
	ret_val = 0.f;
	return ret_val;
    }
    m3add_(wi, ui, vi);
    dwi = v3inv_(ww, wi);
    ret_val = sqrt(sqrt(dui * dvi)) / sqrt(dwi * .125f);
    return ret_val;
} /* ccuij_ */

doublereal suv_(real *u, real *v)
{
    /* System generated locals */
    real ret_val;

    /* Local variables */
    static integer i__, j;
    static real uu[9]	/* was [3][3] */, vv[9]	/* was [3][3] */, ww[9]	/* 
	    was [3][3] */, ueq, veq, uiso[9]	/* was [3][3] */, viso[9]	
	    /* was [3][3] */;
    extern doublereal ccuij_(real *, real *);
    static real suv_bot__, suv_top__;

/*       =============================================================== */
/* 	Normalized correlation coefficient ("similarity") between two */
/* 	anisotropic displacement vectors Uij and Vij. */
/* 	See Merritt (1999) Acta Crystallographica D55, 1997-2004. */
/*       =============================================================== */
    /* Parameter adjustments */
    --v;
    --u;

    /* Function Body */
    uu[0] = u[1];
    uu[4] = u[2];
    uu[8] = u[3];
    uu[3] = u[4];
    uu[1] = u[4];
    uu[6] = u[5];
    uu[2] = u[5];
    uu[7] = u[6];
    uu[5] = u[6];
    vv[0] = v[1];
    vv[4] = v[2];
    vv[8] = v[3];
    vv[3] = v[4];
    vv[1] = v[4];
    vv[6] = v[5];
    vv[2] = v[5];
    vv[7] = v[6];
    vv[5] = v[6];
    ueq = (uu[0] + uu[4] + uu[8]) / 3.f;
    veq = (vv[0] + vv[4] + vv[8]) / 3.f;
    for (i__ = 1; i__ <= 3; ++i__) {
	for (j = 1; j <= 3; ++j) {
	    ww[i__ + j * 3 - 4] = vv[i__ + j * 3 - 4] * ueq / veq;
	    uiso[i__ + j * 3 - 4] = 0.f;
	    viso[i__ + j * 3 - 4] = 0.f;
	}
    }
    for (i__ = 1; i__ <= 3; ++i__) {
	uiso[i__ + i__ * 3 - 4] = ueq;
	viso[i__ + i__ * 3 - 4] = veq;
    }
    suv_top__ = ccuij_(uu, ww);
    suv_bot__ = ccuij_(uu, uiso) * ccuij_(vv, viso);
    if (suv_top__ == 0.f || suv_bot__ == 0.f) {
	ret_val = 0.f;
    } else {
	ret_val = suv_top__ / suv_bot__;
    }
    return ret_val;
} /* suv_ */

/* 	======================================================= */
/* 	The rest of the file is just generic 3x3 matrix algebra */
/* 	======================================================= */
/* Subroutine */ int v3cross_(real *b, real *c__, real *a)
{
/* 	CROSS PRODUCT OF TWO VECTORS */
/* 	======================= */
    /* Parameter adjustments */
    --a;
    --c__;
    --b;

    /* Function Body */
    a[1] = b[2] * c__[3] - c__[2] * b[3];
    a[2] = b[3] * c__[1] - c__[3] * b[1];
    a[3] = b[1] * c__[2] - c__[1] * b[2];
    return 0;
} /* v3cross_ */


doublereal v3dot_(real *a, real *b)
{
    /* System generated locals */
    real ret_val;

/* 	DOT PRODUCT OF TWO VECTORS */
/* 	================ */
    /* Parameter adjustments */
    --b;
    --a;

    /* Function Body */
    ret_val = a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
    return ret_val;
} /* v3dot_ */


/* Subroutine */ int m3add_(real *a, real *b, real *c__)
{
    static integer i__, j;

/* 	A=B+C */
/* 	======================== */
    /* Parameter adjustments */
    c__ -= 4;
    b -= 4;
    a -= 4;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
	for (j = 1; j <= 3; ++j) {
	    a[i__ + j * 3] = b[i__ + j * 3] + c__[i__ + j * 3];
	}
    }
    return 0;
} /* m3add_ */


/* Subroutine */ int m3mul_(real *a, real *b, real *c__)
{
    static integer i__, j, k;
    static real s;

/* 	A=B*C */
/* 	======================== */
    /* Parameter adjustments */
    c__ -= 4;
    b -= 4;
    a -= 4;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
	for (j = 1; j <= 3; ++j) {
	    s = 0.f;
	    for (k = 1; k <= 3; ++k) {
		s += b[i__ + k * 3] * c__[k + j * 3];
	    }
	    a[i__ + j * 3] = s;
	}
    }
    return 0;
} /* m3mul_ */


doublereal v3inv_(real *a, real *b)
{
    /* System generated locals */
    real ret_val;

    /* Local variables */
    static real c__[9]	/* was [3][3] */, d__;
    static integer i__, j;
    extern doublereal v3dot_(real *, real *);
    extern /* Subroutine */ int v3cross_(real *, real *, real *);

/* 	INVERT A GENERAL 3X3 MATRIX AND RETURN DETERMINANT */
/* 	A=(B)-1 */
/* 	====================== */
    /* Parameter adjustments */
    b -= 4;
    a -= 4;

    /* Function Body */
    v3cross_(&b[7], &b[10], c__);
    v3cross_(&b[10], &b[4], &c__[3]);
    v3cross_(&b[4], &b[7], &c__[6]);
    d__ = v3dot_(&b[4], c__);
    if (dabs(d__) <= 1e-30f) {
	ret_val = 0.f;
	return ret_val;
    }
    for (i__ = 1; i__ <= 3; ++i__) {
	for (j = 1; j <= 3; ++j) {
	    a[i__ + j * 3] = c__[j + i__ * 3 - 4] / d__;
	}
    }
    ret_val = d__;
    return ret_val;
} /* v3inv_ */