function G = gainmatrix_TPSP(r,rho,phi)
%function G = gainmatrix_TPSP(r,rho,phi)
%
%--------------------------------------------------------------------------%
%
% Fonction calculant le vecteur de transfert de l'activité du dipôle électrique 
% à la position rho avec comme orientation le vecteur phi, en utilisant un modèle
% de tête à trois sphères concentriques homogènes
%
%--------------------------------------------------------------------------%
%
% ENTREES :
% r     <-->  matrice (3xN+1) des positions des N+1 capteurs en coordonnées
% cartésiennes
% rho   <-->  vecteur colonne position du dipôle considéré
% phi    <-->  vecteur colonne d'orientation du dipôle considéré
%
% SORTIE :
% G   <-->  vecteur de transfert 
%
sigma     = [0.0033]; %(en siemens/cm) [Cosandier-Rimélé, 2006]
mu        = [0.6342,0.9364,1.0362];  % [Cosandier-Rimélé, 2006]
ll        = [0.5979,0.2037,0.0237];  % [Cosandier-Rimélé, 2006]
%--------------------------------------------------------------------------%


%déclaration de variables
N = size(r,2) - 1;

%switching matrix : permet de calculer les différences de potentiels
%V = [diag(ones(1,N)),-ones(N,1)]; 
%On ne calcule pas de différences de potentiels ici ou autrement dit
%On calcule une différence de potentiels avec une référence constante
%arbitraire identique pour toutes les sources

%calcul de H
rhomu = rho * mu;
H = [];
for m = 1 : length(mu)
    RHOm = rhomu(:,m) * ones(1,N+1);
    rRHOm = r - RHOm;
    rRHOm2 = rRHOm'*rRHOm; 
    c1 = 2 * rRHOm' * rhomu(:,m) .* diag(rRHOm2).^(-3/2) + diag(rRHOm2).^(-1/2) - diag(r'*r).^(-1/2); 
    c2 = 2 * diag(rRHOm2).^(-3/2) + (diag(rRHOm2).^(-1/2) + diag(r'*r).^(-1/2)) .* (diag(r'*r).^(1/2) .* ...
        diag(rRHOm2).^(1/2) .* ( diag(r'*r).^(1/2) .* diag(rRHOm2).^(1/2) + diag(r'*r) - r'*rhomu(:,m))).^(-1); 
    h = ((c1 - c2 .* (r'*rhomu(:,m))) * rhomu(:,m).' + rhomu(:,m)'*rhomu(:,m) * c2*ones(1,3) .* r.')...
        /(4*pi*sigma*rhomu(:,m)'*rhomu(:,m));
    H = [H;h];    
end

%calcul de G
%G = V * (kron(ll,eye(N+1))) * H;
G = kron(ll,eye(N+1)) * H;

G = G*phi;