def create_grid(size, height):
    has_height = (height != None and len(height) != 0)
    row_step = (GRID_SIZE - (-GRID_SIZE)) / float(size - 1)
    col_step = (GRID_SIZE - (-GRID_SIZE)) / float(size - 1)

    # first calculate all points
    points = [(GRID_SIZE - col * col_step, has_height and height[row][col] or 0, GRID_SIZE - row * row_step) for
            row in xrange(size) for col in xrange(size)]

    horizontal_curves = []
    vertical_curves = []
    # horizontal lines
    for i in xrange(0, len(points), size):
        for j in xrange(size - 1):
            horizontal_curves.append(curve(frame=grid_frame, pos=[
                points[i + j],
                points[i + j + 1]
                ]))

    # vertical lines
    for i in xrange(0, size):
        #print "i=%d" % i
        for j in xrange(0, len(points) - size, size):
            #print "j=%d" % j
            #print "i+j=%d" % (i+j)
            #print "i+j+size=%d" % (i+j+size)
            if i + j < len(points) and i + j + size < len(points):
                vertical_curves.append(curve(frame=grid_frame, pos=[
                    points[i + j],
                    points[i + j + size]
                    ]))
    
    return (horizontal_curves, vertical_curves)

def main():
    if len(argv) < 3:
        raise Exception('>>> ERROR! Please supply values for black hole mass [>= 1.0] and spin [0.0 - 1.0] <<<')
    m = float(argv[1])
    a = float(argv[2])
    horizon = m * (1.0 + sqrt(1.0 - a * a))
    cauchy = m * (1.0 - sqrt(1.0 - a * a))
    #  set up the scene
    scene.center = (0.0, 0.0, 0.0)
    scene.width = scene.height = 1024
    scene.range = (20.0, 20.0, 20.0)
    inner = 2.0 * sqrt(cauchy**2 + a**2)
    ellipsoid(pos = scene.center, length = inner, height = inner, width = 2.0 * cauchy, color = color.blue, opacity = 0.4)  # Inner Horizon
    outer = 2.0 * sqrt(horizon**2 + a**2)
    ellipsoid(pos = scene.center, length = outer, height = outer, width = 2.0 * horizon, color = color.blue, opacity = 0.3)  # Outer Horizon
    ergo = 2.0 * sqrt(4.0 + a**2)
    ellipsoid(pos = scene.center, length = ergo, height = ergo, width = 2.0 * horizon, color = color.gray(0.7), opacity = 0.2)  # Ergosphere
    if fabs(a) > 0.0:
        ring(pos=scene.center, axis=(0, 0, 1), radius = a, color = color.white, thickness=0.01)  # Singularity
    else:
        sphere(pos=scene.center, radius = 0.05, color = color.white)  # Singularity
    ring(pos=scene.center, axis=(0, 0, 1), radius = sqrt(isco(a)**2 + a**2), color = color.magenta, thickness=0.01)  # ISCO
    curve(pos=[(0.0, 0.0, -15.0), (0.0, 0.0, 15.0)], color = color.gray(0.7))
    #cone(pos=(0,0,12), axis=(0,0,-12), radius=12.0 * tan(0.15 * pi), opacity=0.2)
    #cone(pos=(0,0,-12), axis=(0,0,12), radius=12.0 * tan(0.15 * pi), opacity=0.2)
    #sphere(pos=(0,0,0), radius=3.0, opacity=0.2)
    #sphere(pos=(0,0,0), radius=12.0, opacity=0.1)
    # animate!
    ball = sphere()  # Particle
    counter = 0
    dataLine = stdin.readline()
    while dataLine:  # build raw data arrays
        rate(60)
        if counter % 1000 == 0:
            ball.visible = False
            ball = sphere(radius = 0.2)  # Particle
            ball.trail = curve(size = 1)  #  trail
        data = loads(dataLine)
        e = float(data['v4e'])
        if e < -120.0:
            ball.color = color.green
        elif e < -90.0:
            ball.color = color.cyan
        elif e < -60.0:
            ball.color = color.yellow
        elif e < -30.0:
            ball.color = color.orange
        else:
            ball.color = color.red
        r = float(data['r'])
        th = float(data['th'])
        ph = float(data['ph'])
        ra = sqrt(r**2 + a**2)
        sth = sin(th)
        ball.pos = (ra * sth * cos(ph), ra * sth * sin(ph), r * cos(th))
        ball.trail.append(pos = ball.pos, color = ball.color)
        counter += 1
        dataLine = stdin.readline()

def drawGrid( posn=(0,0,0), sq=1, H=12, W = 1, normal='z', colour= clr.white ) :
    ht = H*sq
    wd = W*sq
    for i in range( 0, ht+1, sq ):  # FOR EACH HORIZONTAL LINE
        if normal == 'x':
            vs.curve( pos=[(posn[1]+i, posn[0], posn[2]+wd),(posn[1]+i, posn[0], posn[2])], color=colour)
            #print (posn[0] + posn[1]+i + posn[2]+wd) #for testing purposes
        else: vs.curve( pos=[(posn[0],    posn[1], posn[2]+i),
                                           (posn[0]+wd, posn[1], posn[2]+i)], color=colour)

 def draw_grid(self, width, length):
     scene.center = (int(width / 2), int(length / 2))
     scene.forward = (0, 0, -1) # Default view direction
     for x in range(-5, width * 10 + 5, 10):
         curve(pos=[(x / 10 + 0.5, -0.5), (x / 10 + 0.5, length - 0.5)],
               color=color.green)
     for y in range(-5, length * 10 + 5, 10):
         curve(pos=[(-0.5, y / 10 + 0.5), (width - 0.5, y / 10 + 0.5)],
               color=color.green)

 def __init__(self, size=100, small_interval=5, big_interval=50):
     self.frame = visual.frame(pos=(0,0,0))
     for i in range(-size, size+small_interval, small_interval):
         if i % big_interval == 0:
             c = visual.color.gray(0.65)
         else:
             c = visual.color.gray(0.25)
         visual.curve(frame=self.frame, pos=[(i,0,size),(i,0,-size)], color=c)
         visual.curve(frame=self.frame, pos=[(size,0,i),(-size,0,i)], color=c)

def axes( frame, colour, sz, posn ): # Make axes visible (of world or frame).
                                     # Use None for world.   
    directions = [vs.vector(sz,0,0), vs.vector(0,sz,0), vs.vector(0,0,sz)]
    texts = ["X","Y","Z"]
    posn = vs.vector(posn)
    for i in range (3): # EACH DIRECTION
       vs.curve( frame = frame, color = colour, pos= [ posn, posn+directions[i]])
       vs.label( frame = frame,color = colour,  text = texts[i], pos = posn+ directions[i],
                                                                    opacity = 0, box = False )

        def iterate(tree, gen, filled, spec_fill=0):
            """Itère sur l'arbre (on commence = gen=1 et filled=[0, 0])
            filled détermine la hauteur maximale 
                où on peut dessiner pour une génération donnée
            gen est la génération"""
            # On stockera dans height les différentes hauteurs des sous-arbres
            height = []
            width = gen * COLWIDTH
            altitude = 0
            current_path = None
            x = width - COLWIDTH / 2
            if len(filled) <= gen:
                filled.append(filled[-1])

            for k in xrange(len(tree[INDEX_SUBS])):
                if k == 0:
                    h, w, alt, new_path = iterate(
                        tree[INDEX_SUBS][k],
                        gen + 1,
                        filled,
                        max(filled[gen], spec_fill) - (len(tree[INDEX_SUBS]) - 1) / 2 * COLWIDTH,
                    )
                else:
                    h, w, alt, new_path = iterate(tree[INDEX_SUBS][k], gen + 1, filled)
                if w > width:
                    width = w
                if alt > altitude:
                    altitude = alt
                height.append(h)
                if new_path != None:
                    current_path = new_path
            if len(height) == 0:
                height.append(max(filled[gen], filled[gen - 1], spec_fill) + COLWIDTH / 2)
            center = (height[0] + height[-1]) / 2
            altitude = max(altitude, center + COLWIDTH / 2)

            # Les traits
            for k in xrange(len(tree[INDEX_SUBS])):
                visual.curve(pos=[(center, -x, 0), (height[k], -x - COLWIDTH, 0)], radius=LINEWIDTH)

            # Le nœud
            point = (center, -x, 0)
            if tree[INDEX_STATE] == 2:
                current_path = tree[INDEX_PATH]
                # Le nœud courant est vert
                visual.sphere(pos=point, radius=WIDTH, color=(0.0, 1.0, 0.0))
            elif tree[INDEX_STATE] == 1:
                # Si c'est lu, on met en rouge
                visual.sphere(pos=point, radius=WIDTH, color=(0.5, 0.0, 0.0))
            else:
                # Sinon, on met en rouge clair
                visual.sphere(pos=point, radius=WIDTH, color=(1.0, 0.3, 0.3))

            self.dots[tree[INDEX_PATH]] = point
            if filled[gen] < height[-1] + COLWIDTH / 2:
                filled[gen] = height[-1] + COLWIDTH / 2
            return center, width, altitude, current_path

    def plate(self):
        if (self.container):
            for obj in self.container.objects:
                obj.visible = 0
        self.container = visual.frame()
        field_length =  self.length() + 2*self.config['display']['radius']
        field_width = 3*self.config['display']['radius']*self.count + 10
        if (self.ordered):
            init = self.conns[0].obj[self.orderby]
            pborder = 0
            t = -1
            i = 0
            for conn in self.conns:
                if (init != conn.obj[self.orderby]):
                    if ((i % 2) == 1):
                        bcolor = self.config['colors']['box_light']
                    else:
                        bcolor = self.config['colors']['box']

                    if type(init) == (type(1)):
                        labeltext = self.orderby + "\n%d" % (init)
                    elif type(init) == (type('str')):
                        labeltext = self.orderby + "\n'%s'" % (init)
                    visual.box(frame=self.container, pos=(field_length/2 - self.level, -(self.config['display']['radius']+1), 3*self.config['display']['radius']*t/2+self.config['display']['radius'] + pborder/2), \
                             width = (3*self.config['display']['radius']*t - pborder), length = field_length, height = 1, color = bcolor)
                    visual.label(frame=self.container, pos = (self.config['display']['radius'] -self.level, 0,  3*self.config['display']['radius']*t/2+self.config['display']['radius'] + pborder/2),\
                        yoffset = 4*self.config['display']['radius'], xoffset = -4* self.config['display']['radius'],text = labeltext)
                    pborder = 3*self.config['display']['radius']*t+self.config['display']['radius']
                    init = conn.obj[self.orderby]
                    i += 1
                t += 1

            if type(init) == (type(1)):
                labeltext = self.orderby + "\n%d" % (init)
            elif type(init) == (type('str')):
                labeltext = self.orderby + "\n'%s'" % (init)
            if ((i % 2) == 1):
                bcolor = self.config['colors']['box_light']
            else:
                bcolor = self.config['colors']['box']
            visual.box(frame=self.container, pos=(field_length/2 - self.level, -(self.config['display']['radius']+1), 3*self.config['display']['radius']*t/2+self.config['display']['radius'] + pborder/2), \
                     width = (3*self.config['display']['radius']*t - pborder), length = field_length, height = 1, color = bcolor)
            visual.label(frame=self.container, pos = (self.config['display']['radius'] - self.level, 0, 3*self.config['display']['radius']*t/2+self.config['display']['radius'] + pborder/2),\
                yoffset = 4*self.config['display']['radius'], xoffset = -4* self.config['display']['radius'], text = labeltext)
        else:
             visual.box(frame=self.container, pos = (field_length/2 - self.level, -(self.config['display']['radius']+1),field_width/2), width = field_width, length = field_length, height = 1, color = self.config['colors']['box'])

        desctext = 'From %s to %s\n' % (time.strftime("%F %H:%M:%S", time.localtime(self.starttime)), \
            time.strftime("%F %H:%M:%S", time.localtime(self.endtime)))
        desctext += self.build_str_filter(" and ", "Filtering on ")
        visual.label(frame=self.container, pos = (field_length/2 - self.level, self.config['display']['radius']+0.5,0), yoffset = 4*self.config['display']['radius'], text = desctext)
        for i in range(self.config['display']['graduation']):
            visual.curve(frame=self.container, pos=[(field_length/self.config['display']['graduation']*i - self.level, -(self.config['display']['radius']+1)+1,0), (field_length/self.config['display']['graduation']*i - self.level,-(self.config['display']['radius']+1)+1,field_width)])
        for i in range(self.config['display']['graduation']/self.config['display']['tick']+1):
            ctime = time.strftime("%H:%M:%S", time.localtime(self.mintime + field_length/self.config['display']['graduation']*self.config['display']['tick']*i))
            visual.label(frame=self.container, pos=(field_length/self.config['display']['graduation']*self.config['display']['tick']*i - self.level, -(self.config['display']['radius']+1)+1,0), text = '%s' % (ctime), border = 5, yoffset = 1.5*self.config['display']['radius'])

    def __init__(self, npts=100000, Rotate=False, title=''):
        ptsize = 3
        mapptsize = 1
        self.scene_angle = 0
        self.localscene = visual.display(title=title + ' Local Coordinate System')
        self.globalscene = visual.display(title=title + ' Global Coordinate System')
        #self.globalscene.background = (1, 1, 1)
        #self.globalscene.up = (0, 0, 1)
        #self.globalscene.forward = (0, 1, -0.4)
        self.currentrow = 0

        self.auto_colour = False

        # Set up the global coordinate frame visualiser
        self.globalscene.select()

        self.globalscene.up =(0, 0, 1.0)

        #self.visualrobot = visual.box(length=3, height=2, width=1.5)
        w = 1
        wid = 0.2
        self.robotpts = np.array(( (0, -wid, 0, w), (0, wid, 0, w), (3*wid, 0, 0, w), (0, -wid, 0, w) ))
        #self.visualrobot = visual.points(pos=self.robotpts[:, :3], size=4, shape='round', color=(0, 1, 1))
        self.visualrobot = visual.curve(pos=self.robotpts[:, :3], color=(0, 1, 1))

        X = np.array([(-1, -1), (-1, 1), (1, 1), (1, -1), (-1, -1)])
        square = visual.curve(pos=50*X)
        self.leftmappts = visual.points(size=ptsize, shape='square', color=(1, 0, 0))
        self.rightmappts = visual.points(size=ptsize, shape='square', color=(0, 1, 0))
        self.spinmappts = visual.points(size=ptsize, shape='square', color=(0, 0, 1))
        X = np.zeros((npts, 3))
        self.mappts = visual.points(pos=X, size=mapptsize, shape='square', color=(1, 1, 1))
        self.trajpts_ind = 0
        self.trajpts = visual.points(pos=np.zeros((10000, 3)), color=(0, 0, 1), size=2)
        visual.scene.show_rendertime = True

        if Rotate:
            # Enable continuous rotation
            RV = RotateVisual(self.globalscene)
            RV.start()
        else:
            # Enable mouse panning
            MT = dispxyz.EnableMouseThread(self.globalscene)
            MT.start()

        # Set up the local coordinate frame visualiser
        if showlocal:
            self.localscene.select()
            self.leftpts = visual.points(color=(1, 0, 0), size=ptsize, shape='square')
            self.rightpts = visual.points(color=(0, 1, 0), size=ptsize, shape='square')
            self.spinpts = visual.points(color=(0, 0, 1), size=ptsize, shape='square')
            visual.scene.show_rendertime = True
        
        self.colourmin = -4
        self.colourmax = 4

def set_scene(R):        # draw scene, ball, trails, spin, info box
    scene = vp.display(background=(.2,.5,1), forward=(-1,-.1,-.1),
                       center=(.5*R,1,0), ambient=.4, fullscreen=1)
    floor = vp.box(pos=(R/2,0,0), length=1.1*R, height=.1, width=8, 
                   color=vp.color.orange, opacity=0.7)  # transparent 
    zone = vp.curve(pos=[(R,0,1),(R,1,1),(R,1,-1),(R,0,-1)], radius=.02)
    ball = vp.sphere(pos=(0,0,0), radius=.2, material=vp.materials.rough)
    trail = vp.curve(pos=(0,0,0), radius=0.04)
    ideal = vp.curve(pos=(0,0,0), radius=0.04, color=vp.color.green)
    spin = vp.arrow(axis=omega,pos=(0,0,0),length=1)   # omega dir
    info = vp.label(pos=(1.1*R,2,-2),text='Any key=repeat')
    return scene, ball, trail, ideal, spin

def calculate(x, y, z, vx, vy, vz, dt, m, g, B2, S0, omega):
    """ Calculate the trajectory of a baseball including air resistance and spin by
repeatedly calling the do_time_step function.  Also draw the trajectory using visual python. """
    t = 0.0
    # Establish lists with initial position and velocity components and time.
    x_list = [x]
    y_list = [y]
    z_list = [z]
    vx_list = [vx]
    vy_list = [vy]
    vz_list = [vz]
    t_list = [t]

    # Set up visual elements.
    mound = visual.box(pos=(0,0,0), length=0.1, width=0.5, height=0.03, color=visual.color.white)
    plate = visual.box(pos=(18,0,0), length=0.5, width=0.5, height=0.03, color=visual.color.white)
    ball = visual.sphere(pos=(x,y,z), radius=0.05, color=visual.color.white)
    ball.trail = visual.curve(color=ball.color)

    while y >= 0.0:
        visual.rate(100) # Limit to no more than 100 iterations per second.
        t, x, y, z, vx, vy, vz = do_time_step(t, dt, x, y, z, vx, vy, vz, m, B2, g, S0, omega)
        x_list.append(x)
        y_list.append(y)
        z_list.append(z)
        vx_list.append(vx)
        vy_list.append(vy)
        vz_list.append(vz)
        t_list.append(t)
        ball.pos = (x,y,z)
        ball.trail.append(pos=ball.pos)

    return t_list, x_list, y_list, z_list, vx_list, vy_list, vz_list

def golf_ball_calc(x,y,z,vx,vy,vz,dt,m,g,B2,S0,w,w_vector):
    t = 0.0
    x_list = [x]
    y_list = [y]
    z_list = [z]
    vx_list = [vx]
    vy_list = [vy]
    vz_list = [vz]
    t_list = [t]
    v_vector = visual.vector(vx,vy,vz)

    tee = visual.box(pos=(0,0,0), length=0.05, width=0.05, height=0.5,color=visual.color.white)
    ball = visual.sphere(pos=(x,y,z), radius = 0.25, color = visual.color.white)
    ball.trail = visual.curve(color = visual.color.red)

    while y > 0.0:
        visual.rate(100)
        t,x,y,z,vx,vy,vz = golfball_step(t,x,y,z,vx,vy,vz,m,g,B2,S0,w,dt,w_vector,v_vector)
        x_list.append(x)
        y_list.append(y)
        z_list.append(z)
        vx_list.append(vx)
        vy_list.append(vy)
        vz_list.append(vz)
        t_list.append(t)
        v_vector = visual.vector(vx,vy,vz)
        ball.pos = (x,y,z)
        ball.trail.append(pos=ball.pos)

    return t_list,x_list,y_list,z_list,vx_list,vy_list,vz_list

 def frame(self, sign):
     self.glass = curve(
         pos=[
             (sign * self.time, -self.time * 1.1, 0),
             (sign * self.time * 0.05, -self.time * 0.1, 0),
             (sign * self.time * 0.05, self.time * 0.1, 0),
             (sign * self.time, self.time * 1.1, 0),
             (sign * self.time, self.time * 1.15, 0),
         ],
         radius=self.time * 0.025,
     )
     self.base = cylinder(
         pos=(0, sign * self.time * 1.15, 0),
         axis=(0, 1, 0),
         length=self.time * 0.1,
         radius=self.time * 1.2,
         color=(0.66, 0.46, 0.13),
     )
     self.pole = cylinder(
         pos=(sign * self.time, -self.time * 1.1, 0),
         axis=(0, 1, 0),
         length=self.time * 2.3,
         radius=self.time * 0.06,
         color=(0.66, 0.46, 0.13),
     )

    def display(self,data):
        #dominant = reduce(lambda y,x: x if x[1] > y[1] else y, zip(range(len(data)), data), [0,0])
        #if (dominant[1] > 10 and dominant[0] != 0.0):
        #    freq = freqs[dominant[0]]
        #    print freq,
        #    note = ftomidi(freq)
        #    print note
        #    PlayNote(note)
        #    dominant = dominant[0]
        #else:
        #    dominant = False

        if (len(self.curves) > self.maxlen):
            todel = self.curves.pop()
            todel.visible = False
            del todel

        for oldcurve in self.curves:
            for point in range(len(oldcurve.pos)):
                oldcurve.pos[point][2] -= self.step

        y = []
        for point in data:
            y.append(point)
            y.append(point)

        curve = visual.curve(color=visual.color.white, display=self.g, radius=0)

        for point in zip(range(len(data)),data):
            r = g = b = 1
            curve.append(pos=(point[0],point[1] * 10,0), color=(r,g,b))
        self.curves.insert(0,curve)

 def run(self):
     '''
     Run the simulation with racers that had been previously added to the world
     by add_racer method.
     '''
     # create the scene with the plane at the top
     visual.scene.center = visual.vector(0,-25,0)
     visual.box(pos=(0,0,0), size=(12,0.2,12), color=visual.color.green)
     # create the visual objects that represent the racers (balls)
     balls = [ visual.sphere(pos=(index,0,0), radius=0.5) for index in xrange(len(self.racers))]
     for ball, racer in zip(balls, self.racers):
         color = visual.color.blue
         try:
             # try to set the color given by a string
             color = getattr(visual.color, racer.color)
         except AttributeError:
             pass
         ball.trail  = visual.curve(color=color)
     while not reduce(lambda x, y: x and y, [racer.result for racer in self.racers]):
         # slow down the looping - allow only self.time_calibration
         # number of loop entries for a second
         visual.rate(self.time_calibration)
         # move the racers
         for racer in self.racers:
             self.__move_racer(racer) 
             
         for ball, racer in zip(balls, self.racers):
             ball.pos.y = -racer.positions[-1][1]
             ball.pos.x = racer.positions[-1][0]
             ball.trail.append(pos=ball.pos)
             
         self.current_time += self.interval
         self.timeline.append(self.current_time)

def extractXCurve(point):
    """after forwardMap has been called, call this 
    on the x parameter to get a visual curve of the elastica"""
    vcurve = vis.curve()
    for part in point.x:
        vcurve.append([part[0], part[1]])
    return vcurve

 def __init__(self, **kwargs):
     super(Board, self).__init__(**kwargs)
     
     visual.box(pos = (9, 9, -1),
                size = (20, 20, 1),
                color = (1, 0.75, 0),
                material = visual.materials.wood,
                frame = self)
     # lines
     params = dict(color = (.5, .5, .5), frame = self)
     for i in range(19):
         visual.curve(pos = [(0, i, 0), (18, i, 0)], radius = .05, **params)
     for i in range(19):
         visual.curve(pos=[(i, 0, 0), (i, 18, 0)], radius = .05, **params)
     
     self.stones, self.star = [], []
     for x in itertools.product((4, 10, 16), (4, 10, 16)): self.star.append(x[:])

 def drawRay(self):
     #print "Drawing Ray"
     #print self.points
     self.ray = curve(pos=self.points,
                      color=(0,0,0),
                      #radius=0.005
                      radius=0.002
                      )

def set_scene(R, r):        # create bodies, velocity arrows
    vp.display(title='Three-body motion', background=(1,1,1))
    body, vel = [], []      # bodies, vel arrows
    c = [(1,0,0), (0,1,0), (0,0,1), (0,0,0)]    # RGB colors
    for i in range(3):
        body.append(vp.sphere(pos=r[i],radius=R,color=c[i],make_trail=1))
        vel.append(vp.arrow(pos=body[i].pos,shaftwidth=R/2,color=c[i]))
    line, com = vp.curve(color=c[3]), vp.sphere(pos=(0,0), radius=R/4.)
    return body, vel, line

def update_grid():
    delete_current_grid()

    if not showSpan:
        return 

    global span_frame

    span_frame = visual.frame() 

    vecs = map(lambda i: visual.vector(vectors[i].axis)*span_fact,selected)
    for v in vecs:
        rest = vecs[0:len(vecs)] # deep copy  
        rest.remove(v)
        all_sums = all_sums_comb(rest)
        for aSum in all_sums:
            visual.curve(frame = span_frame, pos = [aSum,aSum+v], color = span_color)
            if not is_span_positive:
                visual.curve(frame = span_frame, pos = [aSum,aSum-v], color = span_color)