January 20, 2020
By: Kevin
平面星系模拟
内容
演讲的主旨是c的历史在将在clojure上重演, 从不接受到成为业界标准.
大叔的雄辩和文采值得一观, 推荐他的blog和文章:
以前在Bob的基础上做了个cljs版 , 加了些注释, 放到了这里
代码
位置定义
(ns physics.position "定义位置")
(defrecord position [x y])
(defn make
"构造一个position"
([] (->position 0 0))
([x y] (->position x y)))
(defn origin?
"判断是否处于(0, 0)位置"
[p]
(and (zero? (:x p)) (zero? (:y p))))
(defn add
"position位置矢量相加"
[{x1 :x y1 :y} {x2 :x y2 :y}]
(->position (+ x1 x2) (+ y1 y2)))
(defn subtract
"position位置矢量相减"
[{x1 :x y1 :y} {x2 :x y2 :y}]
(->position (- x1 x2) (- y1 y2)))
(defn distance
"position距离计算"
[{x1 :x y1 :y} {x2 :x y2 :y}]
(Math/sqrt
(+
(Math/pow (- x1 x2) 2)
(Math/pow (- y1 y2) 2))))
(defn mean
"a,b两个数字的平均数"
[a b]
(/ (+ a b) 2))
(defn average
"两个位置的中间点"
[{x1 :x y1 :y} {x2 :x y2 :y}]
(->position (mean x1 x2) (mean y1 y2)))
矢量定义
(ns physics.vector
"vector其实是就是一个泛化的position, 因为position是个矢量"
(:require [physics.position :as position :refer [position]]))
(defn make
([] (position. 0 0))
([x y] (position. x y)))
(defn zero_mag?
"是否是零向量"
[v]
(position/origin? v))
(def add position/add)
(def subtract position/subtract)
(defn scale [v s]
(make
(* (:x v) s)
(* (:y v) s)))
(defn magnitude
"矢量的大小"
[v]
(Math/sqrt
(+
(Math/pow (:x v) 2)
(Math/pow (:y v) 2))))
(defn unit
"单位向量, 保留方向, 大小设置为1"
[v]
(scale v (/ 1 (magnitude v))))
(defn rotate90
"垂直向量, 旋转90度"
[{x :x y :y}]
(make (- y) x))
(defn equal
"向量判等"
[{x1 :x y1 :y} {x2 :x y2 :y}]
(and (== x1 x2) (== y1 y2)))
物理对象定义
(ns physics.star
"模拟实体对象:位置, 质量, 速度, 力, 加速度"
(:refer-clojure :exclude (merge))
(:require-macros [cljs.core.async.macros :refer [go go-loop]])
(:require [physics.position :as position]
[physics.vector :as vector]
[cljs.core.async :refer [chan put! >!]]))
(defrecord star [position
mass
velocity
force
name])
(defn make
([]
(->star (position/make) 0 (vector/make) (vector/make) "TILT"))
([position mass velocity force name]
(->star position mass velocity force name)))
(defn gravity
"牛顿引力公式"
[m1 m2 r]
(/ (* m1 m2) (* r r)))
(defn force-between
"计算两个物体之间的力"
[o1 o2]
(let [p1 (:position o1)
p2 (:position o2)
d (position/distance p1 p2)
uv (vector/unit (vector/subtract p2 p1))
g (gravity (:mass o1) (:mass o2) d)]
(vector/scale uv g)))
(defn accumulate-forces
"计算当前object的受力"
([object world]
(assoc object :force (accumulate-forces object world (vector/make))))
([object world f]
(cond
(empty? world) f
(= object (first world)) (recur object (rest world) f)
:else (recur object
(rest world)
(vector/add f (force-between object (first world)))))))
(defn calculate-forces-on-all
"计算所有object的受力"
[world]
(map #(accumulate-forces % world) world))
(defn accelerate
"计算当前object的速度"
[{f :force :as o}]
(let [m (:mass o)
v (:velocity o)
av (vector/add v (vector/scale f (/ 1.0 m)))]
(assoc o :velocity av)))
(defn accelerate-all
"计算所有object的速度"
[world]
(map accelerate world))
(defn reposition
"计算object的位置, 并生成新的object"
[{v :velocity p :position :as o}]
(assoc o :position (position/add p v)))
(defn reposition-all
"重新计算所有的object的位置"
[world]
(map reposition world))
(def audio-chan
"unbuffered channel
发生碰撞的时候, 往这个channel放一个msg"
(chan))
(defn collided?
"判断两个object是否发生碰撞, 如果发生碰撞, 往channel发一个msg"
[o1 o2]
(let [p1 (:position o1)
p2 (:position o2)
mass (+ (:mass o1) (:mass o2))
distance (position/distance p1 p2)
radius (Math/sqrt mass)]
(if (>= (max 3 radius) distance)
(do (go (>! audio-chan "1"))
true))))
(defn center-of-mass
"质量中心计算, 返回一个position, 代表两个object的质量中心"
[{p1 :position, m1 :mass}
{p2 :position, m2 :mass}]
(let [s (/ m1 (+ m1 m2))
uv (vector/unit (vector/subtract p2 p1))
d (vector/scale uv s)]
(position/add p1 d)))
(defn merge
"两个object合并的计算, 合并为一个新的object, 计算新object的位置, 质量, 速度, 受力, 名称"
[{n1 :name, m1 :mass, v1 :velocity f1 :force, :as o1}
{n2 :name, m2 :mass, v2 :velocity f2 :force, :as o2}]
(let [p (center-of-mass o1 o2)
m (+ m1 m2)
mv1 (vector/scale v1 m1)
mv2 (vector/scale v2 m2)
v (vector/scale (vector/add mv1 mv2) (/ 1 m))
f (vector/add f1 f2)
n (if (> m1 m2) (str n1 "." n2) (str n2 "." n1))]
(make p m v f n)))
(defn remove-obj
"移除object"
[o world]
(remove #(= o %) world))
(defn difference-list
"从l1中移除所有l2中的object"
[l1 l2]
(reduce #(remove-obj %2 %1) l1 l2))
(defn collide-all
"计算碰撞关系"
[world]
(loop [colliding-world world collided-world [] collisions []]
(if (empty? colliding-world)
[collisions collided-world]
(let [impactor (first colliding-world)
targets (rest colliding-world)
colliders (filter #(collided? impactor %) targets)
merger (reduce merge impactor colliders)
survivors (difference-list targets colliders)
new-collisions (if (empty? colliders) collisions (conj collisions (:position merger)))]
(recur survivors (conj collided-world merger) new-collisions)))))
(defn update-all [world]
(let [[collisions collided-world] (collide-all world)]
[collisions (-> collided-world
calculate-forces-on-all
accelerate-all
reposition-all)]))
使用reagent定义UI
(ns solar-system.core
"UI定义"
(:require-macros [cljs.core.async.macros :refer [go go-loop]])
(:require [reagent.core :as reagent :refer [atom]]
[physics.vector :as vector]
[physics.star :as star]
[physics.position :as position]
[clojure.pprint :refer [pprint]]
[cljs.core.async :refer [ <!]])
(:refer-clojure :exclude (vector)))
(def SOLAR-SYSTEM-SIZE 1000) ;; 世界大小
(def STAR-NUM 100) ;; 行星的数量
(def CENTER (position/make (/ SOLAR-SYSTEM-SIZE 2.0) (/ SOLAR-SYSTEM-SIZE 2.0)))
(def MAX-STAR-MESS 5) ;; 行星的质量上限
(def SUM-MASS 2000) ;; 太阳质量
(def ORBIT-HISTORY-LENGHT 2000) ;; 轨道历史长度
(def REDIUS-COEFFICIENT 3.0)
(defn random-velocity
"随机速度"
[p sun]
(let [sp (:position sun)
sd (position/distance p sp)
v (Math/sqrt (/ 1 sd))
direction (vector/rotate90 (vector/unit (vector/subtract p sp)))]
(vector/scale direction (+ (rand 0.01) (* v 13.5 3)))))
(defn random-position
"在太阳的周围随机设定位置"
[sun-position]
(let [r (+ (rand 300) 30)
theta (rand (* 2 Math/PI))]
(position/add sun-position (position/make (* r (Math/cos theta)) (* r (Math/sin theta))))))
(defn random-star
"增加n个随机行星"
[sun n]
(let [sp (:position sun)
p (random-position sp)]
(star/make p (rand MAX-STAR-MESS) (random-velocity p sun) (vector/make) (str "r" n))))
(defn create-world
"创世纪"
[]
(let [v0 (vector/make)
sun (star/make CENTER SUM-MASS (vector/make 0 0) v0 "sun")]
(loop [world [sun]
n STAR-NUM]
(if (zero? n)
world
(recur (conj world (random-star sun n)) (dec n))))))
(defonce app-state
(atom {:text "我的太阳系!"
:my-world (create-world)}))
(defonce orbit-history (atom []))
(defn reset-orbit []
(reset! orbit-history []))
(defn redius-for-star
"星球的直径"
[star]
(* REDIUS-COEFFICIENT (Math/pow (get-in star [:mass])
(/ 1 3))))
(defn all-stars
"所有的星球, 以圆形表示"
[]
(for [star (:my-world @app-state)]
^{:key (:name star)}
[:circle {:r (redius-for-star star)
:fill (if (clojure.string/includes? (:name star) "sun")
"red"
"black")
:cx (get-in star [:position :x])
:cy (get-in star [:position :y])}]))
(defn all-history
"所有的轨迹, 每个历史点都用一个小圆圈表示"
[]
(for [orbit-point @orbit-history]
^{:key (str (:x orbit-point) (:y orbit-point))}
[:circle {:r 1
:fill "blue"
:cx (:x orbit-point)
:cy (:y orbit-point)}]))
(defn main-svg
"所有星球的svg"
[]
[:svg {:view-box [0 0 SOLAR-SYSTEM-SIZE SOLAR-SYSTEM-SIZE]
:width 500
:height 500}
(all-stars)])
(defn path-svg
"所有轨迹的svg"
[]
[:svg {:view-box [0 0 SOLAR-SYSTEM-SIZE SOLAR-SYSTEM-SIZE]
:width 500
:height 500}
(all-history)])
(defn scene
"svg全景"
[]
[:svg {:width 500
:height 500}
[main-svg]
[path-svg]])
(defn play-audio []
(.play (.getElementById js/document "play")))
(go-loop []
(<! star/audio-chan)
(play-audio)
(recur))
(defn solar []
[:center
[:h1 (:text @app-state)]
[:audio {:id "play" :src "/blog/img/orbit/explosion.mp3" :crossOrigin "anonymous"}]
[scene]])
(reagent/render-component [solar]
js/klipse-container)
(defn run-loop []
(let [new-world (vec (second (star/update-all (:my-world @app-state))))
orbit-point (map #(get-in % [:position]) new-world)
history (into (take ORBIT-HISTORY-LENGHT @orbit-history) orbit-point)]
(swap! app-state assoc-in [:my-world] new-world)
(reset! orbit-history history)))
(defonce time-updater (js/setInterval run-loop 16))
做些调整
- 调整星球大小
- 关闭或者打开轨迹
- 三体运动模拟
- 简化模型, 忽略各个行星之间的吸引力
- 调整行星数量, 质量
- 移动太阳位置
- etc...
Hava fun
(comment (swap! app-state update-in [:my-world 0 :position :x ] #(+ % 100)))
@app-state