nbody.mlish (5834B)
1 #lang s-exp "../../mlish.rkt" 2 (require "../rackunit-typechecking.rkt") 3 4 (define +pi+ 3.141592653589793) 5 6 (check-type +pi+ : Float) 7 8 (define +days-per-year+ 365.24) 9 10 (define * fl*) 11 12 (define +solar-mass+ (* 4.0 (* +pi+ +pi+))) 13 14 (check-type +solar-mass+ : Float) 15 16 (define +dt+ 0.01) 17 18 (define-type Body (body Float ; x 19 Float ; y 20 Float ; z 21 Float ; vx 22 Float ; vy 23 Float ; vz 24 Float ; mass 25 )) 26 27 (define *sun* 28 (body 0.0 0.0 0.0 0.0 0.0 0.0 +solar-mass+)) 29 30 (define *jupiter* 31 (body 4.84143144246472090 32 -1.16032004402742839 33 -1.03622044471123109e-1 34 (* 1.66007664274403694e-3 +days-per-year+) 35 (* 7.69901118419740425e-3 +days-per-year+) 36 (* -6.90460016972063023e-5 +days-per-year+) 37 (* 9.54791938424326609e-4 +solar-mass+))) 38 39 (define *saturn* 40 (body 8.34336671824457987 41 4.12479856412430479 42 -4.03523417114321381e-1 43 (* -2.76742510726862411e-3 +days-per-year+) 44 (* 4.99852801234917238e-3 +days-per-year+) 45 (* 2.30417297573763929e-5 +days-per-year+) 46 (* 2.85885980666130812e-4 +solar-mass+))) 47 48 (define *uranus* 49 (body 1.28943695621391310e1 50 -1.51111514016986312e1 51 -2.23307578892655734e-1 52 (* 2.96460137564761618e-03 +days-per-year+) 53 (* 2.37847173959480950e-03 +days-per-year+) 54 (* -2.96589568540237556e-05 +days-per-year+) 55 (* 4.36624404335156298e-05 +solar-mass+))) 56 57 (define *neptune* 58 (body 1.53796971148509165e+01 59 -2.59193146099879641e+01 60 1.79258772950371181e-01 61 (* 2.68067772490389322e-03 +days-per-year+) 62 (* 1.62824170038242295e-03 +days-per-year+) 63 (* -9.51592254519715870e-05 +days-per-year+) 64 (* 5.15138902046611451e-05 +solar-mass+))) 65 66 (define *system* (list *sun* *jupiter* *saturn* *uranus* *neptune*)) 67 68 (check-type *system* : (List Body)) 69 70 (define (offset-momentum -> Unit) 71 (let loop-i : Unit 72 ([i *system*] [px 0.0] [py 0.0] [pz 0.0]) 73 (cond 74 [(isnil i) 75 (match (head *system*) with ; sun 76 [body x y z vx vy vz m -> 77 (let ([new 78 (body x y z 79 (fl/ (fl- 0.0 px) +solar-mass+) 80 (fl/ (fl- 0.0 py) +solar-mass+) 81 (fl/ (fl- 0.0 pz) +solar-mass+) 82 m)]) 83 (set! *system* (cons new (tail *system*))))])] 84 [else 85 (match (head i) with 86 [body x y z vx vy vz m -> 87 (loop-i (tail i) 88 (fl+ px (fl* vx m)) 89 (fl+ py (fl* vy m)) 90 (fl+ pz (fl* vz m)))])]))) 91 92 (define (energy [o : (List Body)] -> Float) 93 (let loop-o : Float ([o o] [e 0.0]) 94 (cond 95 [(isnil o) e] 96 [else 97 (match (head o) with 98 [body x y z vx vy vz m -> 99 (let* ([e (fl+ e (fl* 0.5 100 (fl* m 101 (fl+ (fl+ (fl* vx vx) 102 (fl* vy vy)) 103 (fl* vz vz)))))]) 104 (let loop-i : Float ([i (tail o)] [e e]) 105 (if (isnil i) 106 (loop-o (tail o) e) 107 (match (head i) with 108 [body x2 y2 z2 vx2 vy2 vz2 m2 -> 109 (let* ([dx (fl- x x2)] 110 [dy (fl- y y2)] 111 [dz (fl- z z2)] 112 [dist (flsqrt (fl+ (fl+ (fl* dx dx) (fl* dy dy)) 113 (fl* dz dz)))] 114 [e (fl- e (fl/ (fl* m m2) dist))]) 115 (loop-i (tail i) e))]))))])]))) 116 117 (define (advance [bs : (List Body)] -> (List Body)) 118 (if (isnil bs) 119 bs 120 (let ([new (advance2 bs)]) 121 (cons (head new) (advance (tail new)))))) 122 ;; bs is non-null 123 (define (advance2 [bs : (List Body)] -> (List Body)) 124 (match (head bs) with 125 [body o1x o1y o1z vx vy vz om -> 126 (let loop-i : (List Body) 127 ([i (tail bs)] [res (nil {Body})] [vx vx] [vy vy] [vz vz]) 128 (cond 129 [(isnil i) 130 (cons 131 (body 132 (fl+ o1x (fl* +dt+ vx)) 133 (fl+ o1y (fl* +dt+ vy)) 134 (fl+ o1z (fl* +dt+ vz)) 135 vx vy vz om) 136 (reverse res))] 137 [else 138 (match (head i) with 139 [body i1x i1y i1z i1vx i1vy i1vz im -> 140 (let* ([dx (fl- o1x i1x)] 141 [dy (fl- o1y i1y)] 142 [dz (fl- o1z i1z)] 143 [dist2 (fl+ (fl+ (fl* dx dx) (fl* dy dy)) (fl* dz dz))] 144 [mag (fl/ +dt+ (fl* dist2 (flsqrt dist2)))] 145 [dxmag (fl* dx mag)] 146 [dymag (fl* dy mag)] 147 [dzmag (fl* dz mag)]) 148 (loop-i 149 (tail i) 150 (cons (body i1x i1y i1z 151 (fl+ i1vx (fl* dxmag om)) 152 (fl+ i1vy (fl* dymag om)) 153 (fl+ i1vz (fl* dzmag om)) 154 im) res) 155 (fl- vx (fl* dxmag im)) 156 (fl- vy (fl* dymag im)) 157 (fl- vz (fl* dzmag im))))])]))])) 158 159 (check-type (real->decimal-string (energy *system*) 9) 160 : String -> "-0.169289903") 161 162 (offset-momentum) 163 164 (check-type (real->decimal-string (energy *system*) 9) 165 : String -> "-0.169075164") 166 167 (check-type 168 (real->decimal-string 169 (energy (advance *system*)) 170 9) 171 : String -> "-0.169074954") 172 173 (check-type 174 (real->decimal-string 175 (energy (advance (advance *system*))) 9) 176 : String -> "-0.169074743") 177 178 (check-type 179 (real->decimal-string 180 (energy 181 (for/fold ([bs *system*]) 182 ([i (in-range 10)]) 183 (advance bs))) 184 9) 185 : String -> "-0.169073022")