xlm.m

function XM=xlm(Rp,Zp,Rc,Zc,DRc,DZc,iself);
%
% Compute self and mutual inductance
%
% function XM=xlm(Rp,Zp,Rc,Zc,DRc,DZc,iself);
% compute flux (Wb) at the points (Rp,Zp) due to unit currents (A) in (Rc,Zc)
% XM(length(Rp),length(Rc)): each row refers to a point, each column to a coil
%       
if iself~=1
%
% compute mutual inductance
%
DZ = Zp*ones(size(Zc'))-ones(size(Zp))*Zc';
RA = Rp*ones(size(Rc'))+ones(size(Rp))*Rc';
K2 = 4*(Rp*Rc')./(RA.*RA+DZ.*DZ);
[J1,J2]=ellipke(K2);
JW = (1-0.5*K2).*J1-J2;
XM = 4*pi*1.e-7*JW.*sqrt(RA.*RA+DZ.*DZ);
return
%
else
%
% compute self inductance
%
XS    = zeros(length(Rc),1);
RAT   = DZc./DRc;
ALF   = DZc./Rc;  	

icirc = find(RAT < 0.0);
idisk = find(RAT < 1.0 & RAT >= 0.0);
ilong = find(ALF > 1.5);
ishor = find(ALF <=1.5 & RAT >= 1.0);

if isempty(icirc)==0; clear ALFA ALFA2 BETA BETA2 GAMMA GAMMMA2 ALG8 ASP8 AASP8 RATIO;
RMIN  = sqrt(DRc(icirc).*DZc(icirc)/pi);
ASP8  = 8*Rc(icirc)./RMIN; 		AASP8 = log(ASP8);
XS(icirc) = AASP8-1.75+8*(AASP8+1/3)./(ASP8.*ASP8);
end

if isempty(idisk)==0; clear ALFA ALFA2 BETA BETA2 GAMMA GAMMMA2 ALG8 ASP8 AASP8 RATIO;
ALFA  = ALF(idisk);  			 ALFA2  = ALFA.*ALFA;  	
BETA  = DRc(idisk)./ Rc(idisk);  BETA2  = BETA.*BETA;  		ALG8  = log(8./BETA);
GAMMA = DRc(idisk)./DZc(idisk);  GAMMA2 = GAMMA.*GAMMA;  	
XS(idisk) = ALG8.*(1+(1/96)*BETA2)-0.5+(43/576)*BETA2+(1/32)*ALFA2.*(1+ALG8);
XS(idisk) = XS(idisk) -(pi/3+pi/60*ALFA2)./GAMMA+(25/72+(23/1152)*ALFA2+(1/3)*log(GAMMA))./GAMMA2;
end

if isempty(ilong)==0; clear ALFA ALFA2 BETA BETA2 GAMMA GAMMMA2 ALG8 ASP8 AASP8 RATIO;
ALFA  = ALF(ilong); 			ALFA2  = ALFA.*ALFA; 
AM    = 1./(sqrt(1+ALFA2));		AM3    = AM.*AM.*AM; AM5  = AM.*AM.*AM3; AM7 = AM.*AM.*AM5;
BETA  = DRc(ilong)./Rc(ilong);  BETA2  = BETA.*BETA;
GAMMA = DRc(ilong)./DZc(ilong);  	
XS(ilong) = (1./AM-8/3/pi-AM3/8+AM5/16)./ALFA;
XS(ilong) = XS(ilong) - BETA/3+BETA2/12+((log(8./BETA)-23/12)/6/pi+(AM/12-AM3*5/48+AM5/6-AM7*95/256)).*BETA.*GAMMA;
XS(ilong) = pi*XS(ilong)./ALFA;
end

if isempty(ishor)==0; clear ALFA ALFA2 BETA BETA2 GAMMA GAMMMA2 ALG8 ASP8 AASP8 RATIO;
ALFA  = ALF(ishor); 			ALFA2  = ALFA.*ALFA; 		ALFA4  = ALFA2.*ALFA2; 
BETA  = DRc(ishor)./Rc(ishor);  BETA2  = BETA.*BETA;
GAMMA = DRc(ishor)./DZc(ishor); GAMMA2 = GAMMA.*GAMMA;		GAMMA4 = GAMMA2.*GAMMA2;  	
RATIO = RAT(ishor);
XS(ishor) = log(8./ALFA).*(1+ALFA2/32-ALFA4/1024)-0.5+ALFA2/128+ALFA4/1536;
XS(ishor) = XS(ishor)-GAMMA*pi/3+GAMMA2*25/72+BETA2/32+GAMMA4/180+BETA2/96.*log(8./ALFA)+GAMMA2.*log(RATIO)/6;
end

XS = 4e-7 * pi * Rc.*XS;

DZ = Zp*ones(size(Zc'))-ones(size(Zp))*Zc';
RA = Rp*ones(size(Rc'))+ones(size(Rp))*Rc';
K2 = 4*(Rp*Rc')./(RA.*RA+DZ.*DZ);
K2 = triu(K2,1);
[J1,J2]=ellipke(K2);
JW = (1-0.5*K2).*J1-J2;
XM = 4*pi*1.e-7*JW.*sqrt(RA.*RA+DZ.*DZ);
XM = triu(XM,1);
XM = XM+XM'+diag(XS);
return;
end