def dec2sexstr(deci, sfigs=1, hd='h', lead_psign=False):
    """
    Convert decimal degree to a sexagesimal string of the form 'HH:MM:SS.S' by
    default.

    Parameters
    ----------
    sfigs : number
        Number of significant figures after decimal in seconds
    hd : string, ('h', 'd')
        Hour or degree convention
    lead_sign : Boolean
        Whether to include the leading sign +/- in string
    """
    lead_psymb = ''
    dform = {'h': 24 / 360., 'd': 1.}
    (h_frac, h) = _np.modf(deci * dform[hd])
    (m_frac, m) = _np.modf(h_frac * 60)
    s = m_frac * 60
    if (lead_psign is True) & (h >= 0):
        lead_psymb = '+'
    elif (lead_psign is True) & (h < 0):
        lead_psymb = '-'
        h = abs(h)
    coord_string = "{0}{1:0>2d}:{2:0>2d}:{3:0>2d}.{4:0>{sfigs}d}".format(
        lead_psymb, int(h), _np.abs(int(m)), _np.abs(int(s)),
        _np.abs(int(_np.mod(s,1) * 10**sfigs)), sfigs=sfigs)
    return coord_string

def ImSec(imdata, xcntr, ycntr, ex_rad):
    """Find the image section for asymmetry"""
    SizeY, SizeX = imdata.shape
    ExceedSize = 0
    #All are floor to make the size even number
    xmin = np.floor(xcntr - ex_rad)
    ymin = np.floor(ycntr - ex_rad)
    xmax = np.floor(xcntr + ex_rad)
    ymax = np.floor(ycntr + ex_rad)
    if xmin < 0:
        xmin = 0
        cut_xcntr = xcntr
        ExceedSize = 1
    else:
        cut_xcntr = ex_rad + np.modf(xcntr)[0]
    if ymin < 0:
        cut_ycntr = ycntr
        ymin = 0
        ExceedSize = 1
    else:
        cut_ycntr = ex_rad + np.modf(ycntr)[0]
    if xmax > SizeX - 1:
        xmax = SizeX
        ExceedSize = 1
    if ymax > SizeY - 1:
        ymax = SizeY
        ExceedSize = 1
    CutImDa = imdata[ymin:ymax, xmin:xmax].copy()
    SizeY, SizeX = CutImDa.shape
    return CutImDa, cut_xcntr, cut_ycntr, SizeY, SizeX, ymin, \
           ymax, xmin, xmax, ExceedSize

def decimal_to_sexagesimal(decimal):
    '''
    convert decimal hours or degrees to sexagesimal

    Parameters
    ----------
    decimal | float: decimal number to be converted to sexagismal

    Returns
    -------
    tuple:
        int: hours
        int: minutes
        float: seconds
    OR
    tuple:
        int: degrees
        int: arcminutes
        float: arcseconds

    >>> decimal_to_sexagesimal(10.1)
    (10, 5, 59.999999999998721)
    '''
    fractional, integral = np.modf(decimal)
    min_fractional, minutes = np.modf(fractional * 60)
    seconds = min_fractional * 60.

    return integral.astype(int), minutes.astype(int), seconds

def sexa(value):
    sign = np.sign(value, subok=False)
    absolute = np.absolute(value, subok=False)
    fraction, whole = np.modf(absolute)
    fraction, minutes = np.modf(fraction * 60.0)
    seconds = fraction * 60.0
    return sign, whole, minutes, seconds

def jd2gcal(JD=2443716):
    '''From Jean Meeus' Astronomical Algorithms. It starts at 0000.
    '''
    F, Z = np.modf(JD + 0.5)
    
    if Z >= 2299161:
        alpha = (Z - 1867216.25) // 36524.25
        A = Z + 1 + alpha - alpha // 4
    else:
        A = Z

    B = A + 1524
    C = (B - 122.1) // 365.25
    D = np.modf(365.25 * C)[1]
    E = (B - D) // 30.6001
    D = B - D - np.modf(30.6001 * E)[1] + F
    fD, D = np.modf(D)    
    if E < 14:
        M = E - 1
    else:
        M = E - 13
    if M > 2:
        Y = C - 4716
    else:
        Y = C - 4715
    return int(Y), int(M), int(D), fD

def degreDeci2DegreMinut(degdeci=0.0):
    minutes, degres = np.modf(degdeci)
    minutes = abs(minutes) * 60
    secondes, minutes = np.modf(minutes)
    secondes = secondes * 60
    print ("%i° %i' %.2f'' " % (int(degres), int(minutes), secondes))
    return degres, minutes, secondes

 def test_one_argument_two_output_ufunc_inplace(self, value):
     v = 100. * value * u.cm / u.m
     v_copy = v.copy()
     tmp = v.copy()
     check = v
     np.modf(v, tmp, v)  # cannot use out1,out2 keywords with numpy 1.7
     assert check is v
     assert check.unit == u.dimensionless_unscaled
     v2 = v_copy.to(u.dimensionless_unscaled)
     check2 = v2
     np.modf(v2, tmp, v2)
     assert check2 is v2
     assert check2.unit == u.dimensionless_unscaled
     # can also replace in last position if no scaling is needed
     v3 = v_copy.to(u.dimensionless_unscaled)
     check3 = v3
     np.modf(v3, v3, tmp)
     assert check3 is v3
     assert check3.unit == u.dimensionless_unscaled
     # in np<1.13, without __array_ufunc__, one cannot replace input with
     # first output when scaling
     v4 = v_copy.copy()
     if NUMPY_LT_1_13:
         with pytest.raises(TypeError):
             np.modf(v4, v4, tmp)
     else:
         check4 = v4
         np.modf(v4, v4, tmp)
         assert check4 is v4
         assert check4.unit == u.dimensionless_unscaled

def _make_inverse_warp(from_points, to_points, output_region, approximate_grid):
    x_min, y_min, x_max, y_max = output_region
    if approximate_grid is None: approximate_grid = 1
    x_steps = (x_max - x_min) / approximate_grid
    y_steps = (y_max - y_min) / approximate_grid
    x, y = numpy.mgrid[x_min:x_max:x_steps*1j, y_min:y_max:y_steps*1j]

    # make the reverse transform warping from the to_points to the from_points, because we
    # do image interpolation in this reverse fashion
    transform = _make_warp(to_points, from_points, x, y)

    if approximate_grid != 1:
        # linearly interpolate the zoomed transform grid
        new_x, new_y = numpy.mgrid[x_min:x_max+1, y_min:y_max+1]
        x_fracs, x_indices = numpy.modf((x_steps-1)*(new_x-x_min)/float(x_max-x_min))
        y_fracs, y_indices = numpy.modf((y_steps-1)*(new_y-y_min)/float(y_max-y_min))
        x_indices = x_indices.astype(int)
        y_indices = y_indices.astype(int)
        x1 = 1 - x_fracs
        y1 = 1 - y_fracs
        ix1 = (x_indices+1).clip(0, x_steps-1)
        iy1 = (y_indices+1).clip(0, y_steps-1)
        t00 = transform[0][(x_indices, y_indices)]
        t01 = transform[0][(x_indices, iy1)]
        t10 = transform[0][(ix1, y_indices)]
        t11 = transform[0][(ix1, iy1)]
        transform_x = t00*x1*y1 + t01*x1*y_fracs + t10*x_fracs*y1 + t11*x_fracs*y_fracs
        t00 = transform[1][(x_indices, y_indices)]
        t01 = transform[1][(x_indices, iy1)]
        t10 = transform[1][(ix1, y_indices)]
        t11 = transform[1][(ix1, iy1)]
        transform_y = t00*x1*y1 + t01*x1*y_fracs + t10*x_fracs*y1 + t11*x_fracs*y_fracs
        transform = [transform_x, transform_y]
    return transform

def bilinear_interpolate(x, y, bins=None):
    """Returns interpolated density values on points (x, y).
    
    Ref: http://en.wikipedia.org/wiki/Bilinear_interpolation.
    """
    if bins is None:
        bins = int(numpy.sqrt(len(x)))

    z, unused_xedge, unused_yedge = numpy.histogram2d(y, x, bins=[bins, bins],
                                        range=[(numpy.min(y), numpy.max(y)),
                                               (numpy.min(x), numpy.max(x))]
                                        )
    xfrac, xint = numpy.modf((x - numpy.min(x)) /
                             (numpy.max(x) - numpy.min(x)) * (bins - 1))
    yfrac, yint = numpy.modf((y - numpy.min(y)) /
                             (numpy.max(y) - numpy.min(y)) * (bins - 1))

    xint = xint.astype('i')
    yint = yint.astype('i')

    z1 = numpy.zeros(numpy.array(z.shape) + 1)
    z1[:-1, :-1] = z

    # values at corners of square for interpolation
    q11 = z1[yint, xint]
    q12 = z1[yint, xint + 1]
    q21 = z1[yint + 1, xint]
    q22 = z1[yint + 1, xint + 1]

    return q11 * (1 - xfrac) * (1 - yfrac) + q21 * (1 - xfrac) * (yfrac) + \
        q12 * (xfrac) * (1 - yfrac) + q22 * (xfrac) * (yfrac)

 def __new__(cls, arg1, arg2=None):
     # Assume inputs are numerical, could add an extra
     # case to parse strings as input.
     # if it is not a list, convert to a list
     if not hasattr(arg1, 'unit'):
         arg1 = arg1 * u.cycle
     if arg1.shape == ():
         arg1 = arg1.reshape((1,))
     with u.set_enabled_equivalencies(dimensionless_cycles):
         arg1 = arg1.to(u.Unit(""))
         # Since modf does not like dimensioned quantity
         if arg2 is None:
             ff,ii = numpy.modf(arg1)
         else:
             if not hasattr(arg2, 'unit'):
                 arg2 = arg2 * u.cycle
             if arg2.shape == ():
                 arg2 = arg2.reshape((1,))
             arg2 = arg2.to(u.Unit(""))
             arg1S = numpy.modf(arg1)
             arg2S = numpy.modf(arg2)
             ii = arg1S[1]+arg2S[1]
             ff = arg2S[0]
         index = numpy.where(ff < -0.5)
         ff[index] += 1.0
         ii[index] -= 1
         index = numpy.where(ff > 0.5 )
         ff[index] -= 1.0
         ii[index] += 1
         return super(Phase, cls).__new__(cls, ii.to(u.cycle), ff.to(u.cycle))

 def test_modf_array(self):
     v = np.arange(10.) * u.m / (500. * u.cm)
     q = np.modf(v)
     n = np.modf(v.to(1).value)
     assert q[0].unit == u.dimensionless_unscaled
     assert q[1].unit == u.dimensionless_unscaled
     assert all(q[0].value == n[0])
     assert all(q[1].value == n[1])

def subpx_hist(positions):
    """Historgram the decimal parts of x and y. They should be flat.
    If not, you probably do not have good sub-pixel accuracy."""
    x, y = np.array(positions)[:, 0:2].T
    fracx, fracy = np.modf(x)[0], np.modf(y)[0]
    plt.hist(fracx, bins=np.arange(0, 1.1, 0.1), color="#667788", label="x mod 1.0")
    plt.hist(fracy, bins=np.arange(0, 1.1, 0.1), color="#994433", alpha=0.5, label="y mod 1.0")
    plt.legend()
    plt.show()

def moment_xyzt(sig_xyzt, *args):  # rms=None, dc=None, ac=None):
    # ;
    # ; Calculate moments of a 4d signal of (x,y,z,t), i.e,
    # ; -RMS, i.e., a function of (x,y,t)
    # ; -DC (average in z), i.e., a function of (x,y,t)
    # ; -AC (DC subtracted out), i.e., a function of (x,y,z,t)
    # ;-------------------------------------------------------------------

    try:  # return to caller

        d = np.shape(sig_xyzt)
        if np.size(d) != 4:
            print("Error: Variable must be 4D (x,y,z,t)")
            return

        siz = np.shape(sig_xyzt)
        rms = np.zeros((siz[0], siz[1], siz[2]))
        dc = np.zeros((siz[0], siz[1], siz[2]))
        if "AC" in args:
            ac = np.zeros((siz[0], siz[1], siz[2], siz[3]))

        data = sig_xyzt
        if np.modf(np.log2(siz[3]))[0] != 0.0:
            print("WARNING: Expecting a power of 2 in Z direction")

            if np.modf(np.log2(siz[3] - 1))[0] and (siz[3] > 1):
                print(" -> Truncating last point to get power of 2")
                data = data[:, :, 0 : (siz[3] - 2), :]
                siz[3] = siz[3] - 1

        for ix in range(siz[1]):
            for iy in range(siz[2]):
                for it in range(siz[0]):
                    val = RMSvalue(sig_xyzt[it, ix, iy, :])

                    rms[it, ix, iy] = val.valrms
                    dc[it, ix, iy] = val.valav
                    if "AC" in args:
                        ac[it, ix, iy, :] = [val.acvec, val.acvec[0]]

        res = Bunch()

        if "RMS" in args:
            res.rms = rms
        if "DC" in args:
            res.dc = dc
        if "AC" in args:
            res.ac = ac

        if "RMS" not in args and "DC" not in args and "AC" not in args:
            print("Wrong argument")
        return res
    except:
        print("moment_xyz failed")
        return

def deg2dms(lat,lon, printIt=0):
    lat_frac, lat_deg = np.modf(lat)
    lon_frac, lon_deg = np.modf(lon)
    lat_min = lat_frac * 60
    lon_min = lon_frac * 60
    lon_sfrac, lon_min = np.modf(lon_min)
    lat_sfrac, lat_min = np.modf(lat_min)
    lat_sec = lat_sfrac * 60
    lon_sec = lon_sfrac * 60
    if printIt:
        print("lat: %s %s' %s\", lon: %s %s' %s\"" % (lat_deg, lat_min, lat_sec, lon_min, lon_sec))

def deg2iso(lat, lon, printIt=0):
    lat_frac, lat_deg = np.modf(lat)
    lon_frac, lon_deg = np.modf(lon)
    lat_min = lat_frac * 60.
    lon_min = lon_frac * 60.
    iso_lat_str = "%d%06.3f" % (lat_deg, abs(lat_min))
    iso_lon_str = "%d%06.3f" % (lon_deg, abs(lon_min))
    if printIt:
        print('lat: %s, lon: %s' % (iso_lat_str, iso_lon_str))
    else:
        return iso_lat_str, iso_lon_str

def degrees_to_dms(d):
    """
    Convert a floating-point degree value into a ``(degree, arcminute,
    arcsecond)`` tuple.
    """
    sign = np.copysign(1.0, d)

    (df, d) = np.modf(np.abs(d))  # (degree fraction, degree)
    (mf, m) = np.modf(df * 60.)  # (minute fraction, minute)
    s = mf * 60.

    return np.floor(sign * d), sign * np.floor(m), sign * s

def decimal_to_sexagesimal(decimal):
    """Convert decimal hours or degrees to sexagesimal.

    :param decimal: decimal number to be converted to sexagismal.
    :return: tuple of either (hours, minutes, seconds) or
             (degrees, arcminutes, arcseconds).

    """
    fractional, integral = modf(decimal)
    min_fractional, minutes = modf(fractional * 60)
    seconds = min_fractional * 60.0
    return integral.astype(int), minutes.astype(int), seconds

    def interp( self, nH, T ):

        nH = np.array( nH )
        T  = np.array( T )

        if nH.size != T.size:
            raise ValueError(' owls_ion_tables: array size mismatch !!! ')
        
        # field discovery will have nH.size == 1 and T.size == 1
        # in that case we simply return 1.0

        if nH.size == 1 and T.size == 1:
            ionfrac = 1.0
            return ionfrac


        # find inH and fnH
        #-----------------------------------------------------
        x_nH = ( nH - self.nH[0] ) / self.DELTA_nH
        x_nH_clip = np.clip( x_nH, 0.0, self.nH.size-1.001 )
        fnH,inH = np.modf( x_nH_clip )
        inH = inH.astype( np.int32 )


        # find iT and fT
        #-----------------------------------------------------
        x_T = ( T - self.T[0] ) / self.DELTA_T
        x_T_clip = np.clip( x_T, 0.0, self.T.size-1.001 )
        fT,iT = np.modf( x_T_clip )
        iT = iT.astype( np.int32 )
        

        # short names for previously calculated iz and fz
        #-----------------------------------------------------
        iz = self.iz
        fz = self.fz

                   
        # calculate interpolated value
        # use tri-linear interpolation on the log values
        #-----------------------------------------------------

        ionfrac = self.ionbal[inH,   iT,   iz  ] * (1-fnH) * (1-fT) * (1-fz) + \
                  self.ionbal[inH+1, iT,   iz  ] * (fnH)   * (1-fT) * (1-fz) + \
                  self.ionbal[inH,   iT+1, iz  ] * (1-fnH) * (fT)   * (1-fz) + \
                  self.ionbal[inH,   iT,   iz+1] * (1-fnH) * (1-fT) * (fz)   + \
                  self.ionbal[inH+1, iT,   iz+1] * (fnH)   * (1-fT) * (fz)   + \
                  self.ionbal[inH,   iT+1, iz+1] * (1-fnH) * (fT)   * (fz)   + \
                  self.ionbal[inH+1, iT+1, iz]   * (fnH)   * (fT)   * (1-fz) + \
                  self.ionbal[inH+1, iT+1, iz+1] * (fnH)   * (fT)   * (fz)

        return 10**ionfrac

	def __call__(self, ycoords, xcoords):
		if numpy.isscalar(ycoords): ycoords = (ycoords,)
		if numpy.isscalar(xcoords): xcoords = (xcoords,)
		
		yrem, yind = numpy.modf(numexpr.evaluate("(ycoords-yoffset)/yscale", local_dict=dict(ycoords=numpy.require(ycoords), yoffset=self.yoffset, yscale=self.yscale)))
		xrem, xind = numpy.modf(numexpr.evaluate("(xcoords-xoffset)/xscale", local_dict=dict(xcoords=numpy.require(xcoords), xoffset=self.xoffset, xscale=self.xscale)))
		
		yind = yind.astype(numpy.int)
		xind = xind.astype(numpy.int)
		
		res = cy.interpolator2d(yind, xind, yrem, xrem, self.data, self.ylength, self.xlength)
		
		return res

def hours_to_hms(h):
    """
    Convert an floating-point hour value into an ``(hour, minute,
    second)`` tuple.
    """

    sign = np.copysign(1.0, h)

    (hf, h) = np.modf(np.abs(h))  # (degree fraction, degree)
    (mf, m) = np.modf(hf * 60.0)  # (minute fraction, minute)
    s = mf * 60.0

    return (np.floor(sign * h), sign * np.floor(m), sign * s)