function [fi,LJ1] = degen(J,KJ,LJ1,IN,lab,multab,G,ngen,map,nmberg,cind,norder,steer,inel,ch,sirt)
%
% In case of degeneracy in the eigenvalues, function degen calculates
% a common set of mutual independent commuting group elements and specifies
% a unique common eigencolumn, if possible. 
%
pi2i = 2*pi*i;
kloop = 0;
[fi,nvr1] = eigvec(IN,lab,J,kloop,steer,inel,cind,ch,multab,G);
if nvr1 > 0
  D = fix(nvr1/LJ1);
  dfi(1:nvr1,1:G) = fi(1:G,1:nvr1).';
%
% Calculate a set of mutual independent, commuting group elements to
% element IN.
%
  ntry2 = zeros(1,G);
  I2 = 1;
  ntry2(I2) = 1;
  I2 = multab(I2,IN);
  while I2 ~= 1
    ntry2(I2) = 1;
    I2 = multab(I2,IN);
  end
  ntry1(1) = IN;
  I3 = 1;
  for I2 = 2:G
    if ntry2(I2) ~= 1
      KI = 1;
      while KI <= I3
        if multab(ntry1(KI),I2) ~= multab(I2,ntry1(KI));
          break;
        end
        KI = KI + 1;
      end
      if KI > I3
        ntry3 = zeros(1,G);
        K2 = I2;
        ntry3(K2) = 1;
        K2 = multab(K2,I2);
        while K2 ~= 1
          ntry3(K2) = 1;
          K2 = multab(K2,I2);
        end
        I3 = I3 + 1;
        ntry1(I3) = I2;
        K3 = cind(I2);
        K3 = norder(K3) - 1;
        for K2 = 2:G
          if ntry2(K2) ~= 0
            K4 = 1;
            for K5 = 1:K3
              K4 = multab(K4,I2);
              K6 = multab(K2,K4);
              ntry3(K6) = 1;
            end
          end
        end
        for K2 = 2:G
          if ntry3(K2) == 1
            ntry2(K2) = 1;
          end
        end
      end
    end
  end
  K2 = 1;
  for I2 = 1:G
    if ntry2(I2) == 1
      ntry4(K2) = I2;
      K2 = K2 + 1;
    end
  end
  nvru = K2 - 1;
  if I3 ~= 1
    fid = fopen('Irrepres.m', 'a');
    fprintf(fid, 'Commuting group elements:\n\n');
    fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', ntry1(1:I3));
    fprintf(fid, '\n\n');
    fclose(fid);
%
% A set of commuting group elements is stored in
% ntry1(1:I3). We are looking for common 
% eigenvectors to the set ntry1(1:I3).
%
    I2 = 2;
    while I2 <= I3
      K2 = ntry1(I2);
      K3 = cind(K2);
      K3 = norder(K3);
      for K4 = 1:K3
        eigenv(K4) = exp(pi2i*(K4-1)/K3);
      end
      for K4 = 1:K3
        kloop = 0;
        RIN = eigenv(K4);
        [fi,nvr2] = eigvec(K2,RIN,J,kloop,steer,inel,cind,ch,multab,G);
        kloop = 1;
        if nvr2 ~= 0
          [fi,isec] = intsec(nvr1,nvr2,G,fi,dfi);
          if isec ~= 0
            if isec == LJ1
              break;
            end
            nvr1 = isec;
            dfi(1:nvr1,1:G) = fi(1:G,1:nvr1).';
          end
        end
      end
      if isec == LJ1
        break;
      end
      I2 = I2 + 1;
    end
    if I2 > I3
      fid = fopen('Irrepres.m', 'a');
      fprintf(fid, '\n');
      fprintf(fid, 'Error, not enough basis vectors (only %3d) for irrep %3d', LJ1, J);
      fprintf(fid, '\n\n');
      fclose(fid);
    end
  else
    fid = fopen('Irrepres.m', 'a');
    fprintf(fid, '\n');
    fprintf(fid, 'No commuting element for element %3d', IN);
    fprintf(fid, '\n\n');
    fclose(fid);
  end
  if isec ~= LJ1
    [dfi,nvr1] = subsp(nvr1,ntry4,nvru,J,LJ1,G,dfi,inel,multab,steer,cind,ch,norder);
    if nvr1 ~= LJ1
      dprim = fix(nvr1/LJ1);
      mgt(D,dprim,KJ,LJ1,ntry4,nvru,lab,IN,inel,ngen,map,nmberg,dfi,multab,G,steer);
      LJ1 = 0;
      return
    else
      fi(1:G,1:nvr1) = dfi(1:nvr1,1:G).';
    end
  end
%
% dfi has been reduced to a LJ1 times degenerate space.
%
  ntry1 = zeros(1,G);
  K3 = sirt(J,2);
  K4 = IN;
  for K5 = 1:K3
    K4 = multab(K4,IN);
    ntry1(K4) = 1;
  end
  I2 = 1;
  while I2 <= LJ1
    dfi(1,1:G) = fi(1:G,I2).';
    IV = 1;
    I3 = 1;
    while I3 <= G
      if ntry1(I3) ~= 1
%
% Operate on dfi(1,I) with group element I3
%
        K3 = inel(I3);
        vec(1:G) = dfi(1,multab(K3,1:G));
        I4 = 1;
        while I4 <= IV
          rin = vec(1:G)*dfi(I4,1:G)';
          rinab = abs(rin);
          if abs(rinab - 1) >= 0.000001
            if rinab >= 0.000001
              for K4 = 1:G
                vec(K4) = vec(K4) - rin*dfi(I4,K4);
              end
              rnorm = sum(abs(vec(1:G).').^2);
              if rnorm < 0.000001
                break;
              end
              rnorm2 = 1/sqrt(rnorm);
              vec(1:G) = vec(1:G)*rnorm2;
            end
          else
            break;
          end
          I4 = I4 + 1;
        end
        if I4 > IV
          IV = IV + 1;
          dfi(IV,1:G) = vec(1:G);
          if IV >= LJ1
            fi(1:G,1:LJ1) = dfi(1:LJ1,1:G).';
            return;
          end
        end
      end
      I3 = I3 + 1;
    end
    I2 = I2 + 1;
  end
  if I2 > LJ1
    fid = fopen('Irrepres.m', 'a');
    fprintf(fid, '\n');
    fprintf(fid, 'Error, not enough basisvectors (%3d) for irrep %3d\n', LJ1, J);
    fprintf(fid, '\n');
    fclose(fid);
    LJ1 = 0;
  end
else
  fid = fopen('Irrepres.m', 'a');
  fprintf(fid, '\n');
  fprintf(fid, 'Error, not enough basisvectors (%3d) for irrep %3d\n', LJ1, J);
  fprintf(fid, '\n');
  fclose(fid);
  LJ1 = 0;
end
return