function psi=eigenfunction_box(l,L)
%
% psi=eigenfunction_box(l,L)
%
% The function computes gives the eigenfunctions of the Laplace-Beltrami
% operator in a 2D or 3D box. If you are not worrying about the boundary
% value problem (continuity of the dervivative at the boundary), sine basis is sufficient.
% If you are worrying about the boudary, email Moo K. Chung for your
% problem. 
%
%     l         :  degrees of basis along each dimension. l=0 gives 1/sqrt(area) or 1/sqrt(volume) 
%                  as the first basis. 
%    L         :  dimension of box. Ex. L= [10 20] constructs the eigenfunctions of degree l in 10 x 20 image 
%                  L=[10 20 30] constructs the eigenfunctions of degree in 10 x 20 x 30 image 
%
%
%
%
%
%     (C) Moo K. Chung 2008
%     Department of Biostatistics and Medical Informatics
%     Waisman Laboratory for Brain Imaging
%     University of Wisconsin-Maison
%  
%     hptt://www.stat.wisc.edu/~mchung
%     email://mkchung@wisc.edu
%
%  Histroy: Febuary  8, 2008 created; November 22, 2011 documentation added

% This is not optimized properly.

if length(L)==2
    l1=l(1);
    l2=l(2);
    L1=L(1);
    L2=L(2);
    
    
    x1=1:L1;
    %psi1=2*sqrt(2)*sin (pi*l1*x1/L1)';
    psi1=2*sin (pi*l1*x1/L1)';
        
    psi1=repmat(psi1,[1,L2]);
    %figure; imagesc(psi1); colorbar
    
    x2=1:L2;
    %psi2=sin (pi*l2*x2/L2);
    psi2=sin (pi*l2*x2/L2);
    psi2=repmat(psi2,[L1,1]);
    %figure; imagesc(psi2); colorbar
    
    psi=psi1.*psi2;
    %figure; imagesc(psi); colorbar
end;

if length(L)==3
    l1=l(1);
    l2=l(2);
    l3=l(3);
    L1=L(1);
    L2=L(2);
    L3=L(3);
    
    x1=1:L1;
    psi1=2*sqrt(2)*sin (pi*l1*x1/L1)';
    psi1=repmat(psi1,[1,L2,L3]);
    %figure; imagesc(squeeze(psi1(:,50,:)))
    
    x2=1:L2;
    psi2=sin (pi*l2*x2/L2);
    psi2=repmat(psi2,[L1,1,L3]);
    
    x3=1:L3;
    psi3=sin (pi*l3*x3/L3);
    psi3=repmat(psi3,[L2, 1]);
    %imagesc(psi3)
    psi3=reshape(psi3,1,L2*L3);
    psi3=repmat(psi3,[L1,1]);
    psi3=reshape(psi3,L1,L2,L3);
    %figure; imagesc(squeeze(psi3(50,:,:)))
    psi=psi1.*psi2.*psi3;
end;




