vb_multi_shell_param

 Legendre Expansion Weights for EEG multilayer spherical forward model 
 f = vb_multi_shell_param(R, sigma, Nmax)
 
 INPUTS
    R    : Radii of sphere from innermost to outermost    ( NL x 1 )
   sigma : Conductivity    from innermost to outermost    ( NL x 1 )
   Nmax  : # of Legendre Expansion     (Default: 200)      ( scalar )
  
     NL = # of sphere layers; 

 OUTPUTS:
     f: Legendre Expansion Weights for multishell model  ( Nmax x 1 )


 Single Dipole in a multilayer sphere can be approximated 
        by multiple dipoles in a single shell. 
 Z. Zhang,  Phys. Med. Biol., vol.40, pp. 335-349, 1995.

 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL.46, 245-259, 1999
 John C. Mosher, Richard M. Leahy and Paul S. Lewis

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function f = vb_multi_shell_param(R, sigma, Nmax)
0002 % Legendre Expansion Weights for EEG multilayer spherical forward model
0003 % f = vb_multi_shell_param(R, sigma, Nmax)
0004 %
0005 % INPUTS
0006 %    R    : Radii of sphere from innermost to outermost    ( NL x 1 )
0007 %   sigma : Conductivity    from innermost to outermost    ( NL x 1 )
0008 %   Nmax  : # of Legendre Expansion     (Default: 200)      ( scalar )
0009 %
0010 %     NL = # of sphere layers;
0011 %
0012 % OUTPUTS:
0013 %     f: Legendre Expansion Weights for multishell model  ( Nmax x 1 )
0014 %
0015 %
0016 % Single Dipole in a multilayer sphere can be approximated
0017 %        by multiple dipoles in a single shell.
0018 % Z. Zhang,  Phys. Med. Biol., vol.40, pp. 335-349, 1995.
0019 %
0020 % IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL.46, 245-259, 1999
0021 % John C. Mosher, Richard M. Leahy and Paul S. Lewis
0022 %
0023 % Copyright (C) 2011, ATR All Rights Reserved.
0024 % License : New BSD License(see VBMEG_LICENSE.txt)
0025 
0026 % # of sphere layer
0027 NL = length(R);
0028 
0029 if nargin < 3    
0030     % Default # of Legendre expansion
0031     Nmax = 200;  
0032 end
0033 
0034 % Radius of outermost layer (Sensor distance from origin)
0035 Rmax = R(NL);     
0036 
0037 %
0038 % Computes the weights f(n) for a legendre expansion
0039 % f(n) depend only on the multisphere radii and conductivity
0040 
0041 if NL==1
0042    f=ones(Nmax,1);
0043 else
0044     
0045     for k = 1:NL-1 
0046        s(k) = sigma(k)/sigma(k+1);
0047     end
0048     
0049     a    = Rmax./R;
0050     ainv = R/Rmax;
0051     sm1  = s-1;
0052     npow = 2*[1:Nmax]+1;
0053     npow = npow(:);
0054     
0055     f = zeros(Nmax,1);
0056     
0057     for n = 1:Nmax
0058         np1 = n+1;
0059         Mc = eye(2);
0060 
0061         %  Eq.(17)
0062         for k = 2:NL-1
0063            Mc = Mc*[n+np1*s(k),  np1*sm1(k)*a(k)^npow(n);...
0064                     n*sm1(k)*ainv(k)^npow(n) , np1+n*s(k)];
0065         end
0066 
0067         % Eq.(16)
0068         Mc(2,:) = [n*sm1(1)*ainv(1)^npow(n) , np1+n*s(1)]*Mc; 
0069         Mc = Mc/(npow(n))^(NL-1);
0070         f(n) = n/( n*Mc(2,2) + np1*Mc(2,1) );
0071     end
0072 end
0073