diff --git a/Final/python files/couzin_vicsek_simulation.py b/Final/python files/couzin_vicsek_simulation.py
new file mode 100644
index 0000000000000000000000000000000000000000..e6daa17b0c32bb681dae5489058d0aba3d218e85
--- /dev/null
+++ b/Final/python files/couzin_vicsek_simulation.py	
@@ -0,0 +1,250 @@
+from graphics import * # get here: http://mcsp.wartburg.edu/zelle/python/graphics.py
+import time
+import random
+
+# imports from our files
+from help_functions import *
+from classes import Agent
+from check_boundarys import *
+
+
+# --------------------------------------------------------------
+# global variables
+
+    # example zones for couzin, vicsek
+    # radii of zones:           zor,zoo,zoa
+    # swarm:                    10, 20, 200
+    # torus:                    5,  60, 200
+    # dynamic parallel group:   5,  100,200
+    # highly parallel group:    5,  180,200
+
+# window dimensions
+winWidth = 500
+winHeight = 500
+
+window = GraphWin("Window", winWidth, winHeight)
+
+# max time for one simulation
+maxTime = 4000
+
+# whether the boundary is "free" or not
+# <->
+# can they swim through the side of the window y/n?
+free = False
+
+# max velocity of the agents
+maxV = 8
+
+# constant speed of agents
+speed = 4
+
+# number of agents
+agentNum = 50
+
+# angle of blindness
+blind = 80
+
+# maximum turning angle
+maxTurn = 50
+radTurn = math.radians(maxTurn)
+negRadTurn = math.radians(360-maxTurn)
+
+# radii of zones
+zor_r = 10
+zoo_r = 80
+zoa_r = 110
+
+# True = couzin, False = Vicsek
+couzin = True
+
+if not couzin:
+    # vicsek simulation (= couzin with specific parameters)
+    maxTurn = 180
+    radTurn = math.radians(maxTurn)
+    negRadTurn = -math.radians(maxTurn)
+    zor_r = 0
+    zoa_r = 0
+    blind = 0
+    #free = True
+
+
+# ----------------------------------------------------------------
+# couzin
+
+# returns three lists, one for each zone,
+# contaning all other agents in the zone.
+# ignores all agents in the angle behind the current agent defined by blind.
+# the zones zor/zoo/zoa are defined by their radii zor_r/zoo_r/zoa_r
+def neigbour_in_zones(agent, agents):
+    zor = []
+    zoo = []
+    zoa = []
+
+    for neighbour in agents:
+        # find distance in x and y direction,
+        # and use them to calculate the angle between
+        # the movement-vector of agent and
+        # the directiction-vector between agent and neighbour
+        
+        disVecX = neighbour.point.getX() - agent.point.getX()
+        disVecY = neighbour.point.getY() - agent.point.getY()
+        alpha = calc_angle(agent.x_velocity, agent.y_velocity, disVecX, disVecY)
+
+        if (agent == neighbour):
+            True
+        elif alpha < 180 - blind and alpha > 180 + blind:
+            True
+
+        # if neighbour can be seen:
+        # get distance dis between agent and neighbour,
+        # and check if neighbour is in a zone of agent
+        else:
+            dis = absvec(neighbour.point.getX() - agent.point.getX(), neighbour.point.getY() - agent.point.getY())
+            if dis <= zor_r:
+                zor.append(neighbour)
+            elif dis <= zoo_r:
+                zoo.append(neighbour)
+            elif dis <= zoa_r:
+                zoa.append(neighbour)
+
+    return [zor, zoo, zoa]
+
+
+
+# update Velocity a la couzin
+def updateV_couzin(agent, zones):
+    dx=0
+    dy=0
+
+    # zor not empty
+    if zones[0] != []:
+        for neighbour in zones[0]:
+            disX = neighbour.point.getX() - agent.point.getX()
+            disY = neighbour.point.getY() - agent.point.getY()
+
+            rX = disX / absvec(disX, disY)
+            rY = disY / absvec(disX, disY)
+
+            dx += rX / absvec(rX, rY)
+            dy += rY / absvec(rX, rY)
+
+        dx = -dx
+        dy = -dy
+
+    # zoo not empty; zoa empty
+    elif zones[1] != [] and zones[2] == []:
+        for neighbour in zones[1]:
+            dx += neighbour.x_velocity / absvec(neighbour.x_velocity, neighbour.y_velocity)
+            dy += neighbour.y_velocity / absvec(neighbour.x_velocity, neighbour.y_velocity)
+            
+        dx += agent.x_velocity / absvec(agent.x_velocity, agent.y_velocity)
+        dy += agent.y_velocity / absvec(agent.x_velocity, agent.y_velocity)
+
+    # zoo empty; zoa not empty
+    elif zones[1] == []  and zones[2] != []:
+        for neighbour in zones[2]:
+            disX = neighbour.point.getX() - agent.point.getX()
+            disY = neighbour.point.getY() - agent.point.getY()
+
+            rX = disX / absvec(disX, disY)
+            rY = disY / absvec(disX, disY)
+
+            dx += rX / absvec(rX, rY)
+            dy += rY / absvec(rX, rY)
+
+    # zoo and zoa not empty
+    elif zones[1] != []  and zones[2] != []:
+        for neighbour in zones[1]:
+            dx += neighbour.x_velocity / absvec(neighbour.x_velocity, neighbour.y_velocity)
+            dy += neighbour.y_velocity / absvec(neighbour.x_velocity, neighbour.y_velocity)
+            
+        dx += agent.x_velocity / absvec(agent.x_velocity, agent.y_velocity)
+        dy += agent.y_velocity / absvec(agent.x_velocity, agent.y_velocity)
+
+        for neighbour in zones[2]:
+            disX = neighbour.point.getX() - agent.point.getX()
+            disY = neighbour.point.getY() - agent.point.getY()
+
+            rX = disX / absvec(disX, disY)
+            rY = disY / absvec(disX, disY)
+
+            dx += rX / absvec(rX, rY)
+            dy += rY / absvec(rX, rY)
+
+        dx = 0.5 * dx
+        dy = 0.5 * dy
+
+    # all zones empty
+    else:
+        dx = agent.x_velocity
+        dy = agent.y_velocity
+
+    # randomness factor / error
+    dx += random.uniform(-1, 1)
+    dy += random.uniform(-1, 1)
+
+    return [dx, dy]
+
+# update function
+def update_couzin(agent, agents, free):
+        # Velocity update
+        for agent  in agents:
+            neigh_matrix = neigbour_in_zones(agent, agents)
+            agent.set_temp_velocity(updateV_couzin(agent,neigh_matrix))
+            
+        # move, draw
+        for agent in agents:
+            # check if rotation is in range of maxTurn
+            agent = move_agent_couzin(agent, maxTurn, radTurn, negRadTurn)
+
+            # normalise direction vector to 1, and multiply by constant speed
+            x = agent.x_velocity/absvec(agent.x_velocity, agent.y_velocity)  * agent.speed
+            agent.y_velocity = agent.y_velocity/absvec(agent.x_velocity, agent.y_velocity)  * agent.speed
+            agent.x_velocity = x
+            
+            if free:
+                agent = checkBoundary_free(agent, winWidth, winHeight)
+            else:
+                agent = checkBoundary(agent, winWidth, winHeight)
+
+            agent.draw_Line()
+        return agents
+
+
+# --------------------------------------------------------------
+# main
+
+def main():
+    
+    agents = []
+
+    # generate agents    
+    for i in range(agentNum):
+        
+        agent = Agent(Point(random.uniform(0,winWidth), random.uniform(0,winHeight)), window, maxV)
+        agent.speed = speed
+        # give them a random starting direction
+        agent.set_xvelocity(random.uniform(-2,2))
+        agent.set_yvelocity(random.uniform(-2,2))
+        agent.draw_Line()
+
+        agents.append(agent)
+
+    # check which simulation to run
+    if couzin:
+        # main loop couzin
+        for i in range(maxTime):
+            agents = update_couzin(agent, agents, free)
+            time.sleep(0.01)
+
+    else:
+        # main loop vicsek (= couzin with specific parameters)
+        for i in range(maxTime):
+            agents = update_couzin(agent, agents, free)
+            time.sleep(0.01)
+
+    window.getMouse()
+    window.close()
+            
+
+main()