diff --git a/src/main.moon b/src/main.moon
index 90d4fc7..db6dcc3 100644
--- a/src/main.moon
+++ b/src/main.moon
@@ -23,104 +23,112 @@ class Isotopes
     return @randoms[i]
   generate: (z, n) =>
     @[z] = {} unless @[z]
-    isotope = {}
+    isotope = { :z, :n }
 
-    -- isotopes only exist in this range
-    if n/z > 0.9 and n/z < 1.2
-      v = 2 / math.abs(n - z)        -- base stability is proximity to n == z
-      -- v /= @rng\randomNormal 0.25, 1 -- a slight random variation
-      v /= @normal(z, n)
-      local odd_even                 -- based on statistics from real isotopes
-      if z % 2 == 0
-        if n % 2 == 0
-          odd_even = 0.6
-        else
-          odd_even = 0.205
-      else
-        if n % 2 == 0
-          odd_even = 0.18
-        else
-          odd_even = 0.015
-      -- if @rng\random! <= v * odd_even -- more random variation based on stats
-      if @random(z, n) <= v * odd_even
-        v *= odd_even
-      if v > 1                        -- maximum stability should be close to 1
-        v = 1.05
-      v -= z / 1000           -- lighter elements should have a higher stability
-      isotope.stability = v
+    if n/z > 0.88 and n/z < 1.15
+      lower_const = 16.5 -- aiming for 1e-16 (0.1 femtoseconds) minimum half life
+      upper_const = 1e18 -- aiming for 1e17 (3 billion years) maximum half life (besides pure stability)
+      v = math.abs(z - n)^lower_const / upper_const
+      v = math.abs v / @normal(z, n)
+      odd_even = { [0]: { [0]: 0.6, 0.205 }, { [0]: 0.18, 0.015 } }
+      c = odd_even[z % 2][n % 2]
+      if @random(z, n) <= v * c
+        v *= c
+      v = math.abs v - z / 1e17
+      if @random(z, n) < 0.5
+        v *= 2
+      -- if v < 1 and @normal(z, n) < z / 200
+      --   v = v^0.001
+      if v < 1 and z + n > @normal(z, n) * 280
+        v = v^0.3
+      -- if v < 1 and z + n > @normal(n, z) * 400
+      --   v = v^0.1
+      isotope.decay = v
     else
-      isotope.stability = 0   -- this isotope cannot exist
+      isotope.decay = math.huge
 
-    if isotope.stability < 1
-      isotope.half_life = 0.69314718056 / (1 - isotope.stability)
-      -- isotope.half_life = 0.69314718056 / isotope.stability
+    if isotope.decay < math.huge
+      isotope.half_life = 0.69314718056 / isotope.decay
+      -- print isotope.half_life
     else
-      isotope.half_life = math.huge
+      isotope.half_life = 0
+
+    m = z + n
+    e = 14 - 13 / m^(1/3) - 0.585 * z^2 / m^(4/3) - 19.3 * (n - z)^2 / m^2
+    if z % 2 != 0 and n % 2 != 0
+      e -= 33 / m^(7/4)
+    else
+      e += 33 / m^(7/4)
+
+    isotope.binding_energy = math.abs e
+    -- verified "accuracy"
+    -- if isotope.half_life > 1 and m < 10
+    --   print math.floor(e * 100)/100, m, z, n
 
     @[z][n] = isotope
 
--- Z  > 82 -> unstable
--- Z <= 20 -> stable n/p is nearly 1 (0.775 - 1.29)
--- Z  > 20  & Z <= 82 -> stable n/p between 1:1 and 1.5:1
--- Z  > 92 -> spontaneous fission
 -- Z  > 52 -> alpha rays possible
 -- Can induce fission by neutron bombardment.
+
+-- Decay types
+-- Spontaneous: Splits in half + extra neutrons, Z > 92
 -- Too many neutrons -> beta decay, neutron emission
 -- Too few neutrons -> electron capture, positron emission, proton emission
 -- Alpha: Emits 2Z-2N
--- Beta: Emits e, 1N -> 1Z
+-- Beta: Emits e-, 1N -> 1Z
 -- Electron capture & Positron emission: Emits e+, 1Z -> 1N
+--  (difference: positron emission uses its own electron, electron capture steals one)
 
 export showUnstable = true
 w, h = love.graphics.getDimensions!
+scale = 6
 
 isotopes = Isotopes!
 
 love.load = ->
-  above_0 = 0
-  above_1 = 0
-  min_half_life = math.huge
-  max_half_life = 0
-  min_z, min_n = 0, 0
-  max_z, max_n = 0, 0
-  min_s = 0
-  max_s = 0
-  for z = 1, w
-    for n = 0, h
+  shortest = { half_life: math.huge }
+  longest = { half_life: -math.huge }
+  heaviest = { z: 0, n: 0 }
+  for z = 1, w / scale
+    for n = 0, h / scale
       isotope = isotopes\get z, n
-      if isotope.stability > 0
-        above_0 += 1
-        if isotope.half_life < min_half_life
-          min_half_life = isotope.half_life
-          min_z = z
-          min_n = n
-          min_s = isotope.stability
-        if isotope.stability < 1
-          if isotope.half_life > max_half_life
-            max_half_life = isotope.half_life
-            max_z = z
-            max_n = n
-            max_s = isotope.stability
-        else
-          above_1 += 1
-  print "#{above_0} isotopes", "#{above_1} stable isotopes"
-  print "Shortest half life:", "#{min_half_life}", "#{min_z}P#{min_n}N", min_s
-  print "Longest half life:", "#{max_half_life}", "#{max_z}P#{max_n}N", max_s
-
--- decay constant (probability of decay per second): 1 - stability
--- half life: 0.69314718056 / decay constant
--- I have redefined stability to equal the decay constant
+      if isotope.half_life > 1e15 and isotope.z + isotope.n > heaviest.z + heaviest.n
+        heaviest = isotope
+      if isotope.half_life < shortest.half_life and isotope.half_life != 0
+        shortest = isotope
+      elseif isotope.half_life > longest.half_life and isotope.half_life != math.huge
+        longest = isotope
+  print "Shortest half life: #{shortest.half_life}", "#{shortest.z}P#{shortest.n}N", "Decay: #{shortest.decay}"
+  print "Longest half life: #{longest.half_life}", "#{longest.z}P#{longest.n}N", "Decay: #{longest.decay}"
+  print "Heaviest stable: #{heaviest.z}P#{heaviest.n}N", "#{heaviest.z + heaviest.n}"
 
 love.draw = ->
-  scale = 3
-  for z = 1, w
-    for n = 0, h
+  for z = 1, w / scale
+    for n = 0, h / scale
       isotope = isotopes\get z, n
-      s = isotope.stability
-      if s >= 1 or showUnstable
-        -- c = 1 - s -- side effect -> stable isotopes are now invisible
-        love.graphics.setColor s, s, s, 1
-        love.graphics.rectangle "fill", (z - 1) * scale, (h / scale - n - 1) * scale, scale, scale
+      -- s = math.log 1 / isotope.decay
+      -- love.graphics.setColor s, s, s, 1
+      colors = {
+        {0, {0, 0, 0, 1}},              -- non-existent, black
+        {0.01, {0.1, 0.1, 0.1, 1}},     -- hundredth and less, dark grey
+        {1, {0, 0, 1, 1}},              -- hundredth to second, blue
+        {60, {0, 1, 1, 1}},             -- second to minute, teal
+        {60*60, {0, 0.5, 0, 1}},        -- minute to hour, dark green
+        {60*60*24, {0.2, 1, 0.2, 1}},   -- hour to day, light green
+        {60*60*24*7, {1, 1, 0, 1}},     -- day to week, yellow
+        {60*60*24*30, {1, 0.5, 0, 1}},  -- week to month, orange
+        {60*60*24*365, {1, 0, 0, 1}},   -- month to year, red
+        {60*60*24*365*100, {1, 0.5, 0.5, 1}}, -- year to century, pink
+        {1e15, {0.8, 0.8, 0.8, 1}}            -- century to 31 million years, light grey
+        {math.huge, {1, 1, 1, 1}}             -- purely stable, white
+      }
+      for color in *colors
+        if isotope.half_life <= color[1]
+          love.graphics.setColor color[2]
+          break
+      -- if z == n -- NOTE temp
+      --   love.graphics.setColor 1, 0, 0, 1
+      love.graphics.rectangle "fill", (z - 1) * scale, (h / scale - n - 1) * scale, scale, scale
 
 love.keypressed = (key) ->
   love.event.quit! if key == "escape"