private void solveTridiag(float sub[], float diag[], float sup[], float b[], int n){
/*                  solve linear system with tridiagonal n by n matrix a
                    using Gaussian elimination *without* pivoting
                    where   a(i,i-1) = sub[i]  for 2<=i<=n
                            a(i,i)   = diag[i] for 1<=i<=n
                            a(i,i+1) = sup[i]  for 1<=i<=n-1
                    (the values sub[1], sup[n] are ignored)
                    right hand side vector b[1:n] is overwritten with solution 
                    NOTE: 1...n is used in all arrays, 0 is unused */
      int i;
/*                  factorization and forward substitution */
      for(i=2; i<=n; i++){
        sub[i] = sub[i]/diag[i-1];
        diag[i] = diag[i] - sub[i]*sup[i-1];
        b[i] = b[i] - sub[i]*b[i-1];
      }
      b[n] = b[n]/diag[n];
      for(i=n-1;i>=1;i--){
        b[i] = (b[i] - sup[i]*b[i+1])/diag[i];
      }
    }