1 module Simulate where 2 3 import Array 4 import List 5 import Ratio 6 7 import RandomFix 8 import LinearAlgebra 9 import Types 10 11 simulate :: Circuit -> Seed -> Time -> Time -> Temperature -> RBC -> [(Time, Output)] 12 simulate circuit seed dt end_time temp rbc = select trace 13 where 14 select = takeWhile (\(t, _) -> t <= end_time) . map (\(t, s, o, r) -> (t, o)) 15 trace = (concat . iterate step) ((tunnel . correct . solve . sources) (0, initial', initial', randoms)) 16 step = tunnel . correct . solve . sources . integrate . last 17 randoms = map ((/ (2 ^ 30)) . fromInteger) (random (0, 2 ^ 30 - 1) seed) 18 19 e = 1.6021773e-19 20 k = 1.38066e-23 21 dt' = circa dt 22 23 (initial', solve', integrate', correct_state', backdate') = prepare circuit dt 24 25 solve (t, s, o, r) = (t, s, apply solve' s, r) 26 integrate (t, s, o, r) = (t + dt, apply integrate' s, o, r) 27 correct (t, s, o, r) = (t, apply correct_state' s, v_sub o (apply backdate' s), r) 28 29 sources (t, s, o, r) = (t, random_offsets // map elements circuit, o, drop size r) 30 where 31 elements (i, p, m, Conductor _ ) = (i, 0) 32 elements (i, p, m, Resistor _ ) = (i, 0) 33 elements (i, p, m, Capacitor _ ) = (i, s ! i) 34 elements (i, p, m, Inductor _ ) = (i, s ! i) 35 elements (i, p, m, Vsource list ) = (i, source list) 36 elements (i, p, m, Isource list ) = (i, source list) 37 elements (i, p, m, Junction _ ) = (i, s ! i) 38 39 random_offsets = listArray (1, size) (map ((/ dt') . (* e). (* rbc) . (* 2) . (+ (-0.5))) r) 40 (1, size) = bounds s 41 42 t' = circa t 43 44 source [(t0, v0)] = v0 45 source ((t0, v0) : _) | t' < t0 = v0 46 source ((t0, v0) : (t1, v1) : rest) | t0 >= t1 = error "sources must have increasing times" 47 | t' < t1 = v0 + (v1 - v0) * (t' - t0) / (t1 - t0) 48 | otherwise = source ((t1, v1) : rest) 49 50 tunnel (t, s, o, r) = select (zip cumulative_probabilities new_states) 51 where 52 e2 = e * e 53 kT = k * temp 54 55 e0 = energy_state s 56 energy_state s = sum (map (energy_state' s ) circuit) 57 energy_source s o = sum (map (energy_source' s o) circuit) 58 59 events = (concat . map events') circuit 60 new_states = map (\(_,s,_,_) -> s) events 61 delta_energies = map ((\(r,s,o,_) -> (r, energy_state s - e0 + energy_source s o)) . correct) events 62 probabilities = limit_total (map (no_negative . probability) delta_energies) 63 probability (r, de) | r == 0 = error "tunnel resistance must be nonzero (and should be >> 26k anyway)" 64 | temp == 0 = dt' * (-de) / (e2 * circa r) 65 | de == 0 = dt' * kT / (e2 * circa r) 66 | otherwise = dt' * (-de) / (e2 * circa r * (1 - exp(de / kT))) 67 no_negative p = if p < 0 then 0 else p 68 limit_total ps = let total = sum ps in if total > 1 then map (/ total) ps else ps 69 cumulative_probabilities = tail (scanl (+) 0 probabilities) 70 71 select [] = [(t, s, o, r')] 72 select ((p', s') : zs) | p < p' = [(t, s, o, r'), (correct . solve) (t, s', o, r')] 73 | otherwise = select zs 74 75 (p : r') = r 76 77 energy_state' s (i, p, m, Capacitor (c, _) ) | c == 0 = error "capacitance must be nonzero" 78 | otherwise = 1/2 * (s ! i) ^ 2 / circa c 79 energy_state' s (i, p, m, Junction (c, _, _)) | c == 0 = error "junction capacitance must be nonzero" 80 | otherwise = 1/2 * (s ! i) ^ 2 / circa c 81 energy_state' _ _ = 0 82 83 energy_source' s o (i, p, m, Vsource _ ) = (s ! i) * (o ! i) * dt' 84 energy_source' _ _ _ = 0 85 86 events' (i, p, m, Junction (c, r, _)) = [(r, s // [(i, s ! i + e)], v_map (const 0) o, []), 87 (r, s // [(i, s ! i - e)], v_map (const 0) o, [])] 88 events' _ = [] 89 90 prepare :: Circuit -> Time -> (Vector Approx, Matrix Approx, Matrix Approx, Matrix Approx, Matrix Approx) 91 prepare circuit dt = (v_map circa initial, m_map circa solve, m_map circa integrate, m_map circa correct_state, m_map circa backdate) 92 where 93 elements = matrix element 94 derivative = matrix deriv 95 initial = vector init 96 97 (l1, u1, i1) = gauss_jordan elements 98 (l2, u2, i2) = gauss_jordan (m_transpose elements) 99 constraint = m_select_rows l1 i1 100 freedom = m_transpose (m_select_rows l2 i2) 101 102 (l3, u3, i3) = gauss_jordan (m_mul constraint (m_mul derivative freedom)) 103 inverse3 | l3 == [] = i3 104 | otherwise = error "illegal circuit" 105 correct | l1 == [] = m_zero size 106 | otherwise = m_mul freedom (m_mul inverse3 constraint) 107 correct_state = m_sub (m_unit size) (m_mul derivative correct ) 108 correct_output = m_sub (m_unit size) (m_mul correct derivative) 109 110 solve = m_mul correct_output (m_mul i1 correct_state) 111 112 scaled_derivative = m_map (* (dt / 2)) (m_mul derivative solve) 113 (l4, u4, i4) = gauss_jordan (m_sub (m_unit size) scaled_derivative) 114 integrate | l4 == [] = m_mul i4 (m_add (m_unit size) scaled_derivative) 115 | otherwise = error "unlucky choice of dt - try a smaller one" 116 117 backdate = m_map (/ dt) correct 118 119 matrix f = accumArray (+) (0 % 1) ((1, 1), (size, size)) ((filter no_ground . concat . map f) circuit) 120 vector f = accumArray (+) (0 % 1) ( 1 , size ) (( concat . map f) circuit) 121 122 element (i, p, m, Conductor g ) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, m), -g ), ((i, p), g ), ((i, i), -1 % 1)] 123 element (i, p, m, Resistor r ) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, m), -1 % 1), ((i, p), 1 % 1), ((i, i), -r )] 124 element (i, p, m, Capacitor (c, _) ) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, m), -c ), ((i, p), c )] 125 element (i, p, m, Inductor (l, _) ) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, i), l )] 126 element (i, p, m, Vsource _ ) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, m), -1 % 1), ((i, p), 1 % 1)] 127 element (i, p, m, Isource _ ) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, i), 1 % 1)] 128 element (i, p, m, Junction (c, _, _)) = [((p, i), -1 % 1), ((m, i), 1 % 1), ((i, m), -c ), ((i, p), c )] 129 130 deriv (i, p, m, Conductor _ ) = [] 131 deriv (i, p, m, Resistor _ ) = [] 132 deriv (i, p, m, Capacitor (_, _) ) = [((i, i), 1)] 133 deriv (i, p, m, Inductor (_, _) ) = [((i, p), 1), ((i, m), -1)] 134 deriv (i, p, m, Vsource _ ) = [] 135 deriv (i, p, m, Isource _ ) = [] 136 deriv (i, p, m, Junction (_, _, _)) = [((i, i), 1)] 137 138 init (i, p, m, Conductor _ ) = [] 139 init (i, p, m, Resistor _ ) = [] 140 init (i, p, m, Capacitor (_, q) ) = [(i, q)] 141 init (i, p, m, Inductor (_, f) ) = [(i, f)] 142 init (i, p, m, Vsource _ ) = [] 143 init (i, p, m, Isource _ ) = [] 144 init (i, p, m, Junction (_, _, q)) = [(i, q)] 145 146 no_ground ((r, c), _) = r /= 0 && c /= 0 147 size = maximum (concat (map (\(p, m, i, _) -> [p, m, i]) circuit))