/* modsubs.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"

/* Subroutine */ int cross_(real *b, real *c__, real *a)
{
/*     ======================= */


    /* 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;
} /* cross_ */


doublereal dot_(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;
} /* dot_ */

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

/* 	======================== */

/* 	MULTIPLY 2 3X3 MATRICES */

/* 	A=BC */

    /* 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) {
/* L3: */
		s += b[i__ + k * 3] * c__[k + j * 3];
	    }
/* L2: */
	    a[i__ + j * 3] = s;
	}
/* L1: */
    }
    return 0;
} /* matmul_ */

/* Subroutine */ int matvec_(real *v, real *a, real *b)
{
    static integer i__, j;
    static real s;

/* 	======================== */

/* 	POST-MULTIPLY A 3X3 MATRIX BY A VECTOR */

/* 	V=AB */

    /* Parameter adjustments */
    --b;
    a -= 4;
    --v;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
	s = 0.f;
	for (j = 1; j <= 3; ++j) {
/* L2: */
	    s += a[i__ + j * 3] * b[j];
	}
/* L1: */
	v[i__] = s;
    }
    return 0;
} /* matvec_ */

/* Subroutine */ int minv_(real *a, real *b, real *d__)
{
    static real c__[9]	/* was [3][3] */;
    static integer i__, j;
    extern doublereal dot_(real *, real *);
    extern /* Subroutine */ int cross_(real *, real *, real *);

/* 	====================== */

/* 	INVERT A GENERAL 3X3 MATRIX AND RETURN DETERMINANT IN D */

/* 	A=(B)-1 */

    /* Parameter adjustments */
    b -= 4;
    a -= 4;

    /* Function Body */
    cross_(&b[7], &b[10], c__);
    cross_(&b[10], &b[4], &c__[3]);
    cross_(&b[4], &b[7], &c__[6]);
    *d__ = dot_(&b[4], c__);

/* 	TEST DETERMINANT */

    if (dabs(*d__) > 1e-30f) {
	goto L10;
    }
    *d__ = 0.f;
    return 0;

/* 	DETERMINANT IS NON-ZERO */

L10:
    for (i__ = 1; i__ <= 3; ++i__) {
	for (j = 1; j <= 3; ++j) {
/* L20: */
	    a[i__ + j * 3] = c__[j + i__ * 3 - 4] / *d__;
	}
    }
    return 0;
} /* minv_ */


/* Subroutine */ int scalev_(real *a, real *x, real *b)
{
    static integer i__;

/*     ======================== */

/*     SCALE A VECTOR B WITH SCALAR X AND PUT RESULT IN A */

    /* Parameter adjustments */
    --b;
    --a;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
	a[i__] = b[i__] * *x;
/* L10: */
    }
    return 0;
} /* scalev_ */

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

/* 	====================== */

/* 	TRANSPOSE A 3X3 MATRIX */

/* 	A=BT */

    /* Parameter adjustments */
    b -= 4;
    a -= 4;

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

/* Subroutine */ int unit_(real *v)
{
    /* System generated locals */
    real r__1, r__2, r__3;

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

    /* Local variables */
    static integer i__;
    static real vmod;

/*     ================= */


/* ****VECTOR V REDUCED TO UNIT VECTOR */
    /* Parameter adjustments */
    --v;

    /* Function Body */
/* Computing 2nd power */
    r__1 = v[1];
/* Computing 2nd power */
    r__2 = v[2];
/* Computing 2nd power */
    r__3 = v[3];
    vmod = r__1 * r__1 + r__2 * r__2 + r__3 * r__3;
    vmod = sqrt(vmod);
    for (i__ = 1; i__ <= 3; ++i__) {
/* L1: */
	v[i__] /= vmod;
    }
    return 0;
} /* unit_ */

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

/*     ===================== */


/*     SUBTRACT TWO VECTORS */

    /* Parameter adjustments */
    --c__;
    --b;
    --a;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
	a[i__] = b[i__] - c__[i__];
/* L10: */
    }
    return 0;
} /* vdif_ */

/* Subroutine */ int vset_(real *a, real *b)
{
    static integer i__;

/*     ==================== */


/*     MOVE A VECTOR FROM B TO A */

    /* Parameter adjustments */
    --b;
    --a;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
/* L10: */
	a[i__] = b[i__];
    }
    return 0;
} /* vset_ */

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

/*     ====================== */


/*     ADD TWO VECTORS */

    /* Parameter adjustments */
    --c__;
    --b;
    --a;

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