module magnetic_core(t1, t2);
 inout t1,t2;
 magnetic t1,t2;
 parameter real area=1 from (0:inf);
 parameter real length=1 from (0:inf);
 parameter real gapl=0 from [0:inf);
 parameter real hc=1 from (0:inf);
 parameter real bsat=1 from (0:inf);
 parameter real damping=0 from [0:inf);

 real arg, atanh_arg;
 real permeability;

 analog
 begin
   @(initial_step)
     permeability = `M_PI * 4.0e-7;

// linearize atanh() outside the range of -0.99 to 0.99.
   arg = Phi(t1,t2) / (area * bsat);
   if (arg >= 0.99)
      atanh_arg = 50.25 * (arg - 0.99) + 2.647;
   else
     if (arg < -0.99)
       atanh_arg = 50.25 * (arg + 0.99) - 2.647;
     else
       atanh_arg = atanh( arg );

// Contribution due to nonlinear core
   MMF(t1,t2) <+ length * hc * atanh_arg;
// Contribution due to linear gap
   MMF(t1,t2) <+ gapl * Phi(t1,t2) / (area * permeability);
// Contribution due to loss
   MMF(t1,t2) <+ length * ddt( damping*Phi(t1,t2) / area );
 end
endmodule