function [csd_gamma_cort] = CalcCsdCortGamma(i,para_vec)

%	CalcCsdCortGamma.m

%   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.    


%DEFINE RADII OF SHELLS
cort_radius= para_vec(1);
csf_radius= para_vec(3);
skull_radius= para_vec(5);
scalp_radius= para_vec(7);
dipole_radius=para_vec(9);

% DEFINE CONDUCTIVITY VARIABLES (ohms * meters)
o1 = para_vec(2); %ocort;      
o2 = para_vec(4); %ocsf;    
o3 = para_vec(6); %oskull; 
o4 = para_vec(8); %oscalp;       

% 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 = o1 / o2;
k2 = o2 / o3;
k3 = o3 / o4;

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

		
% CALCULATE gamma

	gamma_tmp = (d ^ (2 * i + 1)) * ...                                                 
				((b ^ (2 * i + 1)) * i * (k1 - 1) * (k2 - 1) * (i + 1) +...         
   				(c ^ (2 * i + 1)) * (k1 * i + i + 1) * (k2 * i + i + 1)) *...     
				((k3 * i + i + 1) + (i + 1) * (k3 - 1) * (d ^ (2 * i + 1))) + ... 
				 (i + 1) * (c ^ (2 * i + 1)) *...                                              
				((b ^ (2 * i + 1)) * (k1 - 1) * (k2 * i + k2 + i) +...                
				(c ^ (2 * i + 1)) * (k1 * i + i + 1) * (k2 - 1)) *...                  
				 (i * (k3 - 1) + (k3 * i + k3 + i) * (d ^ (2 * i + 1)));

	% CALCULATE OMEGA

		omega=(d^(2*i+1)*i*(k3*i+i+1) + d^(2*(2*i+1))*i*(i+1)*(k3-1))*...
			  (c^(2*i+1)*(k1-1)*(k2*i+i+1) + b^(2*i+1)*(k1*i+k1+i)*(k2-1)) +...
			  (c^(2*(2*i+1))*i*(i+1)*(k2-1)*(k1-1) + b^(2*i+1)*c^(2*i+1)*(k1*i+k1+i)*(k2*i+k2+i)) *...
			  (c^(2*i+1)*i*(k3-1)+c^(2*i+1)*d^(2*i+1)*(k3*i+k3+i));

		% CALCULATE gamma_CORT
		%cause ÆPn=n*(n+1)*Pn
		csd_gamma_cort=  i*(i+1) *...	
					 (1/(4 * pi * o1 * (scalp_radius^2))) * ...
        			 (1 / (b * f)) * ((f / b) ^ i) * ...
       			     ((omega * (i + 1) / gamma_tmp) + i);
				 

return;	
	

	
