#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;
}