Code covered by the BSD License

# inuctancies for wire loops

### Maxim Vedenyov (view profile)

calculates self inductance and mutual inductance

L=self_inductance(x,y,r)
```function L=self_inductance(x,y,r)

hsf=0.01; % high scale factor, dge*hsf is step in rough integration and in current discretization

ff=0.3; % fine tune factor, ff*r ~ step in finetune (along contour and into deep)
df=5*2; % fintune depth factor to deep, df*r=maximal step in finetune to deep (to inside direction)

sif=ff*r; % ~ step in finetune

ned=length(x); % number of edges

% limits:
x1=min(x);
x2=max(x);
x12=x2-x1;
y1=min(y);
y2=max(y);
y12=y2-y1;

dge=max([x12 y12]); % size estimation, biggest size of box-container

xe=[x x(1)];
ye=[y y(1)];

dx=diff(xe);
dy=diff(ye);
dL=sqrt(dx.^2+dy.^2);
dxn=dx./dL;
dyn=dy./dL;
cL=cumsum(dL);
cLend=cL(end); % total length
cLee=[0 cL];

% find vectors bisctis:

dxn0=[dxn(end) dxn(1:end-1)];
dyn0=[dyn(end) dyn(1:end-1)];

% mean vectors:
xm=dxn+dxn0;
ym=dyn+dyn0;
Lm=sqrt(xm.^2+ym.^2);
xmn=xm./Lm;
ymn=ym./Lm;

% bissectris vectors is perpendicular to mean vectors:
xbn0=-ymn;
ybn0=xmn;

% orient bistectiss vectors to inised:
si=dge*1e-3; % step to inside
xt=x+si*xbn0;
yt=y+si*ybn0; % test points
in = inpolygon(xt,yt,x,y);
sg=2*in-1; % bool to sign
xbn=sg.*xbn0;
ybn=sg.*ybn0;

% means of edges:
xme=(xe(1:end-1)+xe(2:end))/2;
yme=(ye(1:end-1)+ye(2:end))/2;

% vectors perpendicular to edges:
xpn0=-dyn;
ypn0=dxn;
xt=xme+si*xpn0;
yt=yme+si*ypn0; % test points
in = inpolygon(xt,yt,x,y);
sg=2*in-1; % bool to sign
xpn=sg.*xpn0;
ypn=sg.*ypn0;

% si=dge*0.1;
% plot(xe,ye,'r-');
% hold on;
% quiver(xme,yme,si*xpn,si*ypn);
% axis equal;

% dcf distance conversion coefficients
ca2=xbn.*xpn+ybn.*ypn; % cos(alpha/2), where apha angle between to edges
ica2=1./ca2; % 1/cos(alpha/2);

%si=dge*0.02;
%plot(xe,ye,'r.-');
%hold on;
%xi=x+si*ica2.*xbn;
%yi=y+si*ica2.*ybn;
%plot(xi,yi,'b.-');
%quiver(xme,yme,si*xpn,si*ypn);
%axis equal;

% finetune region points and areas:
%xt=interp1(cLee,xe,t);
%yt=interp1(cLee,ye,t);
%zt=interp1(cLee,ze,t);
dca=r:sif:r+df*r;
ddca=dca(2)-dca(1); % height of trapecia is constant
dca1=dca(1:end-1);
dca2=dca(2:end);
dcam=(dca1+dca2)/2;
Ldca=length(dca);
tnpe=(cLend/sif)*Ldca; % total number of finetune points estimation
tnpe=round(tnpe*2); % with margin
xf=zeros(1,tnpe);
yf=zeros(1,tnpe);
Af=zeros(1,tnpe);
fcc=0; % points count
% points
for fc=1:ned % for each edge

% bisctise vectors on end points:
xbn1=xbn(fc);
ybn1=ybn(fc);
ica21=ica2(fc);
xx1=x(fc);
yy1=y(fc);
if fc==ned
xbn2=xbn(1);
ybn2=ybn(1);
ica22=ica2(1);
xx2=x(1);
yy2=y(1);
else
xbn2=xbn(fc+1);
ybn2=ybn(fc+1);
ica22=ica2(fc+1);
xx2=x(fc+1);
yy2=y(fc+1);
end

dLt=dL(fc); % edge length
np=round(dLt/sif); % number of points, np+1 in total, np-1 without limits points
if np==0
np=1;
end
np1=np+1; % total number of points
dvf=linspace(0,1,np1);% dividing factors
dvfc=(dvf(1:end-1)+dvf(2:end))/2; % deviding factors, centers
for dcc=1:Ldca-1
dcamt=dcam(dcc); % distance of paralle shift of edge

xt1=xx1+ica21*dcamt*xbn1;
yt1=yy1+ica21*dcamt*ybn1;

xt2=xx2+ica22*dcamt*xbn2;
yt2=yy2+ica22*dcamt*ybn2;

dxt=xt2-xt1;
dyt=yt2-yt1;

xtt=xt1+dvfc*dxt;
ytt=yt1+dvfc*dyt; % points -centers of areas

xf(fcc+1:fcc+np)=xtt;
yf(fcc+1:fcc+np)=ytt;

dca1t=dca1(dcc); % distance of paralle shift of edge
dca2t=dca2(dcc); % distance of paralle shift of edge
xtt=xt1+dvf*dxt;
ytt=yt1+dvf*dyt; % edges of areas

% from first edge:
xt1=xx1+ica21*dca1t*xbn1;
yt1=yy1+ica21*dca1t*ybn1;

xt2=xx2+ica22*dca1t*xbn2;
yt2=yy2+ica22*dca1t*ybn2;

dxt=xt2-xt1;
dyt=yt2-yt1;
e1L=sqrt(dxt^2+dyt^2);

% from second edge:
xt1=xx1+ica21*dca2t*xbn1;
yt1=yy1+ica21*dca2t*ybn1;

xt2=xx2+ica22*dca2t*xbn2;
yt2=yy2+ica22*dca2t*ybn2;

dxt=xt2-xt1;
dyt=yt2-yt1;
e2L=sqrt(dxt^2+dyt^2);

At=((e1L/np+e2L/np)/2)*ddca;
Af(fcc+1:fcc+np)=At;

fcc=fcc+np;

end

end

xf=xf(1:fcc);
yf=yf(1:fcc);
Af=Af(1:fcc);

% plot(xe,ye,'r-');
% hold on;
% plot(xf,yf,'k.');

% rough points:
dcae=dca(end); % need final poligon
fpx=x+dcae*ica2.*xbn;
fpy=y+dcae*ica2.*ybn;
% plot(fpx,fpy,'b-');

% plot(xr,yr,'g.');

% integration:
% new parametrization:
%dt=hsf*cLend;
dt=hsf*dge;
nt=round(cLend/dt);
%t=0:dt:cLend;
t=linspace(0,cLend,nt); % closed
xi=interp1(cLee,xe,t);
yi=interp1(cLee,ye,t); % current discretization points

dxi=diff(xi);
dyi=diff(yi);
dLi=sqrt(dxi.^2+dyi.^2);
dxin=dxi./dLi;
dyin=dyi./dLi;
% mean points:
xim=(xi(1:end-1)+xi(2:end))/2;
yim=(yi(1:end-1)+yi(2:end))/2;

ni=length(dxin);
Hf=zeros(size(xf)); %field from finetune
Hr=zeros(size(xr)); %field from rough
for cc=1:ni % current pieces counting

ximt=xim(cc);
yimt=yim(cc);

dxit=dxi(cc);
dyit=dyi(cc);

%finetune:

% r0-r:
drx=xf-ximt;
dry=yf-yimt;
drL=sqrt(drx.^2+dry.^2);

% cross product:
cp=dxit*dry-dyit*drx;

% Biot?avart law
Hf=Hf+cp./(drL.^3); % without 1/(4*pi), without I

%rough:

% r0-r:
drx2=xr-ximt;
dry2=yr-yimt;
drL2=sqrt(drx2.^2+dry2.^2);

% cross product:
cp2=dxit*dry2-dyit*drx2;

% Biot?avart law
Hr=Hr+cp2./(drL2.^3); % without 1/(4*pi), without I

end

% flux is sum with areas
Ff=sum(Hf.*Af); % without 1/(4*pi), without I, without mu

mu0=4*pi*1e-7;

Lf=mu0*Ff/(4*pi);