funcprot(0);

mu0 = 1.2566370614e-6;

Af = 1e-4;	// ferrite

lf = 15e-2;

Nt = 30;

//=======================================================================

Bs = 0.38;  //N30
Br = 0.175;
Hc = 12;
g = Bs/Br-1;
Nl = Nt;
lh = lf;

function f = Bup ( H )
  f = Bs * (H+Hc) / (abs(H+Hc)+Hc*g) + mu0*H;
endfunction

function f = Bdn ( H )
  f = Bs * (H-Hc) / (abs(H-Hc)+Hc*g) + mu0*H;
endfunction

function f = Bmag ( H )
  f = (Bup(H)+Bdn(H))/2;
endfunction

function f = H ( I )
  I = abs(I);
  if(I < 0.001) then
    I = 0.001;
  end
  f = Nl * I / lh;
endfunction

function f = muf ( I )
  f = Bmag(H(I))/(H(I)*mu0);
endfunction

//=======================================================================

function f = L ( i )
  Rf = lf / (mu0*muf(i)*Af);
  f = Nt*Nt / Rf;
endfunction

//===============================================================================

Imin=0.01;
Imax=1;

N = 1000;

idata = linspace ( Imin, Imax, N );
l_data = zeros(1,N);
muf_data = zeros(1,N);

for i=1:N
  l_data(i) = L( idata(i) );
  muf_data(i) = muf( idata(i) );
end


clf();
subplot(111); a=get("current_axes"); a.grid=[1,1];
plot ( idata , l_data, "r");
xtitle("L(i)");

//=======================================================================