interch.m

Go to the documentation of this file.

function [OUTS,OUToldS,OUTD]=interch(S,oldS,D,Ncond,fstep);    %03.04.2001, Bjorn G.
%CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
%C Interchanging e-vectors and corresponding e-values so that they become smooth
%C functions of frequency. Considers the old eigenvectors one-by-one and finds out by which of 
%C the new ones the largest dot product is attained. 
 
% Translation: stoerst=largest, hjelp=help, opptatt=taken 
 
if fstep>1 
  for ii=1:Ncond 
    istoerst=0;
    rstoerst=0;
    for j=1:Ncond
      dotprod=0;
      for k=1:Ncond
        dotprod=dotprod+conj(oldS(k,ii))*S(k,j);
      end
      if abs(real(dotprod))>rstoerst 
        rstoerst=abs(real(dotprod));
        istoerst=j;       
      end
    end
    dot(ii)=rstoerst;
    ind(ii)=ii;       
    opptatt(ii)=0;                                       
  end %
 
%C Sorting inner products abs(realdel) in descending order :
  for ii=1:Ncond  
    for j=1:Ncond-1
      if(dot(j)<dot(j+1));
        hjelp(1)=dot(j+1);
        ihjelp  =ind(j+1);
        dot(j+1)=dot(j);
        ind(j+1)=ind(j);
        dot(j)=hjelp(1);
        ind(j)=ihjelp;
      end
    end
  end
 
%C Doing the interchange in a prioritized sequence :
  for l=1:Ncond %
    ii=ind(l);
    istoerst=0;
    rstoerst=0;
 
    for j=1:Ncond 
      if opptatt(j)==0   
        dotprod=0;
        for k=1:Ncond
          dotprod=dotprod+conj(oldS(k,ii))*S(k,j);
        end                        
        if abs(real(dotprod))>rstoerst
          rstoerst=abs(real(dotprod));
          istoerst=j;       
        end
      end 
    end 
 
    opptatt(ii)=1;
 
     hjelp=S(:,ii);
     S(:,ii)=S(:,istoerst);
     S(:,istoerst)=hjelp;
 
     hjelp=D(ii,ii);
     D(ii,ii)=D(istoerst,istoerst);
     D(istoerst,istoerst)=hjelp;
  end 
 
%C Finding out whether the direction of e-vectors are 180deg. in error. Is done by comparing
%C sign of dotproducts of e-vectors, for new and old S-matrix :
  for ii=1:Ncond
    dotprod=0;
    for k=1:Ncond
      dotprod=dotprod+conj(oldS(k,ii))*S(k,ii);            
    end
    if sign(real(dotprod))<0 
      S(:,ii)=-S(:,ii);
    end
  end
 
end %
 
oldS=S;
 
 
OUTS=S; OUToldS=S; OUTD=D;