1 {- A kernel fragment from a program written by 2 Ron Legere -- http://www.its.caltech.edu/~legere 3 Caltech Quantum Optics 4 5 It has the interesting property that Classic Hugs 6 runs it 20x faster than GHC! 7 Reason: runExperiment calls itself with identical parameters, 8 and Hugs commons that up for some reason. 9 10 (Even with that fixed, STG Hugs ran the program a lot 11 slower than Classic Hugs.) 12 13 So it seems like an interesting program. It appears here 14 in the form with the silly self-call, because that makes 15 it run a nice long time. It thrashes floating point 16 multiplication and lists. 17 -} 18 19 module Main where 20 21 import System 22 23 infixl 9 .* 24 infix 9 <*> 25 26 main = do 27 [arg] <- getArgs 28 let n = read arg :: Int 29 putStr (show (take n test)) 30 31 test :: StateStream 32 test = runExperiment testforce 0.02 [1.0] (State [1.0] [0.0]) 33 34 35 testforce :: ForceLaw [Float] 36 testforce k [] = [] 37 testforce k ( (State pos vel):atoms) = (-1.0) .* k * pos: 38 testforce k atoms 39 40 {- 41 The test force: K is a list of spring 42 constants. (But here I am only doing one dimension for the purposes 43 of demonstrating the approach) 44 -} 45 46 47 48 {- 49 ****************** 50 ****************** 51 Module: Numerical classical atom (atom.lhs) 52 ****************** 53 ****************** 54 -} 55 56 57 -- We will want types for the whole simulation (where we can configure 58 -- the dimensions, etc), for the results (a state stream), and the force laws. 59 60 data AtomState = State Position Velocity 61 type Position = [Float] 62 type Velocity = [Float] 63 type Force = [Float] 64 65 66 type StateStream = [AtomState] 67 68 69 70 {- 71 I made AtomState a data type, just so I could play with them. I think 72 I would prefer to keep it just a synonym, because that would 73 be simpler! 74 75 Now we need a function to write out the results to a file in a nice format. 76 I think I would prefer a simple x y z /n x y z /n etc 77 78 NOTE that show AtomState only shows the position! 79 -} 80 81 instance Show AtomState where 82 show (State pos vel) = concat [ (show component) ++ "\t" | component <- pos ] 83 showList states = showString (concat [(show state) ++ "\n" | state <- states]) 84 85 86 {- 87 Note that I used lists for the position and velocity to allow for 88 unknown number of dimensions. I suspect this will have to 89 be optimized into tuples at some point! 90 91 92 93 Ok, so how shall we define the ForceLaw type? 94 -} 95 96 type ForceLaw a = a -> StateStream -> [Force] 97 98 99 100 {- 101 The force law maps a stream of states to a stream of forces so that time 102 dependant forces can be used. The parametric type 'a' is to allow the force law 103 to depend on some parameter, for example (the common case!) a seed for a random number 104 generater, and/or the timestep, or the spring constant 105 -} 106 107 runExperiment :: ForceLaw a -> Float -> a -> AtomState -> StateStream 108 109 {- 110 In this form this program takes 1 min when compiled under ghc-4.05, 111 but takes 3 seconds under hugs.... 112 -} 113 runExperiment law dt param init = init : zipWith (propagate dt) 114 (law param stream) 115 stream 116 where stream = 117 runExperiment law dt param init 118 119 {- 120 -- In this form GHC is (as expected) much faster 121 runExperiment law dt param init = stream 122 where 123 stream = init : zipWith (propagate dt) 124 (law param stream) 125 stream 126 -} 127 128 {- 129 runExperiment forces timestep param initialcondition :: [AtomState] 130 is an infinite stream of atom states. We can then use this to 131 generate necessary averages, temperatures, allen variences , or wtf 132 you want. 133 134 We could for example, start the random number generator with seed param, if 135 the type is int . 136 137 It is an error to have the initial 138 atom state not have the correct number of dimensions. 139 -} 140 141 propagate :: Float -> Force -> AtomState -> AtomState 142 143 {- 144 Ok, I see one problem with this, not general enough! Some better propagators 145 exist that can use previous atom states. Actually, by using previous atom states, 146 we will not even need to seperately track the velocities either. Oh well, for now 147 I will stick with that. 148 -} 149 150 propagate dt aforce (State pos vel) = State newpos newvel 151 where newpos = pos + (dt .* vel) 152 newvel = vel + (dt .* aforce) 153 154 -- Note assumes mass =1 155 156 157 158 {- 159 ******************************************************** 160 ******************************************************** 161 162 Numerical Lists 163 164 ******************************************************** 165 ******************************************************** 166 -} 167 168 instance Num a => Num [a] where 169 negate (f:fs) = (negate f):(negate fs) 170 l + [] = l 171 [] + l = l 172 (f:fs) + (g:gs) = (f+g):fs+gs 173 _ * [] = [] 174 [] * _ = [] 175 (f:fs) * (g:gs) = (f*g):(gs*gs) 176 fromInteger c = fromInteger c : [0] 177 178 179 (.*):: Num a => a-> [a] -> [a] 180 c .* [] = [] 181 182 183 c .* (f:fs) = c*f : c .* fs 184 185 (<*>):: Num a => [a] -> [a] -> a 186 f <*> g = sum (f*g) 187 188