model M

  import SI = Modelica.SIunits;

  record State

    SI.Temperature T;

    SI.AbsolutePressure p;

    Real omega;

  end State;



replaceable function equilibriumPressureAdsorption

    input State s;

    output SI.Pressure p_eq;

  end equilibriumPressureAdsorption;



  function equilibriumLoadingAdsorption "Uses modelica's root finding algorithm to calculate the equilibrium loading based on the equilibrium pressure function"

    input State s;

    output Real omega_eq;

  algorithm

    omega_eq:=Modelica.Math.Nonlinear.solveOneNonlinearEquation(function equilibriumLoadingAdsorptionDummy(s = s), 0, 0.02, 1e-006);

  end equilibriumLoadingAdsorption;



  function equilibriumLoadingAdsorptionDummy

    extends Modelica.Math.Nonlinear.Interfaces.partialScalarFunction;

    input State s;

  algorithm

    y:=equilibriumPressureAdsorption(State(T = s.T, p = s.p, omega = u)) - s.p;

  end equilibriumLoadingAdsorptionDummy;



  Real r = equilibriumLoadingAdsorption(State(time,2,3));

  annotation(uses(Modelica(version="3.2.1")));

end M;

 function MetalHydrideStorage.M.State "Automatically generated record constructor for MetalHydrideStorage.M.State"

  input Real T(quantity = "ThermodynamicTemperature", unit = "K", min = 0.0, start = 288.15, nominal = 300.0, displayUnit = "degC");

  input Real p(min = 0.0, nominal = 100000.0, quantity = "Pressure", unit = "Pa", displayUnit = "bar");

  input Real omega;

  output State res;

end MetalHydrideStorage.M.State;



function MetalHydrideStorage.M.equilibriumLoadingAdsorption "Uses modelica's root finding algorithm to calculate the equilibrium loading based on the equilibrium pressure function"

  input MetalHydrideStorage.M.State s;

  output Real omega_eq;

algorithm

  omega_eq := Modelica.Math.Nonlinear.solveOneNonlinearEquation(function MetalHydrideStorage.M.equilibriumLoadingAdsorptionDummy(.MetalHydrideStorage.M.State(#(s.T), #(s.p), #(s.omega))), 0.0, 0.02, 1e-006);

end MetalHydrideStorage.M.equilibriumLoadingAdsorption;



function MetalHydrideStorage.M.equilibriumLoadingAdsorptionDummy

  input Real u "Independent variable";

  output Real y "Dependent variable y=f(u)";

  input MetalHydrideStorage.M.State s;

algorithm

  y := MetalHydrideStorage.M.equilibriumPressureAdsorption(MetalHydrideStorage.M.State(s.T, s.p, u)) - s.p;

end MetalHydrideStorage.M.equilibriumLoadingAdsorptionDummy;



function MetalHydrideStorage.M.equilibriumPressureAdsorption

  input MetalHydrideStorage.M.State s;

  output Real p_eq(quantity = "Pressure", unit = "Pa", displayUnit = "bar");

end MetalHydrideStorage.M.equilibriumPressureAdsorption;



function Modelica.Math.Nonlinear.solveOneNonlinearEquation "Solve f(u) = 0 in a very reliable and efficient way (f(u_min) and f(u_max) must have different signs)"

  input f<function>(#Real u) => #Real f "Function y = f(u); u is computed so that y=0";

  input Real u_min "Lower bound of search interval";

  input Real u_max "Upper bound of search interval";

  input Real tolerance = 1e-013 "Relative tolerance of solution u";

  output Real u "Value of independent variable u so that f(u) = 0";

  protected constant Real eps = 1e-015 "machine epsilon";

  protected Real c "Intermediate point a <= c <= b";

  protected Real d;

  protected Real e "b - a";

  protected Real m;

  protected Real s;

  protected Real p;

  protected Real q;

  protected Real r;

  protected Real tol;

  protected Real fa "= f(a)";

  protected Real fb "= f(b)";

  protected Real fc;

  protected Boolean found = false;

  protected Real a = u_min "Current best minimum interval value";

  protected Real b = u_max "Current best maximum interval value";

algorithm

  fa := unbox(f(#(u_min)));

  fb := unbox(f(#(u_max)));

  fc := fb;

  if fa > 0.0 and fb > 0.0 or fa < 0.0 and fb < 0.0 then

    Modelica.Utilities.Streams.error("The arguments u_min and u_max provided in the function call

        solveOneNonlinearEquation(f,u_min,u_max)

    do not bracket the root of the single non-linear equation 0=f(u):

      u_min  = " + String(u_min, 0, true, 6) + "

    " + "  u_max  = " + String(u_max, 0, true, 6) + "

    " + "  fa = f(u_min) = " + String(fa, 0, true, 6) + "

    " + "  fb = f(u_max) = " + String(fb, 0, true, 6) + "

    " + "fa and fb must have opposite sign which is not the case");

  end if;

  c := a;

  fc := fa;

  e := b - a;

  d := e;

  while not found loop

    if abs(fc) < abs(fb) then

      a := b;

      b := c;

      c := a;

      fa := fb;

      fb := fc;

      fc := fa;

    end if;

    tol := 2e-015 * abs(b) + tolerance;

    m := 0.5 * (c - b);

    if abs(m) <= tol or fb == 0.0 then

      found := true;

      u := b;

    else

      if abs(e) < tol or abs(fa) <= abs(fb) then

        e := m;

        d := e;

      else

        s := fb / fa;

        if a == c then

          p := 2.0 * m * s;

          q := 1.0 - s;

        else

          q := fa / fc;

          r := fb / fc;

          p := s * (2.0 * m * q * (q - r) - (b - a) * (-1.0 + r));

          q := (-1.0 + q) * (-1.0 + r) * (-1.0 + s);

        end if;

        if p > 0.0 then

          q := -q;

        else

          p := -p;

        end if;

        s := e;

        e := d;

        if 2.0 * p < 3.0 * m * q - abs(tol * q) and p < abs(0.5 * s * q) then

          d := p / q;

        else

          e := m;

          d := e;

        end if;

      end if;

      a := b;

      fa := fb;

      b := b + (if abs(d) > tol then d else if m > 0.0 then tol else -tol);

      fb := unbox(f(#(b)));

      if fb > 0.0 and fc > 0.0 or fb < 0.0 and fc < 0.0 then

        c := a;

        fc := fa;

        e := b - a;

        d := e;

      end if;

    end if;

  end while;

end Modelica.Math.Nonlinear.solveOneNonlinearEquation;



function Modelica.Utilities.Streams.error "Print error message and cancel all actions"

  input String string "String to be printed to error message window";



  external "C" ModelicaError(string);

end Modelica.Utilities.Streams.error;



class MetalHydrideStorage.M

  Real r = MetalHydrideStorage.M.equilibriumLoadingAdsorption(MetalHydrideStorage.M.State(time, 2.0, 3.0));

end MetalHydrideStorage.M;