Rotate a rank 2 tensor of dimension 3 x 3 about a given axis by a given angle

rotate_rank2_tensor.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 3

/*--- assign some values to our input matrix ---*/
double** fill_rank2_tensor(void)
{
  int i, j;
  double** t = (double**) malloc(sizeof(double*) * N);
  for(i = 0; i < N; ++i){
    t[i] = (double*) malloc(sizeof(double) * N);
    for(j = 0; j < N; ++j){
      t[i][j] = i * N + j;
    }
  }
  return t;
}

/*--- Axis vector if not normalized, should be subjected to 
 * normalization ---*/
void normalize_vec(double* n)
{
  double mag = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]) + 1e-12;
  n[0] = n[0]/mag;
  n[1] = n[1]/mag;
  n[2] = n[2]/mag;
  return;
}

/*--- dynamically create a matrix of size 3x3 ---*/
double** alloc_mat(void)
{
  int i, j;
  double** r = (double**) malloc(sizeof(double*) * N);
  for(i = 0; i < N; ++i){
    r[i] = (double*) malloc(sizeof(double) * N);
  }
  return r;
}

/*--- dynamically create a rotation matrix of size 3x3 ---*/
double** alloc_rot_mat(void)
{
  return alloc_mat();
}

/*--- fill the entries for the rotation matrix ---*/
/*--- See https://en.wikipedia.org/wiki/Rotation_matrix --*/
void fill_rot_mat(double** r, double* n, double th)
{
 double c, s, cm, n1, n2, n3; 

 c = cos(th);
 cm = 1 - c;
 s = sin(th);
 n1 = n[0];
 n2 = n[1];
 n3 = n[2];

 r[0][0] = c + (n1 * n1 * cm);
 r[0][1] = (n1 * n2 * cm) - (n3 * s);
 r[0][2] = (n1 * n3 * cm) + (n2 * s);

 r[1][0] = (n1 * n2 * cm) + (n3 * s);
 r[1][1] = c + (n2 * n2 * cm);
 r[1][2] = (n2 * n3 * cm) - (n1 * s);
 
 r[2][0] = (n1 * n3 * cm) - (n2 * s);
 r[2][1] = (n2 * n3 * cm) + (n1 * s);
 r[2][2] = c + (n3 * n3 * cm);
}

/*--- create and initialize the rotation matrix ---*/
double** get_rot_mat(double* n, double th)
{
  double** r = alloc_rot_mat();
  fill_rot_mat(r, n, th);
  return r;
}

/*--- Free the heap memory allocated for the matrix ---*/
void free_mat(double** t)
{
  int i;
  for(i = 0; i < N; ++i){
    free(t[i]);
  }
  free(t);
}

/*--- Free the memory for rotation matrix ---*/
void free_rot_mat(double** rm)
{
  free_mat(rm);
}

/*--- does matrix multiplication ---*/
double** multiply_matrices(double** mat1, double** mat2)
{
  int i, j, k;
  double** temp = alloc_mat();

  for(i = 0; i < N; ++i){
    for(j = 0; j < N; ++j){
      temp[i][j] = 0.0;
      for(k = 0; k < N; ++k){
        temp[i][j] += mat1[i][k] * mat2[k][j];
      }
    }
  }
  return temp;
}

/*--- creates a copy of the transposed matrix ---*/
double** get_transpose(double** mat)
{
  int i, j;
  double** res = alloc_mat();
  for(i = 0; i < N; ++i){
    for(j = 0; j < N; ++j){
      res[i][j] = mat[j][i];
    }
  }
  return res;
}

/*--- copy the values of a src matrix to dest matrix ---*/
void copy_mat(double** src, double** dest)
{
  int i, j;
  for(i = 0; i < N; ++i){
    for(j = 0; j < N; ++j){
      dest[i][j] = src[i][j];
    }
  }
}

/*--- function that oversees the rotation operation ---*/
void rotate_rank2_tensor(double** tns, double* n, double th)
{
  normalize_vec(n);
  double** rm = get_rot_mat(n, th);
  double** rm_t = get_transpose(rm);
  double** temp = multiply_matrices(tns, rm_t);
  double** tns_rotated = multiply_matrices(rm, temp);
  copy_mat(tns_rotated, tns);
  free_mat(tns_rotated);
  free_mat(temp);
  free_rot_mat(rm_t);
  free_rot_mat(rm);
}

/*--- free rank2 tensor heap memory allocated ---*/
void free_rank2_tensor(double** t)
{
  free_mat(t);
}

/*--- prints a matrix ---*/
void print_rank2_tensor(double** t){
  int i, j;
  for(i = 0; i < N; ++i){
    for(j = 0; j < N; ++j){
      printf("%0.2lf ", t[i][j]);
    }
    printf("\n");
  }
}

int main(){
  /*-- axis about which one should rotate the tensor (direction cosines) --*/
  double axis_vec[3] = {1.0, 0.0, 0.0}; 

  /*-- angle of rotation about the axis --*/
  double angle = 30 * M_PI/180; // in radians
                                //
  /*-- tns is our tensor to which we fill some values --*/
  double** tns = (double**) fill_rank2_tensor();

  /*-- print the input tensor --*/
  printf("Initial:\n");
  print_rank2_tensor(tns);

  /*-- rotate the tensor about the given axis by the given angle--*/
  rotate_rank2_tensor(tns, axis_vec, angle);

  /*-- print the output rotated tensor --*/
  printf("\nFinal:\n");
  print_rank2_tensor(tns);

  /*-- free the memory allocated for the tensor --*/
  free_rank2_tensor(tns);

  return 0;
}