function [Te_ISA,Te,Tm_ISA,Tm] = trans_matrix(mode,test_opt,ISA,DIM,s,I,H,P,A,no_megsensor)
%TRANS_MATRIX - calculate BEM transfer matrices for EEG and MEG
% function [Te_ISA,Te,Tm_ISA,Tm] = trans_matrix(mode,test_opt,ISA,DIM,s,I,H,P,A,no_megsensor)
%****************************************************************************
%
% function trans_matrix(mode,N,s,no_megsensor)
% this function calculates transfer matrices for EEG and MEG
% mode        =1 , compute Te
%              2 , compute Tm
%              3 , compute BOTH
% ISA         see bem.m
% DIM         : DIM(i) is the number of triangles(or nodes) on surface i
%             : 3 x 1
% s           : [conductivity of brain, conductivity of skull]
%             : 2 x 1
% P           : EEG scan matrix       no_eegsensor x sum(DIM)
% H,A         :
% no_megsensor: number of MEG sensors
%
% The isolated skull approach is incorporated for THREE shells model.
%*****************************************************************************

%<autobegin> ---------------------- 26-May-2004 11:34:40 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% At Check-in: $Author: Mosher $  $Revision: 14 $  $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:40 -----------------------


s3=s(1);
s2=s(2);
no_elt=sum(DIM); 

%disp(' LU decomposition .......')
% Incoporate the isolated skull approach 
% C3: solid angle matrix of brain layer

if ISA, 
   M=DIM(1)+DIM(2);
   b =  M+1:no_elt;  
   I3 = I(b,b);
   C3 =I3 - H(b,b)*(s3+s2)/(s3-s2);

   for i=1:size(C3,1)
     C3(i,i)=0;
     C3(i,i)= -sum(C3(i,:));
   end
   C3 = C3+1/DIM(3);   % deflation 

   [L3,U3]=lu(C3);
   L3=sparse(L3);  U3=sparse(U3);     
   clear C3 b
end   

H = I- H + 1/no_elt;  
clear I
[L,U] = lu(H);
%
L=sparse(L);    U=sparse(U);
clear H 

Te_ISA = [];   %  transfer matrix for EEG, with isolated skull approach(ISA)
Tm_ISA = [];   %  transfer matrix for MEG, with ISA 
Te = [];   % without isa
Tm = [];

if mode==1 | mode==3,
   %disp(' ')
   %disp('>>>> Compute Te ...')
%   t=clock;
   Te = P/U/L; 
   if ISA   % isolated skull approach
      P3 = P(:,M+1:no_elt);    
      Te_ISA = s2/s3*[Te(:,1:M),Te(:,M+1:no_elt)-2*(Te(:,M+1:no_elt)*I3)/U3/L3+(1+s3/s2)*P3/U3/L3];
    end
%   run_Te_in_min=etime(clock,t)/60
%keyboard
   clear P P3
end; 

if mode==2 | mode==3,
   %disp(' ')
   %disp('>>>> Compute Tm .....')
%   t=clock;
   Tm = A/U/L;
   if ISA
     A3 = A(:,M+1:no_elt);
     clear U L A
     Tm_ISA = s2/s3*[Tm(:,1:M),Tm(:,M+1:no_elt)-2*(Tm(:,M+1:no_elt)*I3)/U3/L3+(1+s3/s2)*A3/U3/L3];
   end
%   run_Tm_in_min=etime(clock,t)/60
end;
%keyboard








