def mjd2gmst(mjd):
    """
    Returns GMST given UT1 in the form of an MJD

    The GMST is returned in radians from 0 to 2pi. This
    was converted from Fortran source code of SOFA
    """

    # Set some constants
    DS2R   = 7.272205216643039903848712e-5
    DJ0    = 2451545e0
    DJ1    = 2400000.5e0
    DAYSEC = 86400e0
    CENDAY = 36525e0
    A      = 24110.54841e0 - DAYSEC/2.
    B      = 8640184.812866e0
    C      = 0.093104e0
    D      = -6.2e-6

    if DJ1 < mjd:
        d1 = DJ1
        d2 = mjd
    else:
        d1 = mjd
        d2 = DJ1

    t = (mjd + (DJ1-DJ0 ))/CENDAY

    f = DAYSEC*(0.5 + m.fmod(mjd,1.))

    return m.fmod(DS2R*((A+(B+(C+D*t)*t)*t)+f),2.*m.pi)

def onMouseMove(x, y):
    global mouseDrag, mouseDragX, mouseDragY, mouseDragMove
    global eyeX, eyeY, eyeZ, upX, upY, upZ, cam_theta, cam_phi, cam_r
    if not mouseDrag:
        return
    mouseDragMove = True
    # Mouse point to angle conversion
    # cam_theta = (360.0/winHeight)*y*3.0#3.0 rotations possible
    # cam_phi = (360.0/winWidth)*x*3.0
    cam_phi += 360.0 * (mouseDragX - x) / winWidth * 2.0
    cam_theta += 360.0 * (mouseDragY - y) / winHeight * 2.0
    mouseDragX = x
    mouseDragY = y
    # Restrict the angles within 0~360 deg (optional)
    if cam_theta > 360:
        cam_theta = fmod(cam_theta, 360.0)
    if cam_phi > 360:
        cam_phi = fmod(cam_phi, 360.0)
    newValues = lookInSphere(cam_r, cam_phi, cam_theta)
    eyeX = newValues[0]
    eyeY = newValues[1]
    eyeZ = newValues[2]
    upX = newValues[3]
    upY = newValues[4]
    upZ = newValues[5]
    glutPostRedisplay()

def arange(start, stop, n):
    start = fmod(start, 2*pi)
    stop = fmod(start, 2*pi)
    if(fabs(start - stop) > pi):
        if start > stop: start -= 2*pi
        else: stop -= 2*pi
    return frange(start, stop, n)

    def round_to_sum(l, r):
        """
        Round a list of numbers while maintaining the sum.

        http://stackoverflow.com/questions/15769948/
        round-a-python-list-of-numbers-and-maintain-the-sum

        :param l: array
        :type l: list(float)

        :param r: decimal place
        :type r: int

        :returns: A list of rounded numbers whose sum is equal to the
            sum of the list of input numbers.
        :rtype: list
        """
        q = 10 ** (-r)
        d = (round(sum(l), r) - sum([round(x, r) for x in l])) * (10 ** r)
        d = int(d)
        if d == 0:
            return [round(x, r) for x in l]
        elif d in [-1, 1]:
            c, _ = max(enumerate(l), key=lambda x: math.copysign(
                1, d) * math.fmod(x[1] - 0.5 * q, q))
            return [round(x, r) + q * math.copysign(1, d) if i == c else round(
                x, r) for (i, x) in enumerate(l)]
        else:
            c = [i for i, _ in heapq.nlargest(abs(d), enumerate(
                l), key=lambda x: math.copysign(1, d) * math.fmod(
                    x[1] - 0.5 * q, q))]
            return [round(x, r) + q * math.copysign(
                1, d) if i in c else round(x, r) for (i, x) in enumerate(l)]

def w_h_from_str(string, source_width, source_height):
    if ':' in string:  # used `w:h` syntax
        w, h = string.split(':')
        w = int(w)
        h = int(h)
        if w < -1 or h < -1:
            raise ValueError('unknown negative component')
        # keep aspect ratio if either component is -1
        if w == -1 and h == -1:  # ffmpeg allows this so we should too
            return source_width, source_height
        else:
            if w == -1:
                w = int(h * source_width / source_height) & ~1  # round to even
            elif h == -1:
                h = int(w * source_height / source_width) & ~1
    elif float(string) != 1.0:  # using a singlular int or float
        # round to nearest even number
        even_pixel = lambda x: \
            int(math.floor(x * float(string) / 2) * 2)
        w = even_pixel(source_width)
        h = even_pixel(source_height)
    else:  # use source w,h by default
        w = source_width
        h = source_height
    # w and h must be divisible by 2 for yuv420p outputs
    # don't auto round when using the `w:h` syntax (with no -1 components)
    # because the user may not expect the changes
    if math.fmod(w, 2) != 0 or math.fmod(h, 2) != 0:
        raise ValueError('components not divisible by two')
    return w, h

def analyze(arrStack, mean, std):
    print 'tracking.analyze ...'
    shape = numpy.shape(std)
    hsig  = TH2F('hsig', '', shape[0], 0, shape[0]+1, shape[1], 0, shape[1]+1)
    hsigw = TH1F('hsigw', '', 100, 0, 50)

    # 2D to 1D
    std = numpy.ravel(std)
    mean = numpy.ravel(mean)
    h3 = TH3F('h3', '', shape[0], 0, shape[0], shape[1], 0, shape[1], len(arrStack), 0, len(arrStack))

    for i, arr in enumerate(arrStack):
        if i%100 == 0:
            print 'processing ... {0}/{1}'.format(i, len(arrStack))
        hsig_tmp  = TH2F('hsig_{0:04d}'.format(i), '', shape[0], 0, shape[0]+1, shape[1], 0, shape[1]+1)
        # print 'hsig_{0:04d}'.format(i)

        arr = numpy.ravel(arr)
        pixels, pixel_indices = findPixels(arr, mean, std, i)

        for index, x in numpy.ndenumerate(pixels):
            hsig.Fill(math.floor(pixel_indices[index[0]]/shape[0]), math.fmod(pixel_indices[index[0]], shape[0]), x)
            hsig_tmp.Fill(math.floor(pixel_indices[index[0]]/shape[0]), math.fmod(pixel_indices[index[0]], shape[0]), x)
        hsig_tmp.Write()

        hits = findHits(pixels, pixel_indices, shape[0], shape[1], i, h3)

    h3.Write()

def _rotate_hcline(obj, objdict, cx, cy, ra):
    _layer = obj.getParent()
    if _layer is None:
        raise RuntimeError, "HCLine parent is None"
    _lp = obj.getLocation()
    if _lp.getParent() is not _layer:
        raise RuntimeError, "HCLine/Point parent object conflict!"
    _x, _y = _calc_coords(_lp, cx, cy, ra)
    _pts = _layer.find('point', _x, _y)
    if len(_pts) == 0:
        _np = Point(_x, _y)
        _layer.addObject(_np)
    else:
        _np = _most_used(_pts)
    _da = ra/_dtr
    if abs(fmod(_da, 180.0)) < 1e-10:
        obj.setLocation(_np)
        if _adjust_dimension(_lp, _np):
            _layer.delObject(_lp)
    elif abs(fmod(_da, 90.0)) < 1e-10:
        _layer.addObject(VCLine(_np))
        if _adjust_dimension(_lp, _np):
            _layer.delObject(obj)
    else:
        _layer.addObject(ACLine(_np, _da))
        _layer.delObject(obj)
    _pid = id(_lp)
    if objdict.get(_pid) is not False:
        objdict[_pid] = False

    def test_deg_km_accuracy(self):
        c = quakelib.Conversion(quakelib.LatLonDepth(0,0))

        # Check that 360 * length of 1 longitude degree is equal to the circumference of the equator
        # Confirm accuracy is within 1 meter
        one_deg_len = c.convert2xyz(quakelib.LatLonDepth(0,1)).mag()
        self.assertAlmostEqual(one_deg_len*360.0/1000, 40075.016, 2)

        # Check that 4 * length of 90 degree vertical arc is equal to the polar circumference
        # Confirm accuracy is within 1 meter
        ninety_deg_len = c.convert2xyz(quakelib.LatLonDepth(90,0)).mag()
        self.assertAlmostEqual(ninety_deg_len*4.0/1000, 40007.860, 2)

        # Check that inverse of conversion results in the same value
        for base_lat in range(-90,91,5):
            for base_lon in range(-180, 180, 5):
                base_pt = quakelib.LatLonDepth(base_lat, base_lon)
                conv = quakelib.Conversion(base_pt)
                test_lat = math.fmod(base_lat+random.uniform(-45,45), 90)
                test_lon = math.fmod(base_lon+random.uniform(-45,45), 180)
                test_pt = quakelib.LatLonDepth(test_lat, test_lon)
                new_xyz = conv.convert2xyz(test_pt)
                rev_pt = conv.convert2LatLon(new_xyz)
                # Ensure accuracy to within 1e-7 degrees (~1 cm)
                self.assertAlmostEqual(test_lat, rev_pt.lat(), 7)
                self.assertAlmostEqual(test_lon, rev_pt.lon(), 7)

def eq2gal(ra, dec):

    rmat = numpy.array(
        [
            [-0.054875539726, -0.873437108010, -0.483834985808],
            [+0.494109453312, -0.444829589425, +0.746982251810],
            [-0.867666135858, -0.198076386122, +0.455983795705],
        ],
        dtype="d",
    )
    cosb = math.cos(dec)
    v1 = numpy.array([math.cos(ra) * cosb, math.sin(ra) * cosb, math.sin(dec)])
    v2 = numpy.dot(rmat, v1)
    x = v2[0]
    y = v2[1]
    z = v2[2]
    r = math.sqrt(x * x + y * y)
    if r == 0.0:
        l = 0.0
    else:
        l = math.atan2(y, x)
    if z == 0.0:
        b = 0.0
    else:
        b = math.atan2(z, r)
    ll = math.fmod(l, 2.0 * math.pi)
    if ll < 0.0:
        ll = ll + 2.0 * math.pi

    bb = math.fmod(b, 2 * math.pi)
    if abs(bb) >= math.pi:
        print "Ugh!"

    return ll, bb

def slice_linear(zmin, zmax, n, scale=None):
    """The size of result is (n-1)."""
    COMBS = [
        None,  # not used
        (0, ),
        (0, 5, ),
        (0, 3, 7, ),
        (0, 3, 5, 7, ),
        (0, 2, 4, 6, 8, ),
        (0, 2, 3, 5, 7, 8, ),
        (0, 1, 3, 4, 6, 7, 9, ),
        (0, 1, 3, 4, 5, 6, 7, 9, ),
        (0, 1, 2, 3, 4, 6, 7, 8, 9, ),
        (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ),
    ]
    if scale is None:
        scale = 10 * ceil_pow10((zmax - zmin) / n)
    step = scale // 10
    def generate(comb, bases):
        if comb is None:
            return []
        result = [base + _ * step for base in bases for _ in comb]
        return [_ for _ in result if zmin < _ < zmax]
    zmin_floored = math.floor(zmin - math.fmod(zmin, scale))  # XXX: these might be numerically unstable
    zmax_ceiled = math.ceil(zmax - math.fmod(zmax, scale))
    multiples = list(range(zmin_floored, zmax_ceiled + scale, scale))
    multiples = [_ for _ in multiples if zmin < _ < zmax]
    if len(multiples) == 0:
        k = (zmax_ceiled - zmin_floored) // step
        if k == 0:
            return slice_linear(zmin, zmax, n, scale=(scale // 10))
        result = generate(COMBS[k], [zmin_floored])
        # TODO: compensate len(result) to be == (n - 1)
        return result
    k = n // len(multiples)
    if k <= 0:
        return []
    if k > 10:
        return slice_linear(zmin, zmax, n, scale=(scale // 10))
    result_body = generate(COMBS[k], multiples[:-1])
    result_head = []
    result_tail = []
    # len(result_head + result_body + result_tail) should be == (n - 1)
    len_needed = (n - 1) - len(result_body)
    assert len_needed >= 0
    head_length = len_needed // 2
    tail_length = (len_needed + 1) // 2
    for comb in COMBS[:(k + 1)]:
        candidate = generate(comb, [multiples[0] - scale])
        if len(candidate) == head_length:
            result_head = candidate
        candidate = generate(comb, [multiples[-1]])
        if len(candidate) == tail_length:
            result_tail = candidate
    result = result_head + result_body + result_tail
    if len(result) != (n - 1):
        # n is too large  (e.g. zmin = 10, zmax = 1000, n = 50)
        # Once again, do the linear slicing instead
        return slice_linear(zmin, zmax, n, scale=(scale // 10))
    return result

    def pattern(self):
        self.trap.send_decay(ALL, 90)
        self.trap.start_pattern(ALL)

        index = 0
        hue_offset = 0.0
        stop = False

        color_rings = [ random_color(), random_color(), random_color() , random_color() ]
        while not stop:
            for bottle, angle in self.geo.enumerate_all_bottles(index % 2 == 0):
                self.trap.set_color(bottle, color_rings[self.geo.get_ring_from_bottle(bottle)])
                sleep(.01)
                if self.stop_thread:
                    stop = True
                    break

            index += 1
            hue_offset = math.fmod(hue_offset + .02, 1.0)
            shift = math.sin(index / self.color_shift_between_rings) / 2.0 + .50
            new_offset = math.fmod(shift, 1.0)
            color_rings.pop()
            color_rings.insert(0, hue_to_color(new_offset))

        self.trap.stop_pattern(ALL)
        if self.transition:
            sleep(.02)
            transition_sweep_out(self.trap)

def schwefel_func(x, Os=None, Mr=None, sr=None, adjust=1.0):
    nx = len(x)
    # z = sr_func(x, Os, Mr, 1000.0 / 100.0)
    adjust = 0.2
    z = sr_func(x, Os, Mr, 100000.0 / 100.0 * adjust)
    f = 0.0
    dim = float(len(x))
    for i in range(nx):
        z[i] += 4.209687462275036e+002
        if z[i] > 500:
            f -= (500.0 - fmod(z[i], 500)) * sin(pow(500.0 - fmod(z[i], 500),
                                                     0.5))
            tmp = (z[i] - 500.0) / 100.0
            f += tmp * tmp / nx
        elif z[i] < -500:
            f -= (-500.0 + fmod(fabs(z[i]), 500)) * sin(
                pow(500.0 - fmod(fabs(z[i]), 500), 0.5))
            tmp = (z[i] + 500.0) / 100
            f += tmp * tmp / nx
        else:
            f -= z[i] * sin(pow(fabs(z[i]), 0.5))
        f += 4.189828872724338e+002 * nx
    # return max(f / 1000.0 - 0.838, 0.0)
    # return max(f / 1000.0, 0.0)
    return max((f - (dim * (dim - 1)) * 4.189e+002) / 1000.0, 0.0)

def removeVars( varList ) :

    # ======================
    #  remove unwanted vars
    # ======================

    done           =  False

    varAddCounter  =  0

    while not done :

        if varList :

            varAddCounter += 1

            curMsg  =  '  curVars   :  '

            for ( i , curVar ) in enumerate( varList ) :

                if i != 0  :

                    curMsg +=  ' ,  '

                    if int( fmod( i , 2 ) ) == 0 :

                        curMsg +=  '  \n               '

                curMsg += curVar.ljust( 14 )

            if int( fmod( varAddCounter , 5 ) ) == 0 :

                print varMsg

            print curMsg + '\n'

        varInput  =  raw_input( inputRequest )

        print ''

        # -------------------------------------
        #  check if varList has been completed
        # -------------------------------------

        if not varInput :

            if varList :

                done = True

            else :

                print '\n You need to insert at least one variable to hunt. \n'
                continue

        varInput =  varInput.replace( ',' , '' )

        varList.extend(  varInput.split()  )

    return varList

def draw_palette(H,height,width): 

    H = fmod(H,360) if H >= 360 else 360+fmod(H,360) if H < 0 else H
    name = "_".join([str(H),str(height),str(width)])
    if name not in lookup: 
        p2=[]
        V_step = 100.0/height
        S_step = 100.0/width
        for row in range(height):
            V = 100 - V_step * row
            if V<=0 : V = 1
            l=[]
            for col in range(width):
                S = 100 - S_step * col
                if S<=0 : S = 1
                hex = rgb2hex(hsv2rgb([H,S,V]))
                l.append(hex)
            p2.append(l)
        lookup[name] = p2        # obmit coping as 
                                 # the list ref will not be changed
                                 # when we use [] to creat a new one
    else:
        p2 = lookup[name]
    
    vcmd("call s:clear_match('pal')")
    for row  in range(height):
        for col in range(width):
            hex = p2[row][col]
            hi_grp  = "".join(["cv_pal_",str(row),"_",str(col)])
            pos_ptn = "".join(["\\%",str(row+H_OFF+1),"l\\%"
                                    ,str(col+W_OFF+1),"c"])
            vcmd("call s:hi_color('{}','{}','{}',' ')".format(hi_grp,hex,hex))
            vcmd("sil! let s:pal_dict['{0}'] = matchadd('{0}','{1}')"
                    .format( hi_grp,pos_ptn))
    vcmd("".join(["let s:pal_clr_list=",str(p2)]))

  def differentialWheelsSetSpeed(self, obj, left, right):
    if self.leftWheelMaxSpeed < left :
      left = self.leftWheelMaxSpeed
    elif left < -self.leftWheelMaxSpeed :
      left = -self.leftWheelMaxSpeed

    if self.leftWheelSpeedUnit > self.Accueacy :
      radiuse = math.fmod(math.fabs(left), self.leftWheelSpeedUnit)
      if  left > 0:
        left -= radiuse
      else:
        left += radiuse
      
    if self.rightWheelMaxSpeed < left :
      left = self.rightWheelMaxSpeed
    elif left < -self.rightWheelMaxSpeed :
      left = -self.rightWheelMaxSpeed

    if self.rightWheelSpeedUnit > self.Accueacy :
      radiuse = math.fmod(math.fabs(left), self.rightWheelSpeedUnit)
      if  right > 0:
        right -= radiuse
      else:
        right += radiuse
      
   
    if  self.leftWheelName:
      obj.setAngularVelocityToParts(self.leftWheelName, left, self.leftMaxForce)
      self.currentLeftWheelSpeed = left

    if  self.rightWheelName:
      obj.setAngularVelocityToParts(self.rightWheelName, right, self.rightMaxForce)
      self.currentRightWheelSpeed = right

    return 

def wrap_angle(angle):
    if angle <= math.pi and angle >= -math.pi:
        return angle
    elif angle < 0.0:
        return math.fmod(angle-math.pi,2.0*math.pi)+math.pi
    else:
        return math.fmod(angle+math.pi,2.0*math.pi)-math.pi;

    def from_jd(cls, jd):
        """
            Convert a Julian day number to a year/month/day tuple
            of this calendar (matching jQuery calendars algorithm)

            @param jd: the Julian day number
        """

        jd = math.floor(jd) + 0.5;

        depoch = jd - cls.to_jd(475, 1, 1)

        cycle = math.floor(depoch / 1029983)
        cyear = math.fmod(depoch, 1029983)

        if cyear != 1029982:
            aux1 = math.floor(cyear / 366)
            aux2 = math.fmod(cyear, 366)
            ycycle = math.floor(((2134 * aux1) + (2816 * aux2) + 2815) / 1028522) + aux1 + 1
        else:
            ycycle = 2820

        year = ycycle + (2820 * cycle) + 474
        if year <= 0:
            year -= 1

        yday = jd - cls.to_jd(year, 1, 1) + 1
        if yday <= 186:
            month = math.ceil(yday / 31)
        else:
            month = math.ceil((yday - 6) / 30)

        day = jd - cls.to_jd(year, month, 1) + 1

        return (int(year), int(month), int(day))

def pCJDNToMilankovic(iCJDN, bDoYear=False, bDoMonth=False, bDoDay=False):
    ''' Returns a Revised Julian date is ISO 8601 YYYY-MM-DD format from a given
        Chronological Julian Day Number (CJDN).
        If only part of the date is required (e.g. just the day), this is returned as an integer.
    '''

    if not isinstance(iCJDN, int):
        return False

    #Perform the calculation
    iK3 = (9 * (iCJDN - 1721120)) + 2
    iX3 = floor(iK3 / 328718)
    iK2 = (100 * floor(int(fmod(iK3, 328718)) / 9)) + 99
    iX2 = floor(iK2 / 36525)
    iK1 = (5 * floor(int(fmod(iK2, 36525)) / 100)) + 2
    iX1 = floor(iK1 / 153)
    iC0 = floor((iX1 + 2) / 12)
    iYear = (100 * iX3) + iX2 + iC0
    iMonth = (iX1 - (12 * iC0)) + 3
    iDay = floor(int(fmod(iK1, 153)) / 5) + 1

    #Return integer part, if only one part of the date is required.
    #Start at year, and rudely ignore any subsequent integers requested.
    if bDoYear:
        return iYear
    if bDoMonth:
        return iMonth
    if bDoDay:
        return iDay

    #Return the ISO 8601 date string
    sISO8601Date = ('0000' + str(iYear))[-4:] + '-' + ('00' + str(iMonth))[-2:] + '-' + ('00' + str(iDay))[-2:]
    return sISO8601Date

    def complementaries(self, base_color):
        # return some other colors that "go" with this one
        hsv = RGBtoHSV(base_color)

        (h,s,v,a) = hsv
        
        # take 2 colors which are almost triads
        h = hsv[0]
        delta = random.gauss(0.0, 0.8)
        h2 = math.fmod(h + 2.5 + delta, 6.0)
        h3 = math.fmod(h + 3.5 - delta, 6.0)

        # take darker and lighter versions
        v = hsv[2]
        vlight = self.clamp(v * 1.5, 0.0, 1.0)
        vdark = v * 0.5

        colors = [
            [h, s, vdark, a],
            [h, s, v, a],
            [h, s, vlight, a],
            [h2, s, vlight, a],
            [h2, s, v, a],
            [h2, s, vdark, a],
            [h3, s, vdark, a],
            [h3, s, v, a],
            [h3, s, vlight, a]]

        colors = [ HSVtoRGB(x) for x in colors]
        return colors

def pCJDNToJulian(iCJDN, bDoYear=False, bDoMonth=False, bDoDay=False):
    ''' Returns a Julian date is ISO 8601 YYYY-MM-DD format from a given
        Chronological Julian Day Number (CJDN).
        If only part of the date is required (e.g. just the day), this is returned as an integer.
    '''

    if not isinstance(iCJDN, int):
        return False

    #Perform the calculation
    iY2 = iCJDN - 1721118
    iK2 = (4 * iY2) + 3
    iK1 = (5 * floor(int(fmod(iK2, 1461)) / 4)) + 2
    iX1 = floor(iK1 / 153)
    iC0 = floor((iX1 + 2) / 12)
    iYear = floor(iK2 / 1461) + iC0
    iMonth = (iX1 - (12 * iC0)) + 3
    iDay = floor(int(fmod(iK1, 153)) / 5) + 1

    #Return integer part, if only one part of the date is required.
    #Start at year, and rudely ignore any subsequent integers requested.
    if bDoYear:
        return iYear
    if bDoMonth:
        return iMonth
    if bDoDay:
        return iDay

    #Return the ISO 8601 date string
    sISO8601Date = ('0000' + str(iYear))[-4:] + '-' + ('00' + str(iMonth))[-2:] + '-' + ('00' + str(iDay))[-2:]
    return sISO8601Date