# Gawk library for 3x3 matrices 
# Last edited on 2003-03-09 22:56:18 by stolfi

function r3x3_norm(M,c, i,j,w)
{ # Scale matrix so that {sqrt(M[0,0]^2 + M[1,0]^2 + M[2,0]^2)} is {c}.
  w = c/sqrt(M[0,0]*M[0,0] + M[1,0]*M[1,0] + M[2,0]*M[2,0]);
  for (i = 0; i < 3; i++) 
    { for (j = 0; j < 3; j++)
        { M[i,j] *= w; }
    }
}

function r3x3_mul(M,N,P,   i,j,k,s)
{ for (i = 0; i < 3; i++) 
    { for (j = 0; j < 3; j++)
        { s = 0;
          for (k = 0; k < 3; k++) { s += M[i,k]*N[k,j]; }
          P[i,j] = s;
        }
    }
}

function r3x3_mul_row(u,M,v,   j,k,s)
{ for (j = 0; j < 3; j++)
    { s = 0;
      for (k = 0; k < 3; k++) { s += u[k]*M[k,j]; }
      v[j] = s;
    }
}

function r3x3_inv(Md,Mi,   i,j,d)
{ r3x3_adj(Md,Mi);
  d = Md[0,0]*Mi[0,0] + Md[0,1]*Mi[1,0] + Md[0,2]*Mi[2,0];
  for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { Mi[i,j] /= d; } }
}

function r3x3_adj (A,R)
  { R[0,0] =  (A[1,1] * A[2,2] - A[1,2] * A[2,1]);
    R[0,1] = -(A[0,1] * A[2,2] - A[0,2] * A[2,1]);
    R[0,2] =  (A[0,1] * A[1,2] - A[0,2] * A[1,1]);
              
    R[1,0] = -(A[1,0] * A[2,2] - A[1,2] * A[2,0]);
    R[1,1] =  (A[0,0] * A[2,2] - A[0,2] * A[2,0]);
    R[1,2] = -(A[0,0] * A[1,2] - A[0,2] * A[1,0]);
              
    R[2,0] =  (A[1,0] * A[2,1] - A[1,1] * A[2,0]);
    R[2,1] = -(A[0,0] * A[2,1] - A[0,1] * A[2,0]);
    R[2,2] =  (A[0,0] * A[1,1] - A[0,1] * A[1,0]);
  }

function r3x3_print(M,   i,j)
{ for (i = 0; i < 3; i++) 
    { printf "  "
      for (j = 0; j < 3; j++)
        { if (j > 0) { printf " "; }
          printf "%15.8e", M[i,j];
        }
      printf "\n";
    }
}

function r3x3_read(fname,M,  i,j,lin,nfld,fld,nlin)
{ nlin = 0;
  i = 0;
  ERRNO = "0";
  # printf "reading %s ...", fname > "/dev/stderr";
  while((i < 3) && ((getline lin < fname) > 0))
    { nlin++;
      # printf "×" > "/dev/stderr";
      if (! match(lin, /^[ \011]*([#]|$)/))
        { gsub(/[#].*$/, "", lin);
          nfld = split(lin, fld, " ");
          if (nfld != 3) r3x3_read_error(fname, nlin, ("bad line = \"" lin "\""));
          for(j = 0; j < 3; j++) { M[i,j] = fld[j+1]; }
          i++;
        }
  }
  # printf "done \n" > "/dev/stderr";
  if (ERRNO == "No such file or directory") return 0;
  if (ERRNO != "0") { r3x3_read_error(fname, nlin, ERRNO); }
  if (i != 3) { r3x3_error(("file \"" fname "\" is truncated")); }
  return 1;
}

function r3x3_read_error(fname, nlin, msg)
{ printf "file %s, line %d: ** %s\n", fname, nlin, msg > "/dev/stderr";
  abort = 1;
  exit 1
}

function r3x3_error(msg)
{ printf "** %s\n", msg > "/dev/stderr";
  abort = 1;
  exit 1
}