#! /usr/bin/gawk -f
# Last edited on 2010-10-07 00:38:50 by stolfilocal

# Tries to generate an orthogonal basis for R^n
# by iterated rotation, under the norm 
# <f,g> = (1/n)*(SUM [f[i]**g[i] : i in 0..n-1])

BEGIN{
  if (n == "") { arg_error("must define {n}"); }
  m = n; # For now.
  if (op == "") { arg_error("must define {op}"); }
  split("", U);    # The first basis element.
  split("", R);    # The basis cyclic rotation matrix.
  pi = 3.1415926;
  for (j = 0; j < n; j++) 
    {
      # Clear out the matrix {R}:
      for (i = 0; i < n; i++) { R[i,j] = 0.0; } 
      # Define {U[j]} and column {j} of {R}:
      if (op == "cart") 
        { # Cartesian basis.
            U[j] = (j == 0); 
          i = (j - 1) % n;
          R[i,j] = 1;
        }
      else if (op == "hart")
        { U[j] = 1.00;
          i0 = j; i1 = (n - j) % n;
          # Note that we may have {i0 == i1}:
          R[i0,j] = R[i0,j] + cos(2*pi*j/n);
          R[i1,j] = R[i1,j] + sin(2*pi*j/n);
        }
      else 
        { arg_error(("invalid {op} = \"" op "\"")); }
    }
  
  split("", phi);  # The first {m} elements of the basis.
  
  for (i = 0; i < m; i++)
    { # Generate element {phi[i,0..n-1]} {i}:
      if (i == 0)
        { # Basis element 0 is given: 
          for (j = 0; j < n; j++) { phi[i,j] = U[j]; }
        }
      else
        { # Basis element {i} is row {phi[i-1]} times the matrix {R}: 
          for (j = 0; j < n; j++) 
            { sum = 0;
              for (k = 0; k < n; k++) { sum += phi[i-1,k]*R[k,j]; }
              phi[i,j] = sum;
            }
        }
    }

  # Write out:
  for (j = 0; j < n; j++) 
    { printf "%4d", j;
      for (i = 0; i < m; i++)
        { printf " %+9.6f", phi[i,j];
        }
      printf "\n";
    }

  check_basis(phi,m,n);
}

function check_basis(phi,m,n,  M,i,j,bug,tmp)
  {
    # Checks whether the basis {phi[0..m-1]} with {n} columns is orthonormal.
    bug = 0;
    split("",M);  # Scalar product matrix:
    printf "Scalar product matrix:\n" > "/dev/stderr";
    for (i = 0; i < n; i++)
      { for (j = 0; j < n; j++) 
          { M[i,j] = dot(phi,n,i,j);
            tmp = (i == j);
            if (abs(M[i,j] - tmp) > 1.0e-6) 
              { printf "** <phi[%d],phi[%d]> = %+9.6f\n", M[i,j] > "/dev/stderr";
                bug=1;
              }
          }
      }

    if (bug)
      { 
        printf "Scalar product matrix:\n" > "/dev/stderr";
        for (i = 0; i < n; i++)
          { for (j = 0; j < n; j++) 
              { M[i,j] = dot(phi,n,i,j);
                printf " %+9.6f", M[i,j] > "/dev/stderr";
              }
            printf "\n" > "/dev/stderr";
          }
      }
  }

function dot(phi,n,i,j,   sum,k)
  { 
    # Dot product of rows {i,j} of the matrix {phi} with {n} columns.
    sum = 0;
    for (k = 0; k < n; k++) { sum += phi[i,k]*phi[j,k]; }
    return sum/n;
  }

function abs(x)
  { 
    if (x < 0) { return -x; } else { return x; }
  }

function num_error(i,msg)
  { 
    printf "** element %d: %s\n", i, msg > "/dev/stderr";
    exit 1;
  }

function arg_error(msg)
  { 
    printf "** argument error: %s\n", msg > "/dev/stderr";
    exit 1;
  }