function bem_gain(bem_grid_mfname,bem_xfer_mfname,ISA,bem_gaingrid_mfname, Verbose)
%BEM_GAIN - Computes the EEG/MEG forward/gain matrix associated with a set of grid points 
% function bem_gain(bem_grid_mfname,bem_xfer_mfname,ISA,bem_gaingrid_mfname, Verbose)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% EEG/MEG BEM 3-D  GRID FORWARD GAIN CALCULATION UTILITY 
%
% This function computes the EEG/MEG forward matrix (kernel) associated with a set of 
% grid points computed using the utility "gridmaker.m".    
%
% INPUTS (Required):
%
%     bem_grid_mfname: can be either 
%                      - a "*.mat" file name containing pre-computed 3-D grid information   
%                      - a structure of specific source locations (3xN  array .Loc)and orientations (3xN  array .Orient)
%                          
%     bem_xfer_mfname: "*.mat" file containing pre-computed BEM transfer matrix   (char string) 
%                 ISA: Isolated Skull Approach (0=disable;1=enable)               (1 x 1)
% bem_gaingrid_mfname: ".mat" output filename used to store final gain matrix and 
%                      associated output parameters                               (char_string)
%             Verbose: Toggles Verbose mode on (1) / off (0)
% INPUTS (Optional):
%                -None-
%
% OUTPUTS: Complete set of output parameters are stored in user specified 
%          "*.mat" file defined by "bem_gaingrid_mfname.mat" 
%
%   Other Notes: 3x3 MEG kernel is dotted with sensor orientation to generate more compact
%                1x3 kernel 
%
% External Functions and Files
%     USC/LANL/CNRS BrainStorm Toolbox
%    
%     - John J. Ermer 3/7/00 
%     - John J. Ermer 4/28/00 (Added new kernel function to interpolate when using linear basis)
%     - John J. Ermer 5/05/00 (Added new logic to generate grids for MEG)
%     - Sylvain Baillet 25/03/02 (Added Verbose mode)  
%     - SB 15-May-2002 Updated Verbose mode with new bst_message_window overwrite features
%     - SB 07-Mar-2003 Renamed to bem_gain.m 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%<autobegin> ---------------------- 14-Jun-2004 17:09:45 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bem_kernel.m
%   toolbox\bst_message_window.m
%
% At Check-in: $Author: Mosher $  $Revision: 18 $  $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:09:45 -----------------------


%
%%%% THIS PART LOADS AND OPENS INPUT DATA FILES %%%%
%
if ischar(bem_grid_mfname) % Locations of grid points are specified in a file
    if ~exist(bem_grid_mfname,'file') 
        errordlg(sprintf('BEM Grid File %s does not exist on this platform - please re-compute it',bem_grid_mfname))
        return
    end
    load(bem_grid_mfname);  % Load Pre-computed3-D Grid Parameters
    
elseif isstruct(bem_grid_mfname) % Point locations are specified in a structure 
  
    Rq_bemgrid = bem_grid_mfname.Loc';
    
else
    errordlg('Location of grid points are passed with wrong argument format',sprintf('Error when calling %s',mfilename))
    return
end

if ~exist(bem_xfer_mfname)
    errordlg(sprintf('BEM Transfer Matrix File %s does not exist on this platform - please re-compute it',bem_xfer_mfname));
    return
end

Xfer = load(bem_xfer_mfname);  % Load Pre-computed Transfer Matrices 

%
%%%% THIS PART CHECKS IF INPUT PARAMETERS ARE VALID %%%%
%
if (Xfer.modality<1)|(Xfer.modality>2)
    error(sprintf('Invalid Modality of %1.0f has been specified!!!',Xfer.modality))
end
if (Xfer.modality==1)& ~isfield(Xfer,'Te')
    error(sprintf('EEG Transfer Matrices are not contained in the file: %s',bem_xfer_mfname))
elseif (Xfer.modality==2)& ~isfield(Xfer,'Tm')
    error(sprintf('MEG Transfer Matrices are not contained in the file: %s',bem_xfer_mfname))
end
%
%%%% THIS PART LOADS THE APPROPRIATE TRANSFER MATRIX: STANDARD OR ISA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ISA==0
    if Xfer.modality==1
        Te_save = Xfer.Te;
        Tm_save = [];  % Dummy value for kernel routine
    else
        Tm_save = Xfer.Tm;
        Te_save = [];  % Dummy value for kernel routine
    end
elseif ISA==1
    if Xfer.modality==1
        Te_save = Xfer.Te_ISA;
        Tm_save = [];
    else
        Tm_save = Xfer.Tm_ISA;
        Te_save = [];
    end
else
    if Verbose
        bst_message_window(sprintf('Invalid ISA option chosen!!!'))
    end
    
end
%
clear Xfer.Te Xfer.Tm Xfer.Te_ISA Xfer.Tm_ISA
%
%%%% THIS PART DETERMINES DIMENSION PARAMETERS AND DIVIDES THE INPUT DIPOLE SET INTO A SET OF BLOCKS %%%% 
%
if Xfer.modality==1  %EEG
    M = size(Xfer.R_eeg,1);      % Number of Sensors
    S = [];
    R_meg = [];
    O_meg = [];
elseif Xfer.modality==2  %MEG
    M = size(Xfer.R_meg,1);
    Xfer.R_eeg = [];
    Xfer.P_wts = [];
end
%
P = size(Rq_bemgrid,1);   % Number of Dipole Locations
%
BLK_SIZE = 50;  % Maximum Number of dipoles to process at a single time
% Set to prevent overstepping memory bounds in Kernel Calculation
blk_tot = floor(P/BLK_SIZE);  % Total Number of blocks
blk_rem = rem(P,BLK_SIZE);    % Residual Number of Dipoles
%
%%%% THIS PART COMPUTES THE GAIN MATRIX KERNEL (PROCESSING A BLOCK AT A TIME) %%%%
%
t0 = clock;  % Start the clock a tickin...
if Verbose
    hwait = waitbar(0,'Creating the gain matrix for selected sources. . .');
end

%

if Xfer.modality==1    % EEG
    Reeg = Xfer.eegsensor(:,2:4);
    Rmeg = [];
    GBEM_grid = zeros(M,3*P);   % Pre-initialize Output
    for nblk = 1:blk_tot
        n1 = (nblk-1)*3*BLK_SIZE + 1;  % equivalent dipole start index
        n2 = nblk*3*BLK_SIZE;          % equivalent dipole end index
        n3 = (nblk-1)*BLK_SIZE + 1;    % dipole start index
        n4 = nblk*BLK_SIZE;            % dipole end index
        [GBEM_grid(:,n1:n2) dummy] = bem_kernel(Rq_bemgrid(n3:n4,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,S);
        telap = etime(clock,t0);
        trem = telap*(P-n4)/n4;
        if Verbose & ~rem(nblk,10)
            %bst_message_window('overwrite', sprintf('Completed %.0f of %.0f grid points in %3.2f Sec.',n4,P,telap))
            %Expected time to complete: %.2f Sec....',n4,P,telap,trem))
            waitbar(nblk/blk_tot, hwait);
        end
        
    end
    if blk_rem > 0
        n1 = blk_tot*3*BLK_SIZE + 1;   % equivalent dipole start index (for residual dipoles)
        n2 = 3*P;                      % equivalent dipole end index (for residual dipoles)
        n3 = blk_tot*BLK_SIZE + 1;     % dipole start index
        n4 = P;                        % dipole end index
        [GBEM_grid(:,n1:n2) dummy] = bem_kernel(Rq_bemgrid(n3:P,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,S);
    end
    
elseif Xfer.modality==2   % MEG
    
    Reeg = [];
    GBEM_grid = zeros(M,3*P);       % Pre-initialize Output
    Gtemp = zeros(M,3*BLK_SIZE);
    for nblk = 1:blk_tot
        n1 = (nblk-1)*3*BLK_SIZE + 1;  % equivalent dipole start index
        n2 = nblk*3*BLK_SIZE;          % equivalent dipole end index
        n3 = (nblk-1)*BLK_SIZE + 1;    % dipole start index
        n4 = nblk*BLK_SIZE;            % dipole end index
        Ptemp = BLK_SIZE;              % number dipoles being processed at current time
        [dummy Gtemp] = bem_kernel(Rq_bemgrid(n3:n4,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,Xfer.R_meg);
        %
        for m=1:M
            GBEM_grid(m,n1:n2) = Xfer.O_meg(m,:)*Gtemp(3*m-2:3*m,:);  % Apply Sensor Orient to each lead field
        end
        %                                   
        telap = etime(clock,t0);
        trem = telap*(P-n4)/n4;
        if Verbose & ~rem(nblk,100)
            %bst_message_window('overwrite', sprintf('Completed %.0f of %.0f grid points in %3.1f Sec.',n4,P,telap));
            %Expected time to complete: %.2f Sec....',n4,P,telap,trem))
            waitbar(nblk/blk_tot, hwait);
        end
        
    end
    if blk_rem > 0
        n1 = blk_tot*3*BLK_SIZE + 1;   % equivalent dipole start index (for residual dipoles)
        n2 = 3*P;                      % equivalent dipole end index (for residual dipoles)
        n3 = blk_tot*BLK_SIZE + 1;     % dipole start index
        n4 = P;                        % dipole end index
        Ptemp = n4-n3+1;               % number dipoles being processed at current time
        Gtemp = zeros(M,3*Ptemp);
        [dummy, Gtemp] = bem_kernel(Rq_bemgrid(n3:P,:),Xfer.basis,Xfer.test,Xfer.geometry,Xfer.nodes,Xfer.cdv,Te_save,Xfer.P_wts,Tm_save,Xfer.R_meg);
        %
        for m=1:M 
            GBEM_grid(m,n1:n2) = Xfer.O_meg(m,:)*Gtemp(3*m-2:3*m,:);   % Apply Sensor Orient to each lead field
        end
    end
end

clear Gtemp

% % If required, apply source orientation
% if isstruct(bem_grid_mfname) 
%     if ~isempty(bem_grid_mfname.Orient) % Orientations are actually specified
%         if Verbose & ~isempty(bem_grid_mfname.Orient)
%             bst_message_window('Applying source orientation. . .')  
%             isrc = 0;
%             Gtmp = GBEM_grid; 
%             clear GBEM_grid
%             for src = 1:size(Gtmp,2)/3
%                 isrc = isrc + 1;
%                 GBEM_grid(:,src) = Gtmp(:, 3*(isrc-1)+1:3*isrc) * bem_grid_mfname.Orient(:,src);
%             end
%         end
%     end
% end

telap_fwdgaingrid = etime(clock, t0);
if Verbose
    close(hwait)
    bst_message_window(sprintf('Total time : %3.1f sec.',telap_fwdgaingrid))
end


%
%%%% THIS PART STORES THE FINAL RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
% str1 = 'daz az_min az_max del el_min el_max dRc dRf NRc_min Rmin AZ EL Rgrid Rgrid_N';
% str2 = 'Rmax Rbound c0_bemgrid Rq_bemgrid_aer Rq_bemgrid_c0 Rq_bemgrid c0_bemgrid';
% str3 = 'Xfer.R_eeg R_meg O_meg bem_grid_mfname bem_xfer_mfname Xfer.modality GBEM_grid dip_azel_indx';
%
dip_azel_indx = [];
modality = Xfer.modality;
if isfield(Xfer,'R_eeg')
    R_eeg = Xfer.R_eeg ;
else
    R_eeg = [] ;
end

if isfield(Xfer,'R_meg')
    R_meg = Xfer.R_meg ;
    O_meg = Xfer.O_meg;
else
    R_meg = [];
    O_meg = [];
end

str1 = 'R_eeg R_meg O_meg bem_grid_mfname bem_xfer_mfname modality GBEM_grid dip_azel_indx';
%eval(sprintf('save %s %s %s %s',bem_gaingrid_mfname,str1,str2,str3))
eval(sprintf('save %s %s %s %s',bem_gaingrid_mfname,str1))





