function [RoDip,ThetaDip,PhiDip] = Sphere4DipPot(ECart,DipCart,DipIter)
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      

%==================================================================
%
%	Sphere4DipPot.m
%
%	Markus Junghoefer	[1995]
%
%
%	
%	Function definition
%	Eingabe:	ECart 		ist eine mx3 Matrix der Elektroden Positionen
%	Eingabe:	DipCart 	ist eine nx3 Matrix der Dipol Positionen
%	Ausgabe:	RoDip 		ist eine nxm Matrix der radialen Leadfields
%	Ausgabe:	ThetaDip 	ist eine nxm Matrix der polaren Leadfields
%	Ausgabe:	PhiDip 		ist eine nxm Matrix der azimutalen Leadfields
%	
%=================================================================

if nargin<1;		%Test
	DipIter=15;
	n_scalp=20;
	total_scalp_locations=n_scalp.^2;
	[X_scalp, Y_scalp, Z_scalp]=sphere(n_scalp-1);
	[ECart]=change_nn_n3(X_scalp,Y_scalp,Z_scalp);
	DipCart=[0 0 .5];
end

fprintf('\n');
[SizeECart1,SizeECart2]=size(ECart);
if SizeECart1==3 & SizeECart2~=3
	ECart=ECart';
	SizeECart1=SizeECart2;
	SizeECart2=3;
elseif SizeECart1~=3 & SizeECart2~=3
	fprintf('Bad size of ECart in Sphere4DipPot')
	return;
end
[SizeDipCart1,SizeDipCart2]=size(DipCart);
if SizeDipCart1==3 & SizeDipCart2~=3
	DipCart=DipCart';
	SizeDipCart1=SizeDipCart2;
	SizeDipCart2=3;
elseif SizeDipCart1~=3 & SizeDipCart2~=3
	fprintf('Bad size of DipCart in Sphere4DipPot')
	return;
end
%=============================================

%DEFINE RADII OF SHELLS 
cort_radius = 0.0815; %radius of shell 1 (brain or cort)
csf_radius=0.0836;%radius of shell 2 (cerebro-spinal fluid)
skull_radius=0.0878; %radius of shell 3 (skull)
scalp_radius=0.092; %radius of shell 4 (scalp)


% DEFINE CONDUCTIVITY VARIABLES (ohms * meters)
ocort=1/3;  %sigma1 of shell 1 (brain or cort)   = 0.333
ocsf=3 * ocort;    %sigma2 of shell 2 (cerebro-spinal fluid)  =1.0
oskull= ocort / 80;   %sigma3 of shell 3 (skull) = 1/240 = 0.0042
oscalp = ocort; % sigma4 of shell 4 (scalp) = 0.333

% VARIABLES DEFINING RATIO OF CONDUCTIVITY SHELLS (unitless)
b = cort_radius / scalp_radius;   
c = csf_radius / scalp_radius;    
d = skull_radius / scalp_radius;   

% CONDUCTIVITY	RATIOS (unitless)
k1 = ocort / ocsf;
k2 = ocsf / oskull;
k3 = oskull / oscalp;

% DIPOLE LOCATION RATIO (unitless)	
%f = dipole_radius / scalp_radius;   

%=============================================================================

[DipSphere] = change_sphere_cart(DipCart,[],-1);
[EPosSphere] = change_sphere_cart(ECart,[],-1);
StdEPosRadius=std(EPosSphere(:,3));
fprintf('Std of radius of electrode positions: %g\n\n',StdEPosRadius);
MeanRadius=mean(EPosSphere(:,3));
RadiusFactor=scalp_radius./MeanRadius;
DipSphere(:,3)=DipSphere(:,3).*RadiusFactor;
EPosSphere(:,3)=EPosSphere(:,3).*RadiusFactor;
if any(DipSphere(:,3)>=cort_radius)
	fprintf('Error: Dipole radius has to be smaller than %g per hundred of scalp radius',round(b.*100))
	return;
end

%=============================================

RoDip=zeros(SizeDipCart1,SizeECart1);
ThetaDip=zeros(size(RoDip));
PhiDip=zeros(size(RoDip));
LegRad=zeros(DipIter,1);
LegTang=zeros(DipIter,1);

%Transformiere kartesische Dipfield Punkte, damit die Elektrodenpositionen auf Nordpol
%Drehung von phi um z-Achse mit der Matrix
%TransMat1=[cos(phi) sin(phi) 0 ; -sin(phi) cos(phi) 0 ; 0 0 1];
%Drehung von theta um y-Achse mit der Matrix 
%TransMat2=[cos(theta) 0 -sin(theta) ; 0 1 0 ; sin(theta) 0 cos(theta)]
	
for i=1:SizeDipCart1
	
	Theta=DipSphere(i,1);
	CosTheta=cos(Theta);
	SinTheta=sin(Theta);
	Phi=DipSphere(i,2);
	CosPhi=cos(Phi);
	SinPhi=sin(Phi);	
	
	TransMat1=[CosPhi SinPhi 0 ;-SinPhi CosPhi 0 ;0 0 1];
	TransMat2=[CosTheta 0 -SinTheta;0 1 0 ;SinTheta 0 CosTheta];
		
	TransMat=TransMat2*TransMat1;
	DipCartTmp=(TransMat*DipCart')';
	[DipSphereTmp] = change_sphere_cart(DipCartTmp,[],-1);
	ECartTmp=(TransMat*ECart')';
	[EPosSphereTmp] = change_sphere_cart(ECartTmp,[],-1);
	
	f=DipSphereTmp(i,3).*RadiusFactor./scalp_radius;  

	for j=1:SizeECart1
		CosEPosTheta=cos(EPosSphereTmp(j,1));
		for k=1:DipIter	
			LegRad(k)= calc_legendre(k,0,CosEPosTheta);
			LegTang(k)=calc_legendre(k,1,CosEPosTheta);

			faktor = (d ^ (2 * k + 1)) * ...                                                 
			((b ^ (2 * k + 1)) * k * (k1 - 1) * (k2 - 1) * (k + 1) +...         
  			(c ^ (2 * k + 1)) * (k1 * k + k + 1) * (k2 * k + k + 1)) *...     
			((k3 * k + k + 1) + (k + 1) * (k3 - 1) * (d ^ (2 * k + 1))) + ... 
			 (k + 1) * (c ^ (2 * k + 1)) *...                                              
			((b ^ (2 * k + 1)) * (k1 - 1) * (k2 * k + k2 + k) +...                
			(c ^ (2 * k + 1)) * (k1 * k + k + 1) * (k2 - 1)) *...                  
			 (k * (k3 - 1) + (k3 * k + k3 + k) * (d ^ (2 * k + 1)));
			FaktorRad(k)  =  (1/(4 * pi * oscalp * scalp_radius^2))*((2 * k) + 1)^4 * f^(k - 1) * (c * d)^(2 * k + 1)  / faktor;
			FaktorTang(k) =  (1/(4 * pi * oscalp * scalp_radius^2))*((2 * k) + 1)^4 * f^(k - 1) * (c * d)^(2 * k + 1)  / (k*faktor);
		end
		RoDip(i,j)=FaktorRad*LegRad;
		TangTmp=FaktorTang*LegTang;
		ThetaDip(i,j)=cos(EPosSphereTmp(j,2)).*TangTmp;
		PhiDip(i,j)=sin(EPosSphereTmp(j,2)).*TangTmp;
	end		
end
if nargin<1; 
	colormap(jet(128));
	[RoDip]=change_n3_nn(RoDip);
	[ThetaDip]=change_n3_nn(ThetaDip);
	[PhiDip]=change_n3_nn(PhiDip);
	subplot(2,2,1)
	surf(X_scalp, Y_scalp, Z_scalp,RoDip);
	title('Ro Leadfield')
	subplot(2,2,2)
	surf(X_scalp, Y_scalp, Z_scalp,ThetaDip);
	title('Theta Leadfield')
	subplot(2,2,3)
	surf(X_scalp, Y_scalp, Z_scalp,PhiDip);
	title('Phi Leadfield')
end

return;
