funcprot(0);

mu0 = 1.2566370614e-6;

Af = 1e-4;	// ferrite

lf = 15e-2;

Nt = 300;

f = 50;
w = 2*%pi*f;
C = 35e-6;

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

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.001;
Imax=0.05;

N = 1000;

idata = linspace ( Imin, Imax, N );
lu_data = zeros(1,N);
cu_data = zeros(1,N);
u_data = zeros(1,N);

for i=1:N
  lu_data(i) = L( idata(i) ) * w * idata(i);
  cu_data(i) = idata(i) / ( C * w );
  u_data(i) = abs(cu_data(i) - lu_data(i));
end


clf();
subplot(111); a=get("current_axes"); a.grid=[1,1];
plot ( idata , lu_data, "r");
plot ( idata , cu_data, "b");
plot ( idata , u_data, "g");
xtitle("u(i)");

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