function [fi,nvec] = eigvec(elem,eval,J,kloop,steer,inel,cind,ch,multab,G)
%
% This function calculates the eigenvectors of group element
% elem, corresponding to eigenvalue eval and projects these 
% eigenvectors on the Jth irreducible subspace. The eigenvectors
% are stored in fi(1:G,1:nvec). 
%
if kloop ~= 1
  [loopl,lpstr,numl] = permu(multab,G,steer,inel,elem);
end
IND = 1;
nvr = 0;
for nml = 1:numl
  if nml > 1
    IND = IND + loopl(nml-1);
  end
  nvr = nvr + 1;
  vec = zeros(1,G);
  K2 = lpstr(IND);
  vec(K2) = 1;
  P = 1;
  P = P*eval;
  LPE = IND + loopl(nml);
  for N = (IND+1):(LPE-1)
    K2 = lpstr(N);
    vec(K2) = P;
    P = P*eval;
  end
%
% Projection SJ*vec:
%
  for K4 = 1:G
    fi(K4,nvr) = 0;
  end
  for K4 = 1:G
    for K5 = 1:G
      K6 = inel(K5);
      K6 = multab(K4,K6);
      K6 = cind(K6);
      fi(K4,nvr) = fi(K4,nvr) + conj(ch(J,K6))*vec(K5);
    end
  end
end
%
% The first non-zero eigencolumn is normalized.
%
nvec = nvr;
nvr = 1;
while nvr <= nvec
%
% In case there are no eigenvectors with eigenvalue eval, belonging to
% the Jth irreducible subspace, set nvec = 0 and return
%
  II = 1;
  while II <= G
    if abs(fi(II,nvr)) >= 0.001
      break;
    end
    II = II + 1;
  end
  if II <= G
    rnorm = sum(abs(fi(1:G,nvr)).^2);
    rnorm = 1/sqrt(rnorm);
    fi(1:G,nvr) = fi(1:G,nvr)*rnorm;
    if nvr > 1
%
% We orthogonalize the eigencolumns of element elem with eigenvalue eval,
% to each other.
%
      K3 = nvr - 1;
      vec(1:K3) = fi(1:G,1:K3)'*fi(1:G,nvr); 
      for I3 = 1:G
        vec2(I3) = fi(I3,nvr);
        for K2 = 1:K3
          vec2(I3) = vec2(I3) - vec(K2)*fi(I3,K2);
        end
      end
      I3 = 1;
      while I3 <= G
        if abs(vec2(I3)) >= 0.001
          break;
        end
        I3 = I3 + 1;
      end
      if I3 <= G
        rnorm = sum(abs(vec2(1:G).').^2);
        rnorm = 1/sqrt(rnorm);
        fi(1:G,nvr) = vec2(1:G).'*rnorm;
        nvr = nvr + 1;
      else
        nvec = nvec - 1;
        fi(1:G,nvr:nvec) = fi(1:G,(nvr+1):(nvec+1));
      end
    else
      nvr = nvr + 1;
    end
  else
    nvec = nvec - 1;
    fi(1:G,nvr:nvec) = fi(1:G,(nvr+1):(nvec+1));
  end
end
%
% The nvec orthonormal eigencolumns to element elem with eigenvalue eval
% are stored in fi(1:G,1:nvec).
%
return  