def heisenT(beta):
    # Following notation of Xiang (http://arxiv.org/pdf/1201.1144v4.pdf)
    W=np.array([[math.sqrt(math.cosh(beta)), math.sqrt(math.sinh(beta))],
               [math.sqrt(math.cosh(beta)), -math.sqrt(math.sinh(beta))]])
    #T=np.einsum("au,ad,al,ar->udlr",W,W,W,W)
    T=np.einsum("al,ar->lr",W,W)
    return T

 def testHyperbolic(self):
     self.assertEqual(math.sinh(5), hyperbolic.sineh_op(5))
     self.assertEqual(math.cosh(5), hyperbolic.cosineh_op(5))
     self.assertEqual(math.tanh(5), hyperbolic.tangenth_op(5))
     self.assertEqual(1. / math.sinh(5), hyperbolic.cosecanth_op(5))
     self.assertEqual(1. / math.cosh(5), hyperbolic.secanth_op(5))
     self.assertEqual(1. / math.tanh(5), hyperbolic.cotangenth_op(5))

def sinh(x):
    """ Return the hyperbolic sine of x. """

    if math.isinf(x.real) or math.isinf(x.imag):
        # need to raise DomainError if y is +/- infinity and x is not
        # a NaN and not -infinity
        if math.isinf(x.imag) and not math.isnan(x.real):
            raise ValueError("math domain error")

        if math.isinf(x.real) and -_INF < x.imag < _INF and x.imag != .0:
            if x.real > 0:
                _real = math.copysign(_INF, cos(x.imag))
                _imag = math.copysign(_INF, sin(x.imag))
            else:
                _real = -math.copysign(_INF, cos(x.imag))
                _imag = math.copysign(_INF, sin(x.imag))
            return complex(_real, _imag)

        return  _SPECIAL_VALUE(x,_sinh_special_values)

    if math.fabs(x.real) > _CM_LOG_LARGE_DOUBLE:
        x_minus_one = x.real - math.copysign(1.0, x.real)
        _real = math.cos(x.imag)*math.sinh(x.imag)*math.e
        _imag = math.sin(x.imag)*math.cosh(x.imag)*math.e
    else:
        _real = math.cos(x.imag)*math.sinh(x.real)
        _imag = math.sin(x.imag)*math.cosh(x.real)

    if math.isinf(_real) or math.isinf(_imag):
        raise OverflowError()

    return complex(_real, _imag)

def cosh(x):
    """Return the hyperbolic cosine of x."""

    # special treatment for cosh(+/-inf + iy) if y is not a NaN
    if isinf(x):
        if -_INF < x.imag < _INF and x.imag != .0:
            if x.real > 0:
                _real = math.copysign(_INF, math.cos(x.imag))
                _imag = math.copysign(_INF, math.sin(x.imag))
            else:
                _real = math.copysign(_INF, math.cos(x.imag))
                _imag = -math.copysign(_INF, math.sin(x.imag))
            return complex(_real,_imag)
        else:
            # need to raise math domain error if y is +/- infinity and x is not a NaN
            if x.imag != .0 and not math.isnan(x.real):
                raise ValueError("math domain error")
            return _SPECIAL_VALUE(x,_cosh_special_values)

    if math.fabs(x.real) > _CM_LOG_LARGE_DOUBLE:
        #  deal correctly with cases where cosh(x.real) overflows but
        #  cosh(z) does not.
        x_minus_one = x.real - math.copysign(1.0, x.real)
        _real = cos(x.imag) * math.cosh(x_minus_one) * math.e
        _imag = sin(x.imag) * math.sinh(x_minus_one) * math.e
    else:
        _real = math.cos(x.imag) * math.cosh(x.real)
        _imag = math.sin(x.imag) * math.sinh(x.real)

    ret = complex(_real, _imag)
    #  detect overflow
    if isinf(ret):
        raise OverflowError()
    return ret

def to_wgs84(x, y, twd67=False):
    '''
    The formula is based on
    "https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system"
    '''

    if twd67:
        x, y = twd67_to_twd97(x, y)

    x /= 1000.0  # m to km
    y /= 1000.0  # m to km

    xi = (y - N0) / (k0 * A)
    eta = (x - E0) / (k0 * A)
    xip = (xi -
           b1 * math.sin(2 * xi) * math.cosh(2 * eta) -
           b2 * math.sin(4 * xi) * math.cosh(4 * eta) -
           b3 * math.sin(6 * xi) * math.cosh(6 * eta))
    etap = (eta -
            b1 * math.cos(2 * xi) * math.sinh(2 * eta) -
            b2 * math.cos(4 * xi) * math.sinh(4 * eta) -
            b3 * math.cos(6 * xi) * math.sinh(6 * eta))
    chi = math.asin(math.sin(xip) / math.cosh(etap))
    lat = math.degrees(
            chi +
            d1 * math.sin(2 * chi) +
            d2 * math.sin(4 * chi) +
            d3 * math.sin(6 * chi))
    lon = lon0 + math.degrees(
            math.atan(math.sinh(etap) / math.cos(xip)))

    return (lat, lon)

def to_latlon(northings, eastings, altitude):
    f = 1/298.257223563
    n = f / (2 - f)
    n_0 = 0.0
    k_0 = 0.9996
    e_0 = 500.
    a = 6378.137 #km
    big_a = (a / (1 + n)) * (1 + (n**2 / 4) + (n**4 / 64)) #approximately
    alpha_1 = .5 * n - (2. / 3.) * n**2 + (5. / 16.) * n**3
    alpha_2 = (13./48.) * n**2 - (3./5.) * n**3
    alpha_3 = (61./240.) * n**3
    beta_1 = (1./2.) * n - (2. / 3.) * n**2 + (37./96.) * n**3
    beta_2 = (1./48.) * n**2 + (1. / 15.) * n**3
    beta_3 = (17. /480.) * n**3
    delta_1 = 2. * n - (2./3.) * n**2 - 2* n**3
    delta_2 = (7./3.) * n**2 - (8. / 5.) * n**3 
    delta_3 = (56./15.) * n**3

    psi = (n - n_0) / (k_0 * big_a)
    eta = (e - e_0) / (k_0 * big_a)
    psi_prime = psi - ((beta_1 * sin(2. * 1 * eta) * cosh(2. * 1 * eta)) + 
        (beta_2 * sin(2. * 2 * eta) * cosh(2. * 2 * eta)) + (beta_3 * sin(2. * 3 * eta) * cosh(2. * 3 * eta)))
    eta_prime = eta - ((beta_1 * cos(2. * 1 * psi) * sinh(2. * 1 * eta)) + 
        (beta_2 * cos(2. * 2 * psi) * sinh(2. * 2 * eta)) + (beta_3 * cos(2. * 3 * psi) * sinh(2. * 3 * eta)))
    sigma_prime = 1. - ((2. * 1 * beta_1 * cos(2. * 1 * psi) * cosh(2. * 1 * eta)) +
        (2. * 2 * beta_2 * cos(2. * 2 * psi) * cosh(2. * 2 * eta)) + (2. * 3 * beta_3 * cos(2. * 3 * psi) * cosh(2. * 3 * eta)))
    tau_prime = ((2. * 1 * beta_1 * sin(2. * 1 * psi) * sinh(2. * 1 * eta)) +
        (2. * 2 * beta_2 * sin(2. * 2 * psi) * sinh(2. * 2 * eta)) + (2. * 3 * beta_3 * sin(2. * 3 * psi) * sinh(2. * 3 * eta)))
    chi = asin (sin(psi_prime)/cosh(eta_prime)) 
    phi = chi + (delta_1 * sin(2. * 1 * chi)) + (delta_2 * sin(2. * 2 * chi)) + (delta_3 * sin(2. * 3 * chi))
    lambda_0 = 0
    lambdu = 0 
    k =  0
    gamma = 0
    return None

def lalo_to_xy(la, lo, lo0, E0, ellipsoid):
  lo0 = math.radians(lo0)  
  la = math.radians(la)
  lo = math.radians(lo)  
  
  e = ellipsoid['e']
  k0 = ellipsoid['k0']
  h1p = ellipsoid['h1p']
  h2p = ellipsoid['h2p']
  h3p = ellipsoid['h3p']
  h4p = ellipsoid['h4p']
  A1 = ellipsoid['A1']  
    
  Q = asinh(math.tan(la)) - e * atanh(e * math.sin(la))
  be = math.atan(math.sinh(Q))
  nnp = atanh(math.cos(be) * math.sin(lo - lo0))
  Ep = math.asin(math.sin(be) * math.cosh(nnp))  
  E1 = h1p * math.sin(2.0 * Ep) * math.cosh(2.0 * nnp)
  E2 = h2p * math.sin(4.0 * Ep) * math.cosh(4.0 * nnp)
  E3 = h3p * math.sin(6.0 * Ep) * math.cosh(6.0 * nnp)
  E4 = h4p * math.sin(8.0 * Ep) * math.cosh(8.0 * nnp)
  nn1 = h1p * math.cos(2.0 * Ep) * math.sinh(2.0 * nnp)
  nn2 = h2p * math.cos(4.0 * Ep) * math.sinh(4.0 * nnp)
  nn3 = h3p * math.cos(6.0 * Ep) * math.sinh(6.0 * nnp)
  nn4 = h4p * math.cos(8.0 * Ep) * math.sinh(8.0 * nnp)
  E = Ep + E1 + E2 + E3 + E4
  nn = nnp + nn1 + nn2 + nn3 + nn4
  
  XY = {}
  XY['N'] = A1 * E * k0
  XY['E'] = A1 * nn * k0 + E0

  return XY

def cosh(x):
    _cosh_special = [
        [inf+nanj, None, inf, complex(float("inf"), -0.0), None, inf+nanj, inf+nanj],
        [nan+nanj, None, None, None, None, nan+nanj, nan+nanj],
        [nan, None, 1, complex(1, -0.0), None, nan, nan],
        [nan, None, complex(1, -0.0), 1, None, nan, nan],
        [nan+nanj, None, None, None, None, nan+nanj, nan+nanj],
        [inf+nanj, None, complex(float("inf"), -0.0), inf, None, inf+nanj, inf+nanj],
        [nan+nanj, nan+nanj, nan, nan, nan+nanj, nan+nanj, nan+nanj]
    ]

    z = _make_complex(x)

    if not isfinite(z):
        if math.isinf(z.imag) and not math.isnan(z.real):
            raise ValueError
        if math.isinf(z.real) and math.isfinite(z.imag) and z.imag != 0:
            if z.real > 0:
                return complex(math.copysign(inf, math.cos(z.imag)),
                               math.copysign(inf, math.sin(z.imag)))
            return complex(math.copysign(inf, math.cos(z.imag)),
                           -math.copysign(inf, math.sin(z.imag)))
        return _cosh_special[_special_type(z.real)][_special_type(z.imag)]

    if abs(z.real) > _LOG_LARGE_DOUBLE:
        x_minus_one = z.real - math.copysign(1, z.real)
        ret = complex(e * math.cos(z.imag) * math.cosh(x_minus_one),
                      e * math.sin(z.imag) * math.sinh(x_minus_one))
    else:
        ret = complex(math.cos(z.imag) * math.cosh(z.real),
                      math.sin(z.imag) * math.sinh(z.real))
    if math.isinf(ret.real) or math.isinf(ret.imag):
        raise OverflowError

    return ret

def from_utm(easting: float, northing: float,
             zone: int, hemisphere: int =1) -> (float, float):
    """
    Convert UTM coordinates to decimal latitude and longitude coordinates.

    Keyword Arguments:
        easting: The easting in m.
        northing: The northing in m.
        zone: The zone, an integer between 1 and 60, inclusive.
        hemisphere: A signed number, where a negative number indicates the
                    coordinates are located in the southern hemisphere.

    Returns:
        latitude: The latitude in decimal degrees.
        longitude: The longitude in deciaml degrees.

    Raises:
        OverflowError: The coordinate does not exist.
    """
    easting, northing = easting/1000, northing/1000
    northing_ = 10000 if hemisphere < 0 else 0

    xi_ = xi = (northing - northing_)/(k0*A)
    eta_ = eta = (easting - easting_)/(k0*A)
    for j in range(1, 4):
        p, q = 2*j*xi, 2*j*eta
        xi_ -= beta[j - 1]*sin(p)*cosh(q)
        eta_ -= beta[j - 1]*cos(p)*sinh(q)

    chi = asin(sin(xi_)/cosh(eta_))
    latitude = chi + sum(delta[j - 1]*sin(2*j*chi) for j in range(1, 4))
    longitude_ = radians(6*zone - 183)
    longitude = longitude_ + atan2(sinh(eta_), cos(xi_))
    return degrees(latitude), degrees(longitude)

def sinh(x):

    _sinh_special = [
        [inf+nanj, None, complex(-float("inf"), -0.0), -inf, None, inf+nanj, inf+nanj],
        [nan+nanj, None, None, None, None, nan+nanj, nan+nanj],
        [nanj, None, complex(-0.0, -0.0), complex(-0.0, 0.0), None, nanj, nanj],
        [nanj, None, complex(0.0, -0.0), complex(0.0, 0.0), None, nanj, nanj],
        [nan+nanj, None, None, None, None, nan+nanj, nan+nanj],
        [inf+nanj, None, complex(float("inf"), -0.0), inf, None, inf+nanj, inf+nanj],
        [nan+nanj, nan+nanj, complex(float("nan"), -0.0), nan, nan+nanj, nan+nanj, nan+nanj]
    ]

    z = _make_complex(x)

    if not isfinite(z):
        if math.isinf(z.imag) and not math.isnan(z.real):
            raise ValueError
        if math.isinf(z.real) and math.isfinite(z.imag) and z.imag != 0:
            if z.real > 0:
                return complex(math.copysign(inf, math.cos(z.imag)),
                               math.copysign(inf, math.sin(z.imag)))
            return complex(-math.copysign(inf, math.cos(z.imag)),
                           math.copysign(inf, math.sin(z.imag)))
        return _sinh_special[_special_type(z.real)][_special_type(z.imag)]

    if abs(z.real) > _LOG_LARGE_DOUBLE:
        x_minus_one = z.real - math.copysign(1, z.real)
        return complex(math.cos(z.imag) * math.sinh(x_minus_one) * e,
                       math.sin(z.imag) * math.cosh(x_minus_one) * e)
    return complex(math.cos(z.imag) * math.sinh(z.real),
                   math.sin(z.imag) * math.cosh(z.real))

def c_cosh(x, y):
    if not isfinite(x) or not isfinite(y):
        if isinf(x) and isfinite(y) and y != 0.:
            if x > 0:
                real = copysign(INF, math.cos(y))
                imag = copysign(INF, math.sin(y))
            else:
                real = copysign(INF, math.cos(y))
                imag = -copysign(INF, math.sin(y))
            r = (real, imag)
        else:
            r = cosh_special_values[special_type(x)][special_type(y)]

        # need to raise ValueError if y is +/- infinity and x is not
        # a NaN
        if isinf(y) and not isnan(x):
            raise ValueError("math domain error")
        return r

    if fabs(x) > CM_LOG_LARGE_DOUBLE:
        # deal correctly with cases where cosh(x) overflows but
        # cosh(z) does not.
        x_minus_one = x - copysign(1., x)
        real = math.cos(y) * math.cosh(x_minus_one) * math.e
        imag = math.sin(y) * math.sinh(x_minus_one) * math.e
    else:
        real = math.cos(y) * math.cosh(x)
        imag = math.sin(y) * math.sinh(x)
    if isinf(real) or isinf(imag):
        raise OverflowError("math range error")
    return real, imag

def c_sinh(x, y):
    # special treatment for sinh(+/-inf + iy) if y is finite and nonzero
    if not isfinite(x) or not isfinite(y):
        if isinf(x) and isfinite(y) and y != 0.:
            if x > 0:
                real = copysign(INF, math.cos(y))
                imag = copysign(INF, math.sin(y))
            else:
                real = -copysign(INF, math.cos(y))
                imag = copysign(INF, math.sin(y))
            r = (real, imag)
        else:
            r = sinh_special_values[special_type(x)][special_type(y)]

        # need to raise ValueError if y is +/- infinity and x is not
        # a NaN
        if isinf(y) and not isnan(x):
            raise ValueError("math domain error")
        return r

    if fabs(x) > CM_LOG_LARGE_DOUBLE:
        x_minus_one = x - copysign(1., x)
        real = math.cos(y) * math.sinh(x_minus_one) * math.e
        imag = math.sin(y) * math.cosh(x_minus_one) * math.e
    else:
        real = math.cos(y) * math.sinh(x)
        imag = math.sin(y) * math.cosh(x)
    if isinf(real) or isinf(imag):
        raise OverflowError("math range error")
    return real, imag

def recN(K, MAX_VALUE):
    if(abs(K[0] + K[1] + K[2] + K[3] + K[4] + K[5] + K[6] + K[7]) > MAX_VALUE
    	or abs(K[0] + K[1] + K[2] + K[3] - K[4] - K[5] - K[6] - K[7]) > MAX_VALUE):
    	return recN0T(K)
    else:
    	return .5 * math.log(math.cosh(K[0] + K[1] + K[2] + K[3] + K[4] + K[5] + K[6] + K[7])
               			  / math.cosh(K[0] + K[1] + K[2] + K[3] - K[4] - K[5] - K[6] - K[7]))

def tri_angles(L):
  ans = [0, 0, 0]
  for i in xrange(3):
    top = math.cosh(L[i])*math.cosh(L[(i+2)%3]) - math.cosh( L[(i+1)%3] )
    bottom = math.sinh(L[i])*math.sinh(L[(i+2)%3])
    ans[i] = math.acos(top/bottom)
  return ans

def cosh(self, other=None):
# Return hyperbolic cosine of interval

	if other != None:
		intv = IReal(self, other)
	else:
		if type(self) == float or type(self) == str:
			intv = IReal(self)
		else:
			intv = self
 
	if math.cosh(intv.inf) > math.cosh(intv.sup):
		inf = max(intv.inf, intv.sup)
		sup = min(intv.inf, intv.sup)
	else:
		inf = intv.inf
		sup = intv.sup

	if 0 in intv:
		inf = 0

	ireal.rounding.set_mode(1)

	ireal.rounding.set_mode(-1)
	ireal.rounding.set_mode(-1)

	new_inf = math.cosh(inf)
	ireal.rounding.set_mode(1)
	new_sup = max(float(IReal('%.16f' % math.cosh(sup)).sup), float(IReal('%.20f' % math.cosh(sup)).sup))

	return IReal(new_inf, new_sup)

  def test_funcs_multi(self):
    pi = math.pi

    # sin family
    self.assertQuantity(data.sin(Quantity( (0,pi/2), (2,2))), (0,1), (2,0)) 
    self.assertQuantity(data.sinh(Quantity((0,1), (2,2))), (0, math.sinh(1)), (2, math.cosh(1)*2)) 
    self.assertQuantity(data.asin(Quantity((0,0.5), (2,2))), (0,math.asin(0.5)), (2,2/math.sqrt(1-0.5**2))) 
    self.assertQuantity(data.asinh(Quantity((0,1), (2,2))), (0,math.asinh(1)), (2,2/math.sqrt(1+1**2))) 

    # cos family
    self.assertQuantity(data.cos(Quantity((0,pi/2), (2,2))), (1,0), (0,-2)) 
    self.assertQuantity(data.cosh(Quantity((0,1), (2,2))), (1,math.cosh(1)), (0,math.sinh(1)*2)) 
    self.assertQuantity(data.acos(Quantity((0,0.5), (2,2))), (math.acos(0),math.acos(0.5)), (-2,-2/math.sqrt(1-0.5**2)))
    self.assertQuantity(data.acosh(Quantity((2,3), (2,2))), (math.acosh(2), math.acosh(3)), (2/math.sqrt(2**2-1),2/math.sqrt(3**2-1)))

    # tan family
    self.assertQuantity(data.tan(Quantity((0,1), (2,2))), (0,math.tan(1)), (2,2/math.cos(1)**2))
    self.assertQuantity(data.tanh(Quantity((0,1), (2,2))), (0,math.tanh(1)), (2,2/math.cosh(1)**2)) 
    self.assertQuantity(data.atan(Quantity((0,1), (2,2))), (0, math.atan(1)), (2,2/(1+1**2))) 
    self.assertQuantity(data.atan2(Quantity((0,1), (2,2)), Quantity((1,1), (0,0))), (0,math.atan(1)), (2,2/(1+1**2))) 
    self.assertQuantity(data.atanh(Quantity((0,0.5), (2,2))), (0,math.atanh(0.5)), (2,2/(1-0.5**2)))

    #misc
    self.assertQuantity(data.sqrt(Quantity((1,4), (2,2))), (1,2), (1,1/2))
    self.assertQuantity(data.exp(Quantity((1,4), (2,2))), (math.e, math.e**4), (2 * math.e,2*math.e**4))
    self.assertQuantity(data.log(Quantity((1,4), (2,2))), (0, math.log(4)), (2,1/2))

def fromwgs84(lat, lng, pkm=False):
    """
    Convert coordintes from WGS84 to TWD97

    pkm true for Penghu, Kinmen and Matsu area
    The latitude and longitude can be in the following formats:
        [+/-]DDD°MMM'SSS.SSSS" (unicode)
        [+/-]DDD°MMM.MMMM' (unicode)
        [+/-]DDD.DDDDD (string, unicode or float)
    The returned coordinates are in meters
    """

    _lng0 = lng0pkm if pkm else lng0

    lat = radians(todegdec(lat))
    lng = radians(todegdec(lng))

    t = sinh((atanh(sin(lat)) - 2*pow(n,0.5)/(1+n)*atanh(2*pow(n,0.5)/(1+n)*sin(lat))))
    epsilonp = atan(t/cos(lng-_lng0))
    etap = atan(sin(lng-_lng0) / pow(1+t*t, 0.5))

    E = E0 + k0*A*(etap + alpha1*cos(2*1*epsilonp)*sinh(2*1*etap) + 
                          alpha2*cos(2*2*epsilonp)*sinh(2*2*etap) +
                          alpha3*cos(2*3*epsilonp)*sinh(2*3*etap))
    N = N0 + k0*A*(epsilonp + alpha1*sin(2*1*epsilonp)*cosh(2*1*etap) +
                              alpha2*sin(2*2*epsilonp)*cosh(2*2*etap) +
                              alpha3*sin(2*3*epsilonp)*cosh(2*3*etap))

    return E*1000, N*1000

    def gridToGeodetic(self, north, east):
        """
        Transformation from grid coordinates to geodetic coordinates.
        @param north (corresponds to X in RT 90 and N in SWEREF 99.)
        @param east (corresponds to Y in RT 90 and E in SWEREF 99.)
        @return (latitude, longitude)
        """
        if (self._initialized == False):
            return None

        deg_to_rad = math.pi / 180
        lambda_zero = self._central_meridian * deg_to_rad
        xi = (north - self._false_northing) / (self._scale * self._a_roof)        
        eta = (east - self._false_easting) / (self._scale * self._a_roof)
        xi_prim = xi - \
                self._delta1*math.sin(2.0*xi) * math.cosh(2.0*eta) - \
                self._delta2*math.sin(4.0*xi) * math.cosh(4.0*eta) - \
                self._delta3*math.sin(6.0*xi) * math.cosh(6.0*eta) - \
                self._delta4*math.sin(8.0*xi) * math.cosh(8.0*eta)
        eta_prim = eta - \
                self._delta1*math.cos(2.0*xi) * math.sinh(2.0*eta) - \
                self._delta2*math.cos(4.0*xi) * math.sinh(4.0*eta) - \
                self._delta3*math.cos(6.0*xi) * math.sinh(6.0*eta) - \
                self._delta4*math.cos(8.0*xi) * math.sinh(8.0*eta)
        phi_star = math.asin(math.sin(xi_prim) / math.cosh(eta_prim))
        delta_lambda = math.atan(math.sinh(eta_prim) / math.cos(xi_prim))
        lon_radian = lambda_zero + delta_lambda
        lat_radian = phi_star + math.sin(phi_star) * math.cos(phi_star) * \
                (self._Astar + \
                 self._Bstar*math.pow(math.sin(phi_star), 2) + \
                 self._Cstar*math.pow(math.sin(phi_star), 4) + \
                 self._Dstar*math.pow(math.sin(phi_star), 6))      
        lat = lat_radian * 180.0 / math.pi
        lon = lon_radian * 180.0 / math.pi
        return (lat, lon)

def gal_weights(Z, R):
    '''
    Returns a weight based on a particular model of the MW.
    For now we will use a two-disk model with the form below.
    This can be expanded at a later time.

    Z - Height above/below galactic plane
    R - Distance from galactic center

    '''

    # Parameters
    thick_s_height = 0.674
    thick_s_length = 2.51
    thin_s_height = 0.233
    thin_s_length = 2.34
    a = 0.1


    weight = ( ( ( math.cosh(Z / 2 / thin_s_height) ) ** (-2) )
        * math.exp(-R / thin_s_length) +
        a * ( ( math.cosh(Z / 2 / thick_s_height) ) ** (-2) )
        * math.exp(-R / thick_s_length))

    return weight

def df(y,t): 
    Vs=y[0]
    ns=y[1]
    hs=y[2]
    
    # SOMATIC FUNCTIONS
    minfs=1/(1+math.exp((Vs-vm) /sm))
    ninfs=1/(1+math.exp((Vs-vn) /sn))
    minfps=1/(1+math.exp((Vs-vmp)/smp))
    hinfs=1/(1+math.exp((Vs-vh) /sh))
    
    tauns=taunb/math.cosh((Vs-vn)/(2*sn))
    tauhs=tauhb/math.cosh((Vs-vh)/(2*sh))
    
    I_nas=gna*math.pow(minfs,3)*(1-ns)*(Vs-vna)
    I_ks=gk*math.pow(ns,4)*(Vs-Vk)
    I_naps=gnaps*minfps*hs*(Vs-vna)
    
    I_ls =gls*(Vs-vleaks)

    # SOMATIC EQUATIONS

    dVs= (-I_ks - I_nas-I_naps-I_ls-Iaps)/Cms
    dns= (ninfs-ns)/tauns
    dhs= (hinfs-hs)/tauhs
    
    return [dVs,dns,dhs]