function delta = zhang_fit(berg,R,f)
%ZHANG_FIT - Compute 'Berg' parameters for EEG multilayer sphere
% function delta = zhang_fit(berg,R,f)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% MINIMIZATION FUNCTION FOR COMPUTING BERG PARAMETERS FOR EEG MULTILAYER 
% SPHERICAL FORWARD MODEL (zhang_fit.m)
%
% This function serves as a minimization function for determining the Berg 
% Eccentricity and Magnitude parameters associated with a series approximation
% of a Single Dipole in a Multilayer Sphere -by- multiple dipoles in a single shell. 
% (Zhang: Eq 5i"). This function is intended to be called by the Matlab "fmins" 
% function. This version has not been optimized for speed.
%
% (Ref: Z. Zhang "A fast method to compute surface potentials generated by dipoles 
% within multilayer anisotropic spheres" (Phys. Med. Biol. 40, pp335-349,1995)   
% 
% INPUTS:
%         berg : Contains                                              (2*J-1) x 1
%         R    : Concentric Sphere Radii(in meters) from INNERMOST to
%                OUTERMOST                                                  NL x 1
%         f: Weighting terms corresponding to sphere/conductivity 
%            profile of interest                                          nmax x 1 
%  
%      where: NL = # of sphere layers; nmax = # of series terms; J = # of Berg Dipoles
%
% OUTPUTS:
%          delta: Minimization function output. Best fit of 
%                 Berg Eccentricity and Magnitude Paramters occurs when
%                 this value is minimized                                    1 x 1
%
% External Functions and Files: None
%    
% By John Ermer 5/5/99 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%<autobegin> ---------------------- 26-May-2004 11:34:43 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% At Check-in: $Author: Mosher $  $Revision: 13 $  $Date: 5/26/04 10:02a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:34:43 -----------------------


% History ------------------------------------------------
% Ermer 5-May-1999 creation
% JCM 08-Sep-2003 Comments
% ---------------------------------------------------------


NL = length(R);           % Number of Sphere Layers
nmax = length(f);         % Number of Legendre Series Terms
J = (length(berg)+1)/2;   % Infer Number of Berg Dipoles from Passed data
mu_berg = berg(1:J);      % Berg Eccentricity Factors 1->J
lam_berg = [0.0 berg(J+1:2*J-1)]; % Berg Magnitude Factors 2->J
%
delta = 0.0;  % Compute Factor to be minimized
%
for n = 2:nmax
   term0 = 0;
   for j=2:J
      term0 = term0 + lam_berg(j)*(mu_berg(j)^(n-1) - mu_berg(1)^(n-1));
   end
   term1 = ( (R(1)/R(NL))^(n-1) )* ...
     ( f(n)-f(1)*mu_berg(1)^(n-1) - term0);
   delta = delta + term1*term1;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
