[Master Index] [Index for Toolbox]

# eeg_sph

## (Toolbox/eeg_sph.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);

### Help Text

    EEG_SPH - Calculate the electric potential , spherical head, arbitrary orientation
     function G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
     function G = eeg_sph(L,Channel,Param,Order);
     L is 3 x nL, each column a source location
     Channel is the channel structure, same for Param
     Order is
      -1 current dipole
       0 focal(magnetic) dipole % NOT SUPPORTED
       1 1st order multipole    % NOT SUPPORTED
     Verbose : toggle Verbose mode

### Cross-Reference Information

This function calls

    dlegpoly                           C:\BrainStorm_2001\Toolbox\dlegpoly.m
    dotprod                            C:\BrainStorm_2001\Toolbox\dotprod.m
    good_channel                       C:\BrainStorm_2001\Toolbox\good_channel.m
    rownorm                            C:\BrainStorm_2001\Toolbox\rownorm.m

This function is called by

    bst_headmodeler                    C:\BrainStorm_2001\Toolbox\bst_headmodeler.m

### Listing of function C:\BrainStorm_2001\Toolbox\eeg_sph.m

    function G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
    %EEG_SPH - Calculate the electric potential , spherical head, arbitrary orientation
    % function G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
    % function G = eeg_sph(L,Channel,Param,Order);
    % L is 3 x nL, each column a source location
    % Channel is the channel structure, same for Param
    % Order is
    %  -1 current dipole
    %   0 focal(magnetic) dipole % NOT SUPPORTED
    %   1 1st order multipole    % NOT SUPPORTED
    % Verbose : toggle Verbose mode

    %<autobegin> ---------------------- 14-Jun-2004 17:10:26 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Forward Modeling
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\dlegpoly.m
    %   toolbox\dotprod.m
    %   toolbox\good_channel.m
    %   toolbox\rownorm.m
    %
    % Subfunctions in this file, in order of occurrence in file:
    %   G = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in,lam_berg_in)
    %
    % At Check-in: $Author: Mosher $  $Revision: 19 $  $Date: 6/14/04 3:37p $
    %
    % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-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> ------------------------ 14-Jun-2004 17:10:26 -----------------------

    % /---Script Authors--------------------------------------\
    % |                                                       |
    % | *** John Ermer, Ph.D.                                 |
    % |    Signal % Image Processing Institute                |
    % |    University of Southern California                  |
    % |    Los Angeles, CA, USA                               |
    % |                                                       |
    % | *** John C. Mosher, Ph.D.                             |
    % |  Biophysics Group                                    |
    % |                                                       |
    % | *** Sylvain Baillet Ph.D.                             |
    % | Cognitive Neuroscience & Brain Imaging Laboratory     |
    % | CNRS UPR640 - LENA                                    |
    % | Hopital de la Salpetriere, Paris, France              |
    % | sylvain.baillet@chups.jussieu.fr                      |
    % |                                                       |
    % \-------------------------------------------------------/
    %
    % Date of creation: October, 25 1999
    %
    % Script History -----------------------------------------------------------------------------------------------------------
    % SB  19-Nov-2002 : Edited Header
    %                  Updated management of EEG reference
    % JCM 20-Nov-2002 : Fixed headers to have only one autocomments block
    % SB  09-Mar-2004 : Added verbose mode
    % --------------------------------------------------------------------------------------------------------------------------

    if nargin < 5
        Verbose = 1; % Default
    end

    %-----------------------------

    nmax = 80;

    %-----------------------------

    % EEG Channels
    EEGndx = good_channel(Channel,[],'EEG');

    % Reference Channel
    REFndx = good_channel(Channel,[],'EEG REF');

    % EEG Channel locations

    % Add EEG reference at the end of Channel and Param structures
    Re = [Channel(EEGndx).Loc,Channel(REFndx).Loc]'; % Electrode location array
    if length(Param) ~= length([EEGndx,REFndx])
        Param(REFndx) = Param(EEGndx(1));
    end

    Param = Param([EEGndx,REFndx]);

    Rq = L';
    center = [Param.Center]';
    Re = Re - center;

    Rq = Rq - repmat(center(1,:),size(Rq,1),1); % Back to origin [0 0 0] for the sensors and the sources

    clear tmp

    switch(Param(1).EEGType)
    case 'EEG_SINGLE'
        R = Param(1).Radii(end);
        sigma = Param(1).Conductivity(end);
    otherwise
        R = Param(1).Radii;
        sigma = Param(1).Conductivity;
    end

    switch(Param(1).EEGType)

    case 'EEG_BERG'
        method = 2;
        mu_berg_in = Param(1).Berg.mu;
        lam_berg_in  = Param(1).Berg.lam;
    case 'EEG_3SHELL'
        method = 1;
        mu_berg_in = [];
        lam_berg_in = [];
    otherwise
        method = 1;
        mu_berg_in = [];
        lam_berg_in  = [];
    end

    if 1 % Projection of the EEG sensors on the sphere
        [theta phi Re_sph] = cart2sph(Re(:,1),Re(:,2),Re(:,3));
        Re_sph = R(end)*ones(size(Re_sph));
        [Re(:,1) Re(:,2) Re(:,3)] = sph2cart(theta,phi,Re_sph);

    end

    Gtmp = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in',lam_berg_in');

    if isempty(REFndx) % Average Reference

        G = NaN * zeros(length(EEGndx),size(Gtmp,2));
        %G(EEGndx,:) = Gtmp - repmat(mean(Gtmp),size(Gtmp,1),1);
        G = Gtmp - repmat(mean(Gtmp),size(Gtmp,1),1);
        clear Gtmp

    else % Specific electrode as reference
        % Remove its lead field from all others
        % Note: reference leadfield is at the end of original G matrix (from gain_sph)
        G = NaN * zeros(length(EEGndx),size(Gtmp,2));
        %G(EEGndx,:) = Gtmp(1:end-1,:) - repmat(Gtmp(end,:),size(Gtmp,1)-1,1);
        %G = G(EEGndx,:);
        G = Gtmp(1:end-1,:) - repmat(Gtmp(end,:),size(Gtmp,1)-1,1);

        %     % Reference lead field is zero
        %     G(REFndx,:) = zeros(1,size(G,2));
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Subfunctions -------------------------------------------------------------------------------
    %
    %
    %
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function G = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in,lam_berg_in)

    %GAINP_SPH6X EEG Multilayer Spherical Forward Model
    % function G = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in,lam_berg_in)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % EEG MULTILAYER SPHERICAL FORWARD MODEL (gainp_sph6x.m)
    %
    % This function computes the voltage potential forward gain matrix for an array of
    % EEG electrodes on the outermost layer of a single/multilayer conductive sphere.
    % Each region of the multilayer sphere is assumed to be concentric with
    % isontropic conductivity.  EEG sensors are assumed to be located on the surface
    % of the outermost sphere.
    %
    % Calculation of the electric potiential is performed using either (user-specified)
    % of the following methods (Ref: Z. Zhang "A fast method to compute surface
    % potentials generated by dipoles within multilayer anisotropic spheres"
    % (Phys. Med. Biol. 40, pp335-349,1995)
    %
    %   1) Closed Form Solution (Single Shell Case Only). See formulas (1H,1H')
    %
    %   2) Series Expansion using Legendre Polynomials. See formulas (1I,2I,3I and 4I)
    %
    %   3) Series Approximiation of a Multilayer Sphere as three dipoles in a
    %      single shell using "Berg/Sherg" parameter approximation.
    %      See formulas (1i',5i" and 6i)
    %
    % Dipole generator(s) are assumed to be interior to the innermost "core" layer. For those
    % dipoles external to the sphere, the dipole "image" is computed and used determine the
    % gain function. The exception to this is the Legendre Method where all dipoles MUST be
    % interior to the innermost "core" layer.
    %
    % INPUTS (Required):
    %       Rq   : dipole location(in meters)                                 P x 3
    %       Re   : EEG sensors(in meters) on the scalp                        M x 3
    %       R    : radii(in meters) of sphere from
    %              INNERMOST to OUTERMOST                                     NL x 1
    %       sigma: conductivity from INNERMOST to OUTERMOST                   NL x 1
    %
    % INPUTS (Optional):
    %       nmax : # of terms used in Truncated Legendre Series               scalar
    %              If not specified, a default value based on outermost
    %              dipole magnitude is computed. (Note: This parameter
    %              is ignored when Berg Approximation is commanded)
    %     method : Method used for computing forward potential
    %              1=Legendre Series Approx; 2=Berg Parameter Approx
    %              (Note: Default and all other values invoke Legendre
    %              Series Approx. Exception is single-shell case where
    %              closed form solution is always used)                       scalar
    %  mu_berg_in: User specified initial value for Berg eccentricity
    %              factors (Required if Berg Method is commanded)             3 x 1
    % lam_berg_in: User specified initial value for Berg magnitude
    %              factors (Required if Berg Method is commanded)             3 x 1
    %
    %              WHERE: M=# of sensors; P=# of dipoles; NL = # of sphere layers
    %
    % OUTPUTS:
    %       G    : EEG forward model gain matrix                              M x (3*P)
    %
    % External Functions and Files:
    %    dlegpoly.m; rownorm.m; dotprod.m: USC/LANL MEG/EEG Toolbox
    %    zhang_fit.m: External Function used to fit Berg Parameters (Zhang Eq# 5i")
    %
    %     - John Ermer 6/3/99
    %        - 8/9/99: Modified to use EM image for dipoles external to brain
    %                  (Applies to single-shell and Berg Methods only) (John Ermer)
    %        - 10/31/99: Optimized Processing Associated with Dipoles falling outside sphere
    %                  (Applies to single-shell and Berg Methods only) (John Ermer)
    %        - 01/26/00: Corrected minor dimension error which caused program to fault when
    %                  external dipoles and multiple sensors were present. (John Ermer)
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%% THIS PART CHECKS INPUT PARAMETERS FOR THEIR DIMENSION AND VALIDITY %%%
    %
    NL = length(R);                          % # of concentric sphere layers
    P = size(Rq,1);
    M = size(Re,1);
    %
    if R(1)~= min(R)
        error('Head radii must be specified from innermost to outmost layer!!! ')
    end
    %
    if size(Rq,2) ~= 3
        error('Dipole location must have three columns!!!')
    end
    %
    if nargin < 6,     % Check # of input terms to see if method is specified
        method = 1;                        % Default Method = Legendre Series Expansion
    else
        if (method>2)|(method<0)
            method = 0;                     % Default Method = Legendre Series Expansion
        end
    end
    %
    %%% This part pre-initializes parameters used in future calculations %%%
    %
    G = zeros(M,3*P);   % Pre-Allocate Gain Matrix
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% This part computes the potential for a dipole contained within a single-layer
    %%% homogeneous sphere using closed form formula
    %%% The EEG single-shell solution uses the vector form (which avoids computationally
    %%% expensive intrinsic functions) described by Mosher et al ("EEG and MEG: Forward
    %%% Solutions for inverse problems" IEEE BME Trans March 1999)
    %
    if NL == 1  % Single Shell Case (Closed Form Solution)
        %
        Re_mag = repmat(R(NL),P,M);           %(PxM)
        Re_mag_sq = repmat(R(NL)*R(NL),P,M);  %(PxM)
        Rq_mag = rownorm(Rq);                 %(Px1)
        %
        Rq1 = Rq;                             %(Px1)
        Rq1_mag = Rq_mag;                     %(Px1)
        Rq1_mag_sq = Rq_mag.*Rq_mag;          %(Px1)
        Re_dot_Rq1 = Rq1*Re';                 %(PxM)
        %
        const = 4.0*pi*sigma(NL);
        term = 1./(const*Rq1_mag_sq);         %(Px1)
        %
        %%% This part checks for the presence of Berg dipoles which are external to
        %%% the sphere. For those dipoles external to the sphere, the dipole parameters
        %%% are replaced with the electrical image (internal to sphere) of the dipole
        %
        nx = find(Rq1_mag > R(NL));
        %
        if nx>0
            Rq1_temp = Rq1(nx,:);
            Rq1(nx,:) = R(NL)*R(NL)*Rq1_temp./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,3);
            Rq1_mag(nx,1) = rownorm(Rq1(nx,:));
            Rq1_mag_sq(nx,1) = Rq1_mag(nx,1).*Rq1_mag(nx,1);
            Re_dot_Rq1(nx,:) = R(NL)*R(NL)*Re_dot_Rq1(nx,:)./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,M);
            term(nx,:) = (R(NL)./rownorm(Rq1_temp)).*term(nx,:);
            %  was : term(nx,:) = (R(NL)/rownorm(Rq1_temp))'.*term(nx,:);

        end
        %
        %%% Calculation of Forward Gain Matrix Contribution due to K-th Berg Dipole
        %
        Rq1_mag = repmat(Rq1_mag,1,M);                          %(PxM)
        Rq1_mag_sq = repmat(Rq1_mag_sq,1,M);                    %(PxM)
        term = repmat(term,1,M);                                %(PxM)
        %
        d_mag = reshape( rownorm(reshape(repmat(Re,1,P)',3,P*M)' ...
            -repmat(Rq1,M,1)) ,P,M);           %(PxM)
        d_mag_cub = d_mag.*d_mag.*d_mag;                        %(PxM)
        F_scalar = d_mag.*(Re_mag.*d_mag+Re_mag_sq-Re_dot_Rq1); %(PxM)
        c1 = term.*(2*( (Re_dot_Rq1-Rq1_mag_sq)./d_mag_cub) ...
            + 1./d_mag - 1./Re_mag);                 %(PxM)
        c2 = term.*((2./d_mag_cub) + (d_mag+Re_mag)./(Re_mag.*F_scalar));
        %
        G = G + reshape(repmat((c1 - c2.*Re_dot_Rq1)',3,1),M,3*P) ...
            .*repmat(reshape(Rq1',1,3*P),M,1) ...
            +  reshape(repmat((c2.*Rq1_mag_sq)',3,1),M,3*P) ...
            .*repmat(Re,1,P);
        %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        %%% This part computes the potential for a dipole contained within a multi-layer
        %%% isontropic sphere using Legendre Polynomial Expansion (Zhang Eqs 1H, 1H')
        %%% (Code based on gainp_sph.m by CCH,  Aug/20/1995)
        %
    elseif (NL>1)&(method==1)  % Multi-shell Case Using Legendre Expansion
        %
        Rq_mag = rownorm(Rq);
        %
        if ~all( Rq_mag < R(1)+eps )    % check if dipoles within the brain
            errordlg('Legendre method assumes all dipoles(s) inside brain layer - please modify dipole locations OR use the "Berg" approach')
            return
        end
        %
        %  compute weights fn. fn depends only on the radii and cdv
        Rq_mag = rownorm(Rq);
        Re_mag = R(NL);         % Radius of outermost layer (Sensor distance from origin
        %
        %
        if nargin < 5,     % check # of inputs to see if nmax was specified
            nmax = fix(10/(1-max(Rq_mag)/Re_mag));    % Default for # Legendre Series Terms
        end
        %
        Ren = Re/Re_mag;
        Rqn = Rq./[Rq_mag,Rq_mag,Rq_mag];
        for k = 1:NL-1
            s(k) = sigma(k)/sigma(k+1);
        end
        a = Re_mag./R;
        ainv = R/Re_mag;
        sm1 = s-1;
        twonp1 = 2*[1:nmax]+1;
        twonp1 = twonp1(:);
        f = zeros(nmax,1);
        %
        for n = 1:nmax
            np1 = n+1;
            Mc = eye(2);
            for k = 2:NL-1,
                Mc = Mc*[n+np1*s(k),  np1*sm1(k)*a(k)^twonp1(n);...
                        n*sm1(k)*ainv(k)^twonp1(n) , np1+n*s(k)];
            end;
            Mc(2,:) = [n*sm1(1)*ainv(1)^twonp1(n) , np1+n*s(1)]*Mc;
            Mc = Mc/(twonp1(n))^(NL-1);
            f(n) = n/(n*Mc(2,2)+np1*Mc(2,1));
        end;
        %
        onevec = ones(M,1);
        wtemp = ((twonp1./[1:nmax]').*f)/(4*pi*sigma(NL)*R(NL)^2);
        n = [1:nmax]';
        nm1 = n-1;

        for i = 1:P,  % loop over all dipoles
            rqn = [Rqn(i,1)*onevec,Rqn(i,2)*onevec,Rqn(i,3)*onevec];
            cosgamma = dotprod(rqn,Ren);
            rqn = rqn(1,:);
            [Pl,dP] = dlegpoly(nmax,cosgamma);  % evaluate legendre poly and its derivative
            ratio = (Rq_mag(i)/Re_mag).^nm1;
            z = Ren- cosgamma*rqn;
            w = wtemp.*ratio;

            Gterm1 = Pl'*(w.*n);
            Gterm2 = dP'*w;
            G(:,3*i-2:3*i) = Gterm1*rqn + [z(:,1).*Gterm2,z(:,2).*Gterm2,z(:,3).*Gterm2];
        end
        %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%----------------------------------------------------------------------------------%%
        %%% This part computes the potential for a dipole contained within a multi-layer
        %%% isontropic sphere using Berg Parameter Approximation (Zhang Eqs 1i',5i" and 6i)
        %%% The EEG single-shell solution uses the vector form (which avoids computationally
        %%% expensive intrinsic functions) described by Mosher et al ("EEG and MEG: Forward
        %%% Solutions for inverse problems" IEEE BME Trans March 1999)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
    elseif (NL>1)&(method==2)
        %
        if (~exist('mu_berg_in')|~exist('lam_berg_in'))
            error('Berg Parameters have not been specified!!!')
        elseif size(mu_berg_in)~=size(lam_berg_in)
            error('Berg Parameters are of Unequal Lengths!!!')
        else
            J = length(mu_berg_in);
            mu_berg = mu_berg_in;
            lam_berg = lam_berg_in;
        end
        %
        Re_mag = repmat(R(NL),P,M);           %(PxM)
        Re_mag_sq = repmat(R(NL)*R(NL),P,M);  %(PxM)
        Rq_mag = rownorm(Rq);                 %(Px1)
        Rq_mag_sq = Rq_mag.*Rq_mag;           %(Px1)
        Re_dot_Rq = Rq*Re';                   %(PxM)
        %
        for k=1:J
            %
            Rq1 = mu_berg(k)*Rq;                            %(Px3)
            Rq1_mag = mu_berg(k)*Rq_mag;                    %(Px1)
            Rq1_mag_sq = (mu_berg(k)*mu_berg(k))*Rq_mag_sq; %(Px1)
            Re_dot_Rq1 = mu_berg(k)*Re_dot_Rq;              %(PxM)
            %
            const = 4.0*pi*sigma(NL);
            const1 = const/lam_berg(k);
            term = 1./(const1*Rq1_mag_sq);                  %(PxM)
            %
            %%% This part checks for the presence of Berg dipoles which are external to
            %%% the sphere. For those dipoles external to the sphere, the dipole parameters
            %%% are replaced with the electrical image (internal to sphere) of the dipole
            %
            nx = find(Rq1_mag > R(NL));
            %
            if nx>0
                Rq1_temp = Rq1(nx,:);
                Rq1(nx,:) = R(NL)*R(NL)*Rq1_temp./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,3);
                Rq1_mag(nx,1) = rownorm(Rq1(nx,:));
                Rq1_mag_sq(nx,1) = Rq1_mag(nx,1).*Rq1_mag(nx,1);
                Re_dot_Rq1(nx,:) = R(NL)*R(NL)*Re_dot_Rq1(nx,:)./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,M);
                term(nx,:) = (R(NL)/rownorm(Rq1_temp))'.*term(nx,:);
            end
            %
            %%% Calculation of Forward Gain Matrix Contribution due to K-th Berg Dipole
            %
            Rq1_mag = repmat(Rq1_mag,1,M);                %(PxM)
            Rq1_mag_sq = repmat(Rq1_mag_sq,1,M);          %(PxM)
            term = repmat(term,1,M);                      %(PxM)
            %
            d_mag = reshape( rownorm(reshape(repmat(Re,1,P)',3,P*M)' ...
                -repmat(Rq1,M,1)) ,P,M);        %(PxM)
            d_mag_cub = d_mag.*d_mag.*d_mag;                        %(PxM)
            %
            F_scalar = d_mag.*(Re_mag.*d_mag+Re_mag_sq-Re_dot_Rq1); %(PxM)
            %
            c1 = term.*(2*( (Re_dot_Rq1-Rq1_mag_sq)./d_mag_cub) ...
                + 1./d_mag - 1./Re_mag);                 %(PxM)
            c2 = term.*((2./d_mag_cub) + (d_mag+Re_mag)./(Re_mag.*F_scalar));
            %
            G = G +    reshape(repmat((c1 - c2.*Re_dot_Rq1)',3,1),M,3*P) ...
                .*repmat(reshape(Rq1',1,3*P),M,1) ...
                +  reshape(repmat((c2.*Rq1_mag_sq)',3,1),M,3*P) ...
                .*repmat(Re,1,P);
        end
        %
        %---------------------------------------------------------------------------------------
        %
    end    % End of Check for Forward Model Calculation Method
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON