http://www.pisces-conservation.com/growthhelp/janoschek.htm
    
    Parameters
    ----------
    x: int
    a: float
    beta: float
    k: float
    delta: float

    Returns
    -------
    float
        a - (a - beta) * np.exp(-k*x**delta)
506656
def janoschek(x, a, beta, k, delta):
    
    return a - (a - beta) * np.exp(-k*x**delta)
636815
def actDerivASIG(self, x):
        
        def act(v):
            if   v < -15.0: return 0.0
            elif v >  15.0: return 1.0
            else: return 1.0 / (1.0 + Numeric.exp(-v))
        return (act(x) * (1.0 - act(x))) + self.sigmoid_prime_offset
391467
def feed(self, x):
        
        self._n += 1
        delta = x - self._mean
        self._mean += delta * (1.0 / self._n)
        delta2 = x - self._mean
        self._M2 += delta * delta2
1072159
def _F(self, X):
        
        if isinstance(X, int) or isinstance(X, float):
            if X < 1 and X > 0:
                a = 1/(X**2-1)*(1-2/np.sqrt(1-X**2)*np.arctanh(np.sqrt((1-X)/(1+X))))
            elif X == 1:
                a = 1./3
            elif X > 1:
                a = 1/(X**2-1)*(1-2/np.sqrt(X**2-1)*np.arctan(np.sqrt((X-1)/(1+X))))
            else:  # X == 0:
                c = 0.0000001
                a = 1/(-1)*(1-2/np.sqrt(1)*np.arctanh(np.sqrt((1-c)/(1+c))))

        else:
            a = np.empty_like(X)
            x = X[(X < 1) & (X > 0)]
            a[(X < 1) & (X > 0)] = 1/(x**2-1)*(1-2/np.sqrt(1-x**2)*np.arctanh(np.sqrt((1-x)/(1+x))))

            a[X == 1] = 1./3.

            x = X[X > 1]
            a[X > 1] = 1/(x**2-1)*(1-2/np.sqrt(x**2-1)*np.arctan(np.sqrt((x-1)/(1+x))))
            # a[X>y] = 0

            c = 0.0000001
            a[X == 0] = 1/(-1)*(1-2/np.sqrt(1)*np.arctanh(np.sqrt((1-c)/(1+c))))
        return a
134293
def Sine(x, a, omega, phi, y0):
    
    return a * np.sin(x * omega + phi) + y0
273126
def exp_comp_(t, alpha=1, beta=1):
    """"""
    return beta * (1 - np.exp(-alpha * np.maximum(0, t - 10 * beta)))
1083599
def _Fgamma(self, x, Ep):
        
        L = np.log(Ep)
        B = 1.30 + 0.14 * L + 0.011 * L ** 2  # Eq59
        beta = (1.79 + 0.11 * L + 0.008 * L ** 2) ** -1  # Eq60
        k = (0.801 + 0.049 * L + 0.014 * L ** 2) ** -1  # Eq61
        xb = x ** beta

        F1 = B * (np.log(x) / x) * ((1 - xb) / (1 + k * xb * (1 - xb))) ** 4
        F2 = (
            1.0 / np.log(x)
            - (4 * beta * xb) / (1 - xb)
            - (4 * k * beta * xb * (1 - 2 * xb)) / (1 + k * xb * (1 - xb))
        )

        return F1 * F2
652366
def _B(self, x, a, b):
        
        return special.betainc(a, b, x) * special.beta(a, b)
646314
def bbox(self):
        """"""
        # a(t) = radians(self.theta + self.delta*t)
        #      = (2*pi/360)*(self.theta + self.delta*t)
        # x'=0: ~~~~~~~~~
        # -rx*cos(phi)*sin(a(t)) = ry*sin(phi)*cos(a(t))
        # -(rx/ry)*cot(phi)*tan(a(t)) = 1
        # a(t) = arctan(-(ry/rx)tan(phi)) + pi*k === atan_x
        # y'=0: ~~~~~~~~~~
        # rx*sin(phi)*sin(a(t)) = ry*cos(phi)*cos(a(t))
        # (rx/ry)*tan(phi)*tan(a(t)) = 1
        # a(t) = arctan((ry/rx)*cot(phi))
        # atanres = arctan((ry/rx)*cot(phi)) === atan_y
        # ~~~~~~~~
        # (2*pi/360)*(self.theta + self.delta*t) = atanres + pi*k
        # Therfore, for both x' and y', we have...
        # t = ((atan_{x/y} + pi*k)*(360/(2*pi)) - self.theta)/self.delta
        # for all k s.t. 0 < t < 1
        from math import atan, tan

        if cos(self.phi) == 0:
            atan_x = pi/2
            atan_y = 0
        elif sin(self.phi) == 0:
            atan_x = 0
            atan_y = pi/2
        else:
            rx, ry = self.radius.real, self.radius.imag
            atan_x = atan(-(ry/rx)*tan(self.phi))
            atan_y = atan((ry/rx)/tan(self.phi))

        def angle_inv(ang, k):  # inverse of angle from Arc.derivative()
            return ((ang + pi*k)*(360/(2*pi)) - self.theta)/self.delta

        xtrema = [self.start.real, self.end.real]
        ytrema = [self.start.imag, self.end.imag]

        for k in range(-4, 5):
            tx = angle_inv(atan_x, k)
            ty = angle_inv(atan_y, k)
            if 0 <= tx <= 1:
                xtrema.append(self.point(tx).real)
            if 0 <= ty <= 1:
                ytrema.append(self.point(ty).imag)
        xmin = max(xtrema)
        return min(xtrema), max(xtrema), min(ytrema), max(ytrema)
349443
protected static double regularizedIncBetaCF(double alpha, double beta, double x) {
    final double FPMIN = Double.MIN_VALUE / NUM_PRECISION;
    double qab = alpha + beta;
    double qap = alpha + 1.0;
    double qam = alpha - 1.0;
    double c = 1.0;
    double d = 1.0 - qab * x / qap;
    if(Math.abs(d) < FPMIN) {
      d = FPMIN;
    }
    d = 1.0 / d;
    double h = d;
    for(int m = 1; m < 10000; m++) {
      int m2 = 2 * m;
      double aa = m * (beta - m) * x / ((qam + m2) * (alpha + m2));
      d = 1.0 + aa * d;
      if(Math.abs(d) < FPMIN) {
        d = FPMIN;
      }
      c = 1.0 + aa / c;
      if(Math.abs(c) < FPMIN) {
        c = FPMIN;
      }
      d = 1.0 / d;
      h *= d * c;
      aa = -(alpha + m) * (qab + m) * x / ((alpha + m2) * (qap + m2));
      d = 1.0 + aa * d;
      if(Math.abs(d) < FPMIN) {
        d = FPMIN;
      }
      c = 1.0 + aa / c;
      if(Math.abs(c) < FPMIN) {
        c = FPMIN;
      }
      d = 1.0 / d;
      double del = d * c;
      h *= del;
      if(Math.abs(del - 1.0) <= NUM_PRECISION) {
        break;
      }
    }
    return h;
  }
817142
def _tilt(omega, b, psi):
    
    return np.cosh(psi/2.0)**b * np.exp(-psi**2/2.0 * omega)
296107
def binomialInterval(n, k, alpha = 0.68):
    
    e = float(k)/n
    delta_e = 1/float(n) * np.sqrt(e * (1 - e) * float(n)) * alpha/0.68
    return e, [e - delta_e, e + delta_e]
793846
def fisher(x,k):
    
    return k/(2*np.sinh(k)) * np.exp(k*np.cos(x))*np.sin(x)
1310267
def _calcOnAxisFactor(self, aLocation, deltaAxis, deltasOnSameAxis, deltaLocation):
        
        if deltaAxis == "origin":
            f = 0
            v = 0
        else:
            f = aLocation[deltaAxis]
            v = deltaLocation[deltaAxis]
        i = []
        iv = {}
        for value in deltasOnSameAxis:
            iv[Location(value)[deltaAxis]]=1
        i = sorted(iv.keys())
        r = 0
        B, M, A = [], [], []
        mA, mB, mM = None, None, None
        for value in i:
            if value < f: B.append(value)
            elif value > f: A.append(value)
            else: M.append(value)
        if len(B) > 0:
            mB = max(B)
            B.sort()
        if len(A) > 0:
            mA = min(A)
            A.sort()
        if len(M) > 0:
            mM = min(M)
            M.sort()
        if mM is not None:
            if ((f-_EPSILON <  v) and (f+_EPSILON > v)) or f==v: r = 1
            else: r = 0
        elif mB is not None and mA is not None:
            if v < mB or v > mA: r = 0
            else:
                if v == mA:
                    r = float(f-mB)/(mA-mB)
                else:
                    r = float(f-mA)/(mB-mA)
        elif mB is None and mA is not None:
            if v==A[1]:
                r = float(f-A[0])/(A[1]-A[0])
            elif v == A[0]:
                r = float(f-A[1])/(A[0]-A[1])
            else:
                r = 0
        elif  mB is not None and mA is None:
            if v == B[-2]:
                r = float(f-B[-1])/(B[-2]-B[-1])
            elif v == mB:
                r = float(f-B[-2])/(B[-1]-B[-2])
            else:
                r = 0
        return r
829951
def psit(t, xp, x):
    
    if t == 0:
        return -0.5 / sigma0**2 + \
            (0.5 * (1. - phi0**2) / sigma0**4) * (x - mu0)**2
    else:
        return -0.5 / sigma0**2 + (0.5 / sigma0**4) * \
            ((x - mu0) - phi0 * (xp - mu0))**2
1252067
def get_delta(D, k):
    ''''''
    if k < 0:
        raise Exception('k must be at least 0th order.')
    result = D
    for i in range(k):
        result = D.T.dot(result) if i % 2 == 0 else D.dot(result)
    return result
373458
def calculate_deltat(year, month):
    

    plw = 'Deltat is unknown for years before -1999 and after 3000. ' \
          'Delta values will be calculated, but the calculations ' \
          'are not intended to be used for these years.'

    try:
        if np.any((year > 3000) | (year < -1999)):
            warnings.warn(plw)
    except ValueError:
        if (year > 3000) | (year < -1999):
            warnings.warn(plw)
    except TypeError:
        return 0

    y = year + (month - 0.5)/12

    deltat = np.where(year < -500,

                      -20+32*((y-1820)/100)**2, 0)

    deltat = np.where((-500 <= year) & (year < 500),

                      10583.6-1014.41*(y/100)
                      + 33.78311*(y/100)**2
                      - 5.952053*(y/100)**3
                      - 0.1798452*(y/100)**4
                      + 0.022174192*(y/100)**5
                      + 0.0090316521*(y/100)**6, deltat)

    deltat = np.where((500 <= year) & (year < 1600),

                      1574.2-556.01*((y-1000)/100)
                      + 71.23472*((y-1000)/100)**2
                      + 0.319781*((y-1000)/100)**3
                      - 0.8503463*((y-1000)/100)**4
                      - 0.005050998*((y-1000)/100)**5
                      + 0.0083572073*((y-1000)/100)**6, deltat)

    deltat = np.where((1600 <= year) & (year < 1700),

                      120-0.9808*(y-1600)
                      - 0.01532*(y-1600)**2
                      + (y-1600)**3/7129, deltat)

    deltat = np.where((1700 <= year) & (year < 1800),

                      8.83+0.1603*(y-1700)
                      - 0.0059285*(y-1700)**2
                      + 0.00013336*(y-1700)**3
                      - (y-1700)**4/1174000, deltat)

    deltat = np.where((1800 <= year) & (year < 1860),

                      13.72-0.332447*(y-1800)
                      + 0.0068612*(y-1800)**2
                      + 0.0041116*(y-1800)**3
                      - 0.00037436*(y-1800)**4
                      + 0.0000121272*(y-1800)**5
                      - 0.0000001699*(y-1800)**6
                      + 0.000000000875*(y-1800)**7, deltat)

    deltat = np.where((1860 <= year) & (year < 1900),

                      7.62+0.5737*(y-1860)
                      - 0.251754*(y-1860)**2
                      + 0.01680668*(y-1860)**3
                      - 0.0004473624*(y-1860)**4
                      + (y-1860)**5/233174, deltat)

    deltat = np.where((1900 <= year) & (year < 1920),

                      -2.79+1.494119*(y-1900)
                      - 0.0598939*(y-1900)**2
                      + 0.0061966*(y-1900)**3
                      - 0.000197*(y-1900)**4, deltat)

    deltat = np.where((1920 <= year) & (year < 1941),

                      21.20+0.84493*(y-1920)
                      - 0.076100*(y-1920)**2
                      + 0.0020936*(y-1920)**3, deltat)

    deltat = np.where((1941 <= year) & (year < 1961),

                      29.07+0.407*(y-1950)
                      - (y-1950)**2/233
                      + (y-1950)**3/2547, deltat)

    deltat = np.where((1961 <= year) & (year < 1986),

                      45.45+1.067*(y-1975)
                      - (y-1975)**2/260
                      - (y-1975)**3/718, deltat)

    deltat = np.where((1986 <= year) & (year < 2005),

                      63.86+0.3345*(y-2000)
                      - 0.060374*(y-2000)**2
                      + 0.0017275*(y-2000)**3
                      + 0.000651814*(y-2000)**4
                      + 0.00002373599*(y-2000)**5, deltat)

    deltat = np.where((2005 <= year) & (year < 2050),

                      62.92+0.32217*(y-2000)
                      + 0.005589*(y-2000)**2, deltat)

    deltat = np.where((2050 <= year) & (year < 2150),

                      -20+32*((y-1820)/100)**2
                      - 0.5628*(2150-y), deltat)

    deltat = np.where(year >= 2150,

                      -20+32*((y-1820)/100)**2, deltat)

    deltat = deltat.item() if np.isscalar(year) & np.isscalar(month)\
        else deltat

    return deltat