root/DPWP/trunk/DPWP/diwasp_1_1GD/private/IMLM.m

function [S]=IMLM(xps,trm,kx,Ss,pidirs,miter,displ)

%this code was edited significantly on 5/2010 --- see comments for details

gamma=.1;
beta=1.0;
alpha=0.1;

szd=size(xps,1);
ffreqs=size(xps,3);
ddirs=size(trm,3);

ddir=8*atan(1.0)/ddirs;

if(displ<2)
   warning off;
end


for ff=1:ffreqs
   if(displ>=1)
      disp(['calculating for frequency' blanks(1) num2str(ff) ' of' blanks(1) num2str(ffreqs)]);
  end
   
   for m=1:szd
      for n=1:szd
        H(1:ddirs,m,n)=trm(n,ff,1:ddirs);
      Hs(1:ddirs,m,n)=conj(trm(m,ff,1:ddirs));
      expx(1:ddirs,m,n)=exp(i*kx(m,n,ff,1:ddirs));
      iexpx(1:ddirs,m,n)=exp(-i*kx(m,n,ff,1:ddirs));
      Htemp(:,m,n)=H(:,m,n).*Hs(:,m,n).*expx(:,m,n);
      iHtemp(:,m,n)=H(:,m,n).*Hs(:,m,n).*iexpx(:,m,n);
   end
   end
   
   I=eye(m)*sqrt(eps);
   invcps=pinv(xps(:,:,ff)+I); %changed on 5/5/10 to pinv and +I
   Sftmp=zeros(ddirs,1);
   for m=1:szd
      for n=1:szd
                xtemp=invcps(m,n)*Htemp(:,m,n);
                        Sftmp(:)=Sftmp(:)+xtemp;
                end
   end
   Eo=(1./real(Sftmp(:))); %added 5/5/10 to take just the real part
   [Eo]=normalize(Eo,ddir); %added 5/11/10 to remove negative values before normalizing
   E=Eo;
   T=Eo;

   for it=1:miter   
      for m=1:szd
        for n=1:szd
                expG(m,n,:)=iHtemp(:,m,n).*E(:);   
                        ixps(m,n)=sum(expG(m,n,:))*ddir;
                end
      end
      invcps=pinv(ixps+I); % changed on 5/5/10 to pinv and +I 
      Sftmp=zeros(ddirs,1);
      for m=1:szd
        for n=1:szd
                xtemp=invcps(m,n)*Htemp(:,m,n);
                        Sftmp(:)=Sftmp(:)+xtemp;
                        end
        end
        Told=T;
      T=(1./real(Sftmp(:))); %added on 5/5/10 to take just real part
      [T]=normalize(T,ddir); %added 5/11/10 to remove negative values before normalizing 
      
      %lambda=ones(size(T))-(T./Eo)
      %ei=gamma*lambda.*E;
      ei=gamma*((Eo-T)+alpha*(T-Told));
      %lambda=Eo-T;      %edited on 5/5/10 to calculate ei slightly differently 
      %ei=(abs(lambda).^(beta+1))./(gamma.*lambda);
      E=E+ei;
      [E]=normalize(E,ddir);%added 5/11/10 to remove negative values before normalizing

            
   end
         
   S(ff,:)=Ss(1,ff)*E';
end

warning on;