diff --git a/clock_disparity/distributions.py b/clock_disparity/distributions.py
new file mode 100644
index 0000000..6412e3d
--- /dev/null
+++ b/clock_disparity/distributions.py
@@ -0,0 +1,37 @@
+import random, sys
+
+
+def normal_distribution(mean, standev):
+    def f():
+        return int(random.normalvariate(mean, standev))
+
+    return f
+
+
+def exponential_distribution(mean):
+    def f():
+        total = 0
+        while 1:
+            total += 1
+            if not random.randrange(32):
+                break
+        return int(total * 0.03125 * mean)
+
+    return f
+
+
+def convolve(*args):
+    def f():
+        total = 0
+        for arg in args:
+            total += arg()
+        return total
+
+    return f
+
+
+def transform(dist, xformer):
+    def f():
+        return xformer(dist())
+
+    return f
diff --git a/clock_disparity/networksim.py b/clock_disparity/networksim.py
new file mode 100644
index 0000000..2415262
--- /dev/null
+++ b/clock_disparity/networksim.py
@@ -0,0 +1,75 @@
+from distributions import transform, normal_distribution
+import random
+
+
+class NetworkSimulator():
+
+    def __init__(self, latency=50):
+        self.agents = []
+        self.latency_distribution_sample = transform(normal_distribution(latency, (latency * 2) // 5), lambda x: max(x, 0))
+        self.time = 0
+        self.objqueue = {}
+        self.peers = {}
+        self.reliability = 0.9
+
+    def generate_peers(self, num_peers=5):
+        self.peers = {}
+        for a in self.agents:
+            p = []
+            while len(p) <= num_peers // 2:
+                p.append(random.choice(self.agents))
+                if p[-1] == a:
+                    p.pop()
+            self.peers[a.id] = self.peers.get(a.id, []) + p
+            for peer in p:
+                self.peers[peer.id] = self.peers.get(peer.id, []) + [a]
+
+    def tick(self):
+        if self.time in self.objqueue:
+            for recipient, obj in self.objqueue[self.time]:
+                if random.random() < self.reliability:
+                    recipient.on_receive(obj)
+            del self.objqueue[self.time]
+        for a in self.agents:
+            a.tick()
+        self.time += 1
+
+    def run(self, steps):
+        for i in range(steps):
+            self.tick()
+
+    def broadcast(self, sender, obj):
+        for p in self.peers[sender.id]:
+            recv_time = self.time + self.latency_distribution_sample()
+            if recv_time not in self.objqueue:
+                self.objqueue[recv_time] = []
+            self.objqueue[recv_time].append((p, obj))
+
+    def direct_send(self, to_id, obj):
+        for a in self.agents:
+            if a.id == to_id:
+                recv_time = self.time + self.latency_distribution_sample()
+                if recv_time not in self.objqueue:
+                    self.objqueue[recv_time] = []
+                self.objqueue[recv_time].append((a, obj))
+
+    def knock_offline_random(self, n):
+        ko = {}
+        while len(ko) < n:
+            c = random.choice(self.agents)
+            ko[c.id] = c
+        for c in ko.values():
+            self.peers[c.id] = []
+        for a in self.agents:
+            self.peers[a.id] = [x for x in self.peers[a.id] if x.id not in ko]
+
+    def partition(self):
+        a = {}
+        while len(a) < len(self.agents) / 2:
+            c = random.choice(self.agents)
+            a[c.id] = c
+        for c in self.agents:
+            if c.id in a:
+                self.peers[c.id] = [x for x in self.peers[c.id] if x.id in a]
+            else:
+                self.peers[c.id] = [x for x in self.peers[c.id] if x.id not in a]
diff --git a/clock_disparity/node.py b/clock_disparity/node.py
new file mode 100644
index 0000000..5ca4a7b
--- /dev/null
+++ b/clock_disparity/node.py
@@ -0,0 +1,151 @@
+import os
+from binascii import hexlify
+from Crypto.Hash import keccak
+import random
+
+def to_hex(s):
+    return hexlify(s).decode('utf-8')
+
+memo = {}
+
+def sha3(x):
+    if x not in memo:
+        memo[x] = keccak.new(digest_bits=256, data=x).digest()
+    return memo[x]
+
+def hash_to_int(h):
+    o = 0
+    for c in h:
+        o = (o << 8) + c
+    return o
+
+NOTARIES = 40
+BASE_TS_DIFF = 1
+SKIP_TS_DIFF = 6
+SAMPLE = 9
+MIN_SAMPLE = 5
+POWDIFF = 50 * NOTARIES
+SHARDS = 12
+
+def checkpow(work, nonce):
+    # Discrete log PoW, lolz
+    # Quadratic nonresidues only
+    return pow(work, nonce, 65537) * POWDIFF < 65537 * 2 and pow(nonce, 32768, 65537) == 65536
+
+class MainChainBlock():
+    def __init__(self, parent, pownonce, ts):
+        self.parent_hash = parent.hash if parent else (b'\x00' * 32)
+        assert isinstance(self.parent_hash, bytes)
+        self.hash = sha3(self.parent_hash + str(pownonce).encode('utf-8'))
+        self.ts = ts
+        if parent:
+            assert checkpow(parent.pownonce, pownonce)
+            assert self.ts >= parent.ts
+        self.pownonce = pownonce
+        self.number = 0 if parent is None else parent.number + 1
+            
+
+main_genesis = MainChainBlock(None, 59049, 0)
+
+class Node():
+
+    def __init__(self, _id, network, sleepy=False, careless=False, ts=0):
+        self.blocks = {
+            main_genesis.hash: main_genesis
+        }
+        self.main_chain = [main_genesis.hash]
+        self.timequeue = []
+        self.parentqueue = {}
+        self.children = {}
+        self.ts = ts
+        self.id = _id
+        self.network = network
+        self.used_parents = {}
+        self.processed = {}
+        self.sleepy = sleepy
+        self.careless = careless
+
+    def broadcast(self, x):
+        if self.sleepy and self.ts:
+            return
+        self.network.broadcast(self, x)
+        self.on_receive(x)
+
+    def log(self, words, lvl=3, all=False):
+        #if "Tick:" != words[:5] or self.id == 0:
+        if (self.id == 0 or all) and lvl >= 2:
+            print(self.id, words)
+
+    def on_receive(self, obj, reprocess=False):
+        if obj.hash in self.processed and not reprocess:
+            return
+        self.processed[obj.hash] = True
+        if isinstance(obj, MainChainBlock):
+            return self.on_receive_main_block(obj)
+
+    def add_to_timequeue(self, obj):
+        i = 0
+        while i < len(self.timequeue) and self.timequeue[i].ts < obj.ts: 
+            i += 1
+        self.timequeue.insert(i, obj)
+
+    def add_to_multiset(self, _set, k, v):
+        if k not in _set:
+            _set[k] = []
+        _set[k].append(v)
+
+    def change_head(self, chain, new_head):
+        chain.extend([None] * (new_head.number + 1 - len(chain)))
+        i, c = new_head.number, new_head.hash
+        while c != chain[i]:
+            chain[i] = c
+            c = self.blocks[c].parent_hash
+            i -= 1
+        for i in range(len(chain)):
+            assert self.blocks[chain[i]].number == i
+            assert self.blocks[chain[i]].ts <= self.ts
+
+    def process_children(self, h):
+        if h in self.parentqueue:
+            for b in self.parentqueue[h]:
+                self.on_receive(b, reprocess=True)
+            del self.parentqueue[h]
+
+    def on_receive_main_block(self, block):
+        # Parent not yet received
+        if block.parent_hash not in self.blocks:
+            self.add_to_multiset(self.parentqueue, block.parent_hash, block)
+            return None
+        if block.ts > self.ts:
+            self.add_to_timequeue(block)
+            return None
+        self.log("Processing main chain block %s" % to_hex(block.hash[:4]))
+        self.blocks[block.hash] = block
+        # Reorg the main chain if new head
+        if block.number > self.blocks[self.main_chain[-1]].number:
+            reorging = (block.parent_hash != self.main_chain[-1])
+            self.change_head(self.main_chain, block)
+        # Add child record
+        self.add_to_multiset(self.children, block.parent_hash, block.hash)
+        # Final steps
+        self.process_children(block.hash)
+        self.network.broadcast(self, block)
+
+    def is_descendant(self, a, b):
+        a, b = self.blocks[a], self.blocks[b]
+        while b.number > a.number:
+            b = self.blocks[b.parent_hash]
+        return a.hash == b.hash
+
+    def tick(self):
+        self.ts += 0.1
+        self.log("Tick: %.1f" % self.ts, lvl=1)
+        # Process time queue
+        while len(self.timequeue) and self.timequeue[0].ts <= self.ts:
+            self.on_receive(self.timequeue.pop(0))
+        # Attempt to mine a main chain block
+        pownonce = random.randrange(65537)
+        mchead = self.blocks[self.main_chain[-1]]
+        if checkpow(mchead.pownonce, pownonce):
+            assert self.ts >= mchead.ts
+            self.broadcast(MainChainBlock(mchead, pownonce, self.ts))
diff --git a/clock_disparity/test.py b/clock_disparity/test.py
new file mode 100644
index 0000000..a942e9a
--- /dev/null
+++ b/clock_disparity/test.py
@@ -0,0 +1,46 @@
+from networksim import NetworkSimulator
+from node import Node, NOTARIES, MainChainBlock, main_genesis
+from distributions import normal_distribution
+
+net = NetworkSimulator(latency=12)
+notaries = [Node(i, net, ts=max(normal_distribution(50, 50)(), 0)) for i in range(NOTARIES)]
+net.agents = notaries
+net.generate_peers()
+for i in range(4000):
+    net.tick()
+for n in notaries:
+    print("Main chain head: %d" % n.blocks[n.main_chain[-1]].number)
+    print("Total main chain blocks received: %d" % (len([b for b in n.blocks.values() if isinstance(b, MainChainBlock)]) - 1))
+
+import matplotlib.pyplot as plt
+import networkx as nx
+import random
+
+G=nx.Graph()
+
+#positions = {main_genesis.hash: 0, beacon_genesis.hash: 0}
+#queue = [
+
+for b in n.blocks.values():
+    if b.number > 0:
+        if isinstance(b, MainChainBlock):
+            G.add_edge(b.hash, b.parent_hash, color='b')
+
+cache = {main_genesis.hash: 0}
+
+def mkoffset(b):
+    if b.hash not in cache:
+        cache[b.hash] = cache[b.parent_hash] + random.randrange(35)
+    return cache[b.hash]
+        
+pos={b.hash: (b.ts + mkoffset(b), b.ts) for b in n.blocks.values()}
+edges = G.edges()
+colors = [G[u][v]['color'] for u,v in edges]
+nx.draw_networkx_nodes(G,pos,node_size=10,node_shape='o',node_color='0.75')
+
+nx.draw_networkx_edges(G,pos,
+                width=2,edge_color=colors)
+
+plt.axis('off')
+# plt.savefig("degree.png", bbox_inches="tight")
+plt.show() 
diff --git a/randao_analysis/paths.py b/randao_analysis/paths.py
index dd8aa8f..b7bb4a3 100644
--- a/randao_analysis/paths.py
+++ b/randao_analysis/paths.py
@@ -7,7 +7,7 @@ SKIP = 6
 # The score of a node in the search graph, for A* search
 def score_node(n):
     time, height = n
-    return time - height * (LATENCY + 0.0000001)
+    return time - height * (LATENCY * 5.0 + 0.0000001)
 
 # `alpha`: percentage of the network controlled by the actor
 # `max_height`: stop at this height
@@ -120,7 +120,7 @@ def test():
     # Assume attacker has (0.40, 0.38, 0.36, 0.34, 0.32, 0.30, 0.28); check race probs for handicap 1-8
     print("Basic tests\n")
     for alpha in (0.40, 0.38, 0.36, 0.34, 0.32, 0.30, 0.28):
-        race_results = [race(alpha, 1-alpha, 100) for j in range(100)]
+        race_results = [race(alpha, 1-alpha, 1000) for j in range(100)]
         probs = []
         for i in range(1, 9):
             probs.append(len([x for x in race_results if x >= i]) / 100)
@@ -130,7 +130,7 @@ def test():
     # if either side has share p, they effectively "really" have share p**3
     print("\nRequiring 2/2 notarization\n")
     for alpha in (0.45, 0.44, 0.43, 0.42, 0.41, 0.40, 0.39):
-        race_results = [race(alpha ** 3, (1-alpha) ** 3, 100) for j in range(100)]
+        race_results = [race(alpha ** 3, (1-alpha) ** 3, 1000) for j in range(100)]
         probs = []
         for i in range(1, 9):
             probs.append(len([x for x in race_results if x >= i]) / 100)