#include <math.h>
#include <stdlib.h>

/*normalizes anisotropy axes*/
void axes(double* ap, double* scf, int N)
{
    double sclpdt, r, eps = 1.E-6;
    int i, j, k;

    for (k = 0; k < N; k++) {

        r = sqrt(ap[9 * k] * ap[9 * k] + ap[9 * k + 1] * ap[9 * k + 1] + ap[9 * k + 2] * ap[9 * k + 2]);
        ap[9 * k] /= r;
        ap[9 * k + 1] /= r;
        ap[9 * k + 2] /= r;
        sclpdt = ap[9 * k] * ap[9 * k + 3] + ap[9 * k + 1] * ap[9 * k + 4] + ap[9 * k + 2] * ap[9 * k + 5];
        if (sclpdt > eps) {
            printf("Non orthogonal axes");
            exit;
        } else {
            r = sqrt(ap[9 * k + 3] * ap[9 * k + 3] + ap[9 * k + 4] * ap[9 * k + 4] + ap[9 * k + 5] * ap[9 * k + 5]);
            ap[9 * k + 3] /= r;
            ap[9 * k + 4] /= r;
            ap[9 * k + 5] /= r;
            ap[9 * k + 6] = ap[9 * k + 1] * ap[9 * k + 5] - ap[9 * k + 2] * ap[9 * k + 4];
            ap[9 * k + 7] = ap[9 * k + 2] * ap[9 * k + 3] - ap[9 * k] * ap[9 * k + 5];
            ap[9 * k + 8] = ap[9 * k] * ap[9 * k + 4] - ap[9 * k + 1] * ap[9 * k + 3];
            for (i = 0; i < 3; i++) {
                for (j = 0; j < 3; j++) {
                    if (scf[3 * k + j] == 0.)
                        scf[3 * k + j] = 1.;
                    ap[9 * k + 3 * j + i] /= scf[3 * k + j];
                }
            }
        }
    }

    return;
}