Tangent 4 years ago
parent
commit
5ff3c2f57b
1 changed files with 82 additions and 66 deletions
  1. 82
    66
      src/main.moon

+ 82
- 66
src/main.moon View File

@@ -23,104 +23,112 @@ class Isotopes
23 23
     return @randoms[i]
24 24
   generate: (z, n) =>
25 25
     @[z] = {} unless @[z]
26
-    isotope = {}
26
+    isotope = { :z, :n }
27 27
 
28
-    -- isotopes only exist in this range
29
-    if n/z > 0.9 and n/z < 1.2
30
-      v = 2 / math.abs(n - z)        -- base stability is proximity to n == z
31
-      -- v /= @rng\randomNormal 0.25, 1 -- a slight random variation
32
-      v /= @normal(z, n)
33
-      local odd_even                 -- based on statistics from real isotopes
34
-      if z % 2 == 0
35
-        if n % 2 == 0
36
-          odd_even = 0.6
37
-        else
38
-          odd_even = 0.205
39
-      else
40
-        if n % 2 == 0
41
-          odd_even = 0.18
42
-        else
43
-          odd_even = 0.015
44
-      -- if @rng\random! <= v * odd_even -- more random variation based on stats
45
-      if @random(z, n) <= v * odd_even
46
-        v *= odd_even
47
-      if v > 1                        -- maximum stability should be close to 1
48
-        v = 1.05
49
-      v -= z / 1000           -- lighter elements should have a higher stability
50
-      isotope.stability = v
28
+    if n/z > 0.88 and n/z < 1.15
29
+      lower_const = 16.5 -- aiming for 1e-16 (0.1 femtoseconds) minimum half life
30
+      upper_const = 1e18 -- aiming for 1e17 (3 billion years) maximum half life (besides pure stability)
31
+      v = math.abs(z - n)^lower_const / upper_const
32
+      v = math.abs v / @normal(z, n)
33
+      odd_even = { [0]: { [0]: 0.6, 0.205 }, { [0]: 0.18, 0.015 } }
34
+      c = odd_even[z % 2][n % 2]
35
+      if @random(z, n) <= v * c
36
+        v *= c
37
+      v = math.abs v - z / 1e17
38
+      if @random(z, n) < 0.5
39
+        v *= 2
40
+      -- if v < 1 and @normal(z, n) < z / 200
41
+      --   v = v^0.001
42
+      if v < 1 and z + n > @normal(z, n) * 280
43
+        v = v^0.3
44
+      -- if v < 1 and z + n > @normal(n, z) * 400
45
+      --   v = v^0.1
46
+      isotope.decay = v
51 47
     else
52
-      isotope.stability = 0   -- this isotope cannot exist
48
+      isotope.decay = math.huge
53 49
 
54
-    if isotope.stability < 1
55
-      isotope.half_life = 0.69314718056 / (1 - isotope.stability)
56
-      -- isotope.half_life = 0.69314718056 / isotope.stability
50
+    if isotope.decay < math.huge
51
+      isotope.half_life = 0.69314718056 / isotope.decay
52
+      -- print isotope.half_life
57 53
     else
58
-      isotope.half_life = math.huge
54
+      isotope.half_life = 0
55
+
56
+    m = z + n
57
+    e = 14 - 13 / m^(1/3) - 0.585 * z^2 / m^(4/3) - 19.3 * (n - z)^2 / m^2
58
+    if z % 2 != 0 and n % 2 != 0
59
+      e -= 33 / m^(7/4)
60
+    else
61
+      e += 33 / m^(7/4)
62
+
63
+    isotope.binding_energy = math.abs e
64
+    -- verified "accuracy"
65
+    -- if isotope.half_life > 1 and m < 10
66
+    --   print math.floor(e * 100)/100, m, z, n
59 67
 
60 68
     @[z][n] = isotope
61 69
 
62 70
 -- Z  > 52 -> alpha rays possible
63 71
 -- Can induce fission by neutron bombardment.
72
+
73
+-- Decay types
74
+-- Spontaneous: Splits in half + extra neutrons, Z > 92
64 75
 -- Too many neutrons -> beta decay, neutron emission
65 76
 -- Too few neutrons -> electron capture, positron emission, proton emission
66 77
 -- Alpha: Emits 2Z-2N
78
+-- Beta: Emits e-, 1N -> 1Z
67 79
 -- Electron capture & Positron emission: Emits e+, 1Z -> 1N
80
+--  (difference: positron emission uses its own electron, electron capture steals one)
68 81
 
69 82
 export showUnstable = true
70 83
 w, h = love.graphics.getDimensions!
84
+scale = 6
71 85
 
72 86
 isotopes = Isotopes!
73 87
 
74 88
 love.load = ->
75
-  above_0 = 0
76
-  above_1 = 0
77
-  min_half_life = math.huge
78
-  max_half_life = 0
79
-  min_z, min_n = 0, 0
80
-  max_z, max_n = 0, 0
81
-  min_s = 0
82
-  max_s = 0
83
-  for z = 1, w
84
-    for n = 0, h
89
+  shortest = { half_life: math.huge }
90
+  longest = { half_life: -math.huge }
91
+  heaviest = { z: 0, n: 0 }
92
+  for z = 1, w / scale
93
+    for n = 0, h / scale
85 94
       isotope = isotopes\get z, n
86
-      if isotope.stability > 0
87
-        above_0 += 1
88
-        if isotope.half_life < min_half_life
89
-          min_half_life = isotope.half_life
90
-          min_z = z
91
-          min_n = n
92
-          min_s = isotope.stability
93
-        if isotope.stability < 1
94
-          if isotope.half_life > max_half_life
95
-            max_half_life = isotope.half_life
96
-            max_z = z
97
-            max_n = n
98
-            max_s = isotope.stability
99
-        else
100
-          above_1 += 1
101
-  print "#{above_0} isotopes", "#{above_1} stable isotopes"
102
-  print "Shortest half life:", "#{min_half_life}", "#{min_z}P#{min_n}N", min_s
103
-  print "Longest half life:", "#{max_half_life}", "#{max_z}P#{max_n}N", max_s
104
-
95
+      if isotope.half_life > 1e15 and isotope.z + isotope.n > heaviest.z + heaviest.n
96
+        heaviest = isotope
97
+      if isotope.half_life < shortest.half_life and isotope.half_life != 0
98
+        shortest = isotope
99
+      elseif isotope.half_life > longest.half_life and isotope.half_life != math.huge
100
+        longest = isotope
101
+  print "Shortest half life: #{shortest.half_life}", "#{shortest.z}P#{shortest.n}N", "Decay: #{shortest.decay}"
102
+  print "Longest half life: #{longest.half_life}", "#{longest.z}P#{longest.n}N", "Decay: #{longest.decay}"
103
+  print "Heaviest stable: #{heaviest.z}P#{heaviest.n}N", "#{heaviest.z + heaviest.n}"
105 104
 
106 105
 love.draw = ->
107
-  scale = 3
108
-  for z = 1, w
109
-    for n = 0, h
106
+  for z = 1, w / scale
107
+    for n = 0, h / scale
110 108
       isotope = isotopes\get z, n
111
-      s = isotope.stability
112
-      if s >= 1 or showUnstable
113
-        -- c = 1 - s -- side effect -> stable isotopes are now invisible
114
-        love.graphics.setColor s, s, s, 1
115
-        love.graphics.rectangle "fill", (z - 1) * scale, (h / scale - n - 1) * scale, scale, scale
109
+      -- s = math.log 1 / isotope.decay
110
+      -- love.graphics.setColor s, s, s, 1
111
+      colors = {
112
+        {0, {0, 0, 0, 1}},              -- non-existent, black
113
+        {0.01, {0.1, 0.1, 0.1, 1}},     -- hundredth and less, dark grey
114
+        {1, {0, 0, 1, 1}},              -- hundredth to second, blue
115
+        {60, {0, 1, 1, 1}},             -- second to minute, teal
116
+        {60*60, {0, 0.5, 0, 1}},        -- minute to hour, dark green
117
+        {60*60*24, {0.2, 1, 0.2, 1}},   -- hour to day, light green
118
+        {60*60*24*7, {1, 1, 0, 1}},     -- day to week, yellow
119
+        {60*60*24*30, {1, 0.5, 0, 1}},  -- week to month, orange
120
+        {60*60*24*365, {1, 0, 0, 1}},   -- month to year, red
121
+        {60*60*24*365*100, {1, 0.5, 0.5, 1}}, -- year to century, pink
122
+        {1e15, {0.8, 0.8, 0.8, 1}}            -- century to 31 million years, light grey
123
+        {math.huge, {1, 1, 1, 1}}             -- purely stable, white
124
+      }
125
+      for color in *colors
126
+        if isotope.half_life <= color[1]
127
+          love.graphics.setColor color[2]
128
+          break
129
+      -- if z == n -- NOTE temp
130
+      --   love.graphics.setColor 1, 0, 0, 1
131
+      love.graphics.rectangle "fill", (z - 1) * scale, (h / scale - n - 1) * scale, scale, scale
116 132
 
117 133
 love.keypressed = (key) ->
118 134
   love.event.quit! if key == "escape"