docs: private double[][] determineBasis(double[] alpha) {
      final int dim = alpha.length;
      // Primary vector:
      double[] nn = new double[dim + 1];
      for(int i = 0; i < nn.length; i++) {
        double alpha_i = i == alpha.length ? 0 : alpha[i];
        nn[i] = ParameterizationFunction.sinusProduct(0, i, alpha) * FastMath.cos(alpha_i);
      }
      timesEquals(nn, 1. / euclideanLength(nn)); // Normalize
      // Find orthogonal system, in transposed form:
      double[][] basis = new double[dim][];
      int found = 0;
      for(int i = 0; i < nn.length && found < dim; i++) {
        // ith unit vector.
        final double[] e_i = new double[nn.length];
        e_i[i] = 1.0;
        minusTimesEquals(e_i, nn, scalarProduct(e_i, nn));
        double len = euclideanLength(e_i);
        // Make orthogonal to earlier (normal) basis vectors:
        for(int j = 0; j < found; j++) {
          if(len < 1e-9) { // Disappeared, probably linear dependent
            break;
          }
          minusTimesEquals(e_i, basis[j], scalarProduct(e_i, basis[j]));
          len = euclideanLength(e_i);
        }
        if(len < 1e-9) {
          continue;
        }
        timesEquals(e_i, 1. / len); // Normalize
        basis[found++] = e_i;
      }
      if(found < dim) {
        // Likely some numerical instability, should not happen.
        for(int i = found; i < dim; i++) {
          basis[i] = new double[nn.length]; // Append zero vectors
        }
      }
      return transpose(basis);
    }