"""
Ported by SCL, 27 Oct 2012.
...from combination of predprey.m and predprey_rh.m
The original top-of-file comments appear in function docstrings.
"""

import numpy as np


def predprey_rh(x, t, rnom, predprey_K, predprey_xd, predprey_ud):
    """
    predprey_rh.m - feedback function for predator prey system
    RMM, 13 Sep 06

    This function computes a simple linear control law, with saturation, for
    the predator prey system
    """
    # Compute the feedback control law
    rh = (-np.dot(predprey_K, (x - predprey_xd)) + predprey_ud)[0]

    # Saturate the control input between 0 and 4*rnom
    if rh < 0:
        rh = 0
    if rh > 4*rnom:
        rh = 4*rnom
    return rh


def predprey(x, t,
             r = 1.6, d = 0.56, b = 0.6, k = 125, a = 3.2, c = 50,
             control_args=None):
    """
    predprey.m - updated predator prey model for CDS 101/110
    RMM, Oct 07

    This file contains a simple model for a predator prey system.  The 
    model is developed following the treatment in the book by Murray
    (the other one).

    The states for these equations are

        x(1)	# of rabbits (prey)
        x(2)	# of foxes (predator)

    The equations allow for a quadratic term in both the predator and 
    prey species that cause the population to fall off if you get too 
    many of a given species.  These terms are turned off by default.
    """
    # Extract the states
    x1 = x[0]
    x2 = x[1]

    # Check to see if a feedback law has been defined (when r < 0)
    if r < 0 and (control_args is not None):
      # Assume that the user has defined a function that we can call
      # Add this value to the nominal value of the parameter (renegated)
      r = -r + predprey_rh(x, t, -r, *control_args)

    dx1 = r * x1 * (1 - x1/k) - a * x2 * x1/(c + x1)
    dx2 = b *a *x2 * x1 / (c + x1) - d * x2

    return np.array([dx1, dx2])