function [L, R, B] = anm_coupling_2D(keff, xm, xp, dely, diff, sigr, chi, nusigf)
% Compute the 2D ANM coupling matrices for a single node. Flat transverse
% leakage approximation.
% function [L, R, B] = anm_coupling_2D(keff, xm, xp, dely, diff, sigr, chi, nusigf)
% input parameters:
% keff    effective multiplication factor
% xm, xp  node support along X-axis
% dely    node size along Y-axis
% diff    diffusion coefficient array (cm)
% sigr    removal cross section array (cm^-1)
% chi     fission spectrum array
% nusigf  nu * fission cross section array (cm^-1)
% output parameters:
% L, R, B nodal coupling matrices
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
  G=size(diff,2) ;
  F=inv(diag(diff))*(chi'*nusigf/keff-sigr) ; [T,Lambda]=eig(F) ;
  delx=xp-xm ;
  Mm=zeros(2*G,2*G) ; Mp=zeros(2*G,2*G) ; Nm=zeros(G,2*G) ; Np=zeros(G,2*G) ;
  B=zeros(G,2*G) ; Y=zeros(G,G) ;
  for g=1:G
    sqla=sqrt(abs(Lambda(g,g))) ;
    if Lambda(g,g) >= 0
      Mm(g,g)=(sin(sqla*xp)-sin(sqla*xm))/(delx*sqla) ;
      Mm(g,G+g)=-(cos(sqla*xp)-cos(sqla*xm))/(delx*sqla) ;
      Mm(G+g,g)=sqla*sin(sqla*xm) ; Mm(G+g,G+g)=-sqla*cos(sqla*xm) ;
      Mp(G+g,g)=sqla*sin(sqla*xp) ; Mp(G+g,G+g)=-sqla*cos(sqla*xp) ;
      Nm(g,g)=cos(sqla*xm) ; Nm(g,G+g)=sin(sqla*xm) ;
      Np(g,g)=cos(sqla*xp) ; Np(g,G+g)=sin(sqla*xp) ;
    else
      Mm(g,g)=(sinh(sqla*xp)-sinh(sqla*xm))/(delx*sqla) ;
      Mm(g,G+g)=(cosh(sqla*xp)-cosh(sqla*xm))/(delx*sqla) ;
      Mm(G+g,g)=-sqla*sinh(sqla*xm) ; Mm(G+g,G+g)=-sqla*cosh(sqla*xm) ;
      Mp(G+g,g)=-sqla*sinh(sqla*xp) ; Mp(G+g,G+g)=-sqla*cosh(sqla*xp) ;
      Nm(g,g)=cosh(sqla*xm) ; Nm(g,G+g)=sinh(sqla*xm) ;
      Np(g,g)=cosh(sqla*xp) ; Np(g,G+g)=sinh(sqla*xp) ;
    end
    Mp(g,g)=Mm(g,g) ; Mp(g,G+g)=Mm(g,G+g) ;
    B(g,g)=-1.0/(delx*Lambda(g,g)) ; B(g,G+g)=1.0/(delx*Lambda(g,g)) ;
    Y(g,g)=1.0/(dely*Lambda(g,g)) ;
  end
  warning off
  GAR1=Nm*inv(Mm) ; GAR2=Np*inv(Mp) ;
  warning on
  MAT1=[GAR1,-Y+GAR1(:,1:G)*Y,Y-GAR1(:,1:G)*Y] ;
  MAT2=[GAR2,-Y+GAR2(:,1:G)*Y,Y-GAR2(:,1:G)*Y] ;
  S=inv(T)*inv(diag(diff)) ;
  L=T*MAT1*blkdiag(inv(T),S,S,S) ; R=T*MAT2*blkdiag(inv(T),S,S,S) ;
  B=T*B*blkdiag(S,S) ;