前言
近日拜讀了Daniel Shiffman先生的著作——《代碼本色：用編程模擬自然系統(tǒng)》，決定做一組習(xí)作间狂，來(lái)對(duì)書(shū)中提到的隨機(jī)行為及牛頓運(yùn)動(dòng)學(xué)進(jìn)行理解暂吉，并對(duì)一些實(shí)例進(jìn)行拓展學(xué)習(xí)陆蟆，從而提升自己相關(guān)方面的知識(shí)水平和實(shí)踐能力。
《代碼本色》第2章概述
在第2章中凤藏，作者向我們介紹了牛頓的力學(xué)三大定律奸忽。分別為：
牛頓第一運(yùn)動(dòng)定律：除非有不均衡外力的作用，否則物體始終保持靜止或勻速直線運(yùn)動(dòng)狀態(tài)清笨。
牛頓第二運(yùn)動(dòng)定律：力等于質(zhì)量乘以加速度月杉。
牛頓第三運(yùn)動(dòng)定律：每個(gè)作用力都有一個(gè)大小相等、方向相反的反作用力抠艾。
隨后苛萎，作者通過(guò)實(shí)例告訴我們?nèi)绾卧赑rocessing里模擬力，這里包括了制造質(zhì)量检号，制造各種外力腌歉，如引力，摩擦力齐苛，空氣阻力和流體阻力等等翘盖。
在這一章中，我印象比較深的點(diǎn)是各種力的代碼模擬凹蜂，通過(guò)面向?qū)ο蟮乃枷脞裳保瑢⒘Φ谋憩F(xiàn)用代碼實(shí)現(xiàn)阁危。下面，我將介紹我根據(jù)本章內(nèi)容而創(chuàng)作的Processing編程習(xí)作汰瘫。
習(xí)作
創(chuàng)作過(guò)程
對(duì)于浮力的模擬狂打，其實(shí)說(shuō)起來(lái)容易做起來(lái)難。首先混弥，我利用了書(shū)中的一個(gè)例子趴乡，即2-5，該例子中已經(jīng)實(shí)現(xiàn)了重力和流體阻力蝗拿。書(shū)中定義了一個(gè)新的類(lèi)型——Liquid晾捏，并在物體進(jìn)入液體時(shí)給他一個(gè)drag的力，此外哀托，阻力的大小等于阻力系數(shù)乘以對(duì)象速度的平方惦辛，且方向與Mover對(duì)象的速度方向相反。這樣一來(lái)就定義了這個(gè)阻力萤捆，實(shí)現(xiàn)了流體阻力的效果裙品。
首先俗批，我們來(lái)看一下沒(méi)有浮力的情況：
而我要做的是在這個(gè)基礎(chǔ)上在Liquid類(lèi)中添加浮力的效果俗或，首先，我們來(lái)復(fù)習(xí)一下初中物理的知識(shí)——浮力計(jì)算公式岁忘。
公式1. F感廖俊＝F向上－F向下
“F向上”指下表面受到的向上的力，F(xiàn)向下指上表面受到的向下的力干像，這是浮力的最原始的計(jì)算公式帅腌。
公式2. F浮＝G排＝ρ液gV排
這是根據(jù)阿基米德原理得到的麻汰，V排指排出液體的體積速客，ρ液指液體密度。
公式3. F肝弼辍＝G物
即ρ液gV物溺职，利用二力平衡，即根據(jù)漂浮位喂、懸浮的物體浮力與自重相等浪耘。
公式4. F浮＝G物－F拉
測(cè)量浮力時(shí)根據(jù)此公式計(jì)算塑崖，F(xiàn)拉指的是彈簧測(cè)力計(jì)的拉力七冲。
公式5. F浮=ρgh
h指的是物體全部浸入液體中時(shí)表面距離液面的高度。
我們比較常用的是公式2规婆。浮力是什么時(shí)候開(kāi)始起作用的呢澜躺？答案是當(dāng)物體接觸到液體表面的那一刻開(kāi)始蝉稳。根據(jù)公式2，物體浸入液體的面積（其實(shí)應(yīng)該是體積掘鄙，在二維中稱(chēng)之為面積）就是我們說(shuō)的V排颠区。
對(duì)此，我畫(huà)了個(gè)草圖幫助理解：
那么通铲，我們?cè)诔绦蛑腥菀椎玫降闹凳乔蛐牡淖鴺?biāo)和液體表面的高度（即y坐標(biāo)）毕莱，那么我們?cè)撊绾斡?jì)算S的面積呢？
我分了兩種情況颅夺，分別畫(huà)了圖幫助理解：
情況1：
這種情況是圓心坐標(biāo)已經(jīng)低于水平面了朋截，換言之，就是一半以上的圓已經(jīng)進(jìn)入了液體吧黄。我們假設(shè)：
R為圓的半徑
h為水面高度-圓心的高度
那么S這一塊的面積就是紅色的扇形面積加上黃色的這個(gè)等腰三角形的面積部服。
其中，扇形的面積又可以通過(guò)圖中標(biāo)出的θ角來(lái)求出拗慨，公式為（1-2θ/2π）×S圓廓八。
三角形的面積可以根據(jù)勾股定理算出d，面積為dh2/2赵抢。
情況2
接下來(lái)剧蹂，我們來(lái)看第二種情況。這種情況是圓心坐標(biāo)高于水平面了烦却，換言之宠叼，進(jìn)入液體的面積不到一半。
需要注意的是其爵，這里的h應(yīng)該是-h冒冬，原因是圓心還沒(méi)有過(guò)液體表面，所以這一段距離為-h摩渺。其他的計(jì)算方式和情況1相似简烤。唯一不同的是，此時(shí)的面積S為扇形面積-三角形面積：
這時(shí)摇幻，我在進(jìn)行公式推導(dǎo)時(shí)發(fā)現(xiàn)了一個(gè)問(wèn)題横侦，根據(jù)arccos()函數(shù)的性質(zhì)，arccos(-x)=2π-arccos(x)囚企。運(yùn)用在這個(gè)問(wèn)題中丈咐，就是arccos(h/R)=2π-arccos(-h/R)。而這個(gè)2π又剛好讓扇形面積的正反兩面抵消龙宏。經(jīng)過(guò)化簡(jiǎn)棵逊，我們可以將兩個(gè)公式合并成同一個(gè)。公式如下：
至此银酗，V排的公式就已經(jīng)得出辆影，接下來(lái)就可以編寫(xiě)代碼了：
首先徒像，在Liquid類(lèi)中，我添加了一個(gè)獲取浮力的方法：
PVector buoyancy (Mover m)
{//浮力
PVector l = m. position;
float h=y-l. y;//水面高度和球心高度的差
float R=m. mass*16/2;//半徑
float Sf= (R*R) * (2*PI-acos (h/R)) +sqrt (R*R-h*h) *h;//V排
PVector buoyancyForce;
buoyancyForce=new PVector (0, -0.1*Sf, 0);
buoyancyForce.mult(0.001);//由于整個(gè)程序中的力都比較小蛙讥，浮力也不能太大
return buoyancyForce;
}
接著锯蛀，在主程序中，將這個(gè)力應(yīng)用次慢。
if(liquid.y-movers[i].position.y<movers[i].mass*8)//如果圓接觸液體表面
{
PVector buoyancy=liquid.buoyancy(movers[i]);
movers[i].applyForce(buoyancy);
}
運(yùn)行結(jié)果如下：
為什么球突然消失了呢旁涤，這使我非常困惑。我檢查了很久迫像，也沒(méi)有發(fā)現(xiàn)問(wèn)題在哪里劈愚。最后我利用processing的調(diào)試器，進(jìn)行逐步調(diào)試并仔細(xì)觀察變量變化闻妓，終于發(fā)現(xiàn)了問(wèn)題所在：
可以看到菌羽，消失的小球的y變成了NaN。
我去查閱了資料由缆，發(fā)現(xiàn)會(huì)出現(xiàn)NaN的最常見(jiàn)情況就是用把0作為了除數(shù)注祖，或者是其他的計(jì)算錯(cuò)誤。我開(kāi)始認(rèn)真分析我的代碼中哪里會(huì)出現(xiàn)這種錯(cuò)誤均唉。最后是晨，鎖定了一個(gè)嫌疑人——acos（）函數(shù)，也就是arccos函數(shù)浸卦。
我發(fā)現(xiàn)在h進(jìn)行變化的時(shí)候署鸡，會(huì)無(wú)限接近于R案糙，也就是球體逐漸整個(gè)進(jìn)入液體的過(guò)程限嫌。而arccos的函數(shù)圖像是這樣的：
也就是說(shuō)，arccos(1)和arccos(-1)的結(jié)果是無(wú)窮时捌。而h/R的值是會(huì)到達(dá)這個(gè)取值范圍的怒医！
所以我做了個(gè)判斷，如果h的絕對(duì)值逼近R了奢讨，這種情況也就是整個(gè)球就快要全部浸入了稚叹，這是，我更換浮力的計(jì)算公式拿诸。代碼如下：
PVector buoyancy (Mover m)
{//浮力
PVector l = m. position;
float h=y-l. y;//水面高度和球心高度的差
float R=m. mass*16/2;//半徑
float Sf= (R*R/2) * (2*PI-acos (h/R)) +sqrt (R*R-h*h) *h;//V排
PVector buoyancyForce;
if(-h-R<0)
{
buoyancyForce=new PVector (0, -0.1*Sf, 0);
buoyancyForce.mult(0.001);//由于整個(gè)程序中的力都比較小扒袖，浮力也不能太大
println(buoyancyForce);
}
else
{
buoyancyForce=new PVector(0,-density *m.mass,0);
}
return buoyancyForce;
}
其中density為液體的密度。
這樣一來(lái)亩码，就可以保證不會(huì)觸碰到這個(gè)“危險(xiǎn)地帶”季率。最終效果如下：
最后，我又根據(jù)第0章的知識(shí)描沟，使用了Perlin噪聲飒泻，模擬了海面的效果：
Perlin噪聲的這段代碼如下：
void display() {
stroke(#09AAE3);
fill(50);
//rect(x, y, w, h);
float x = 0;
while (x < width) {
line(x, 180 + 50 * noise(x / 100, t), x, height);
x = x + 1;
}
t = t + 0.01;
}
最后鞭光，貼出全部代碼：
Mover類(lèi)：
class Mover {
// position, velocity, and acceleration
PVector position;
PVector velocity;
PVector acceleration;
// Mass is tied to size
float mass;
Mover(float m, float x, float y) {
mass = m;
position = new PVector(x, y);
velocity = new PVector(0, 0);
acceleration = new PVector(0, 0);
}
// Newton's 2nd law: F = M * A
// or A = F / M
void applyForce(PVector force) {
// Divide by mass
PVector f = PVector.div(force, mass);
// Accumulate all forces in acceleration
acceleration.add(f);
}
void update() {
// Velocity changes according to acceleration
velocity.add(acceleration);
// position changes by velocity
position.add(velocity);
// We must clear acceleration each frame
acceleration.mult(0);
}
// Draw Mover
void display() {
stroke(0);
strokeWeight(2);
fill(255);
ellipse(position.x, position.y, mass*16, mass*16);
}
// Bounce off bottom of window
void checkEdges() {
if (position.y > height) {
velocity.y *= -0.9; // A little dampening when hitting the bottom
position.y = height;
}
}
}
Liquid類(lèi)：
// Liquid class
class Liquid {
float t=0;
// Liquid is a rectangle
float x, y, w, h;
// Coefficient of drag
float c;
float density;
Liquid(float x_, float y_, float w_, float h_, float c_,float density_) {
x = x_;
y = y_;
w = w_;
h = h_;
c = c_;
density=density_;
}
// Is the Mover in the Liquid?
boolean contains(Mover m) {
PVector l = m.position;
return l.x > x && l.x < x + w && l.y > y && l.y < y + h;
}
// Calculate drag force
PVector drag(Mover m) {
// Magnitude is coefficient * speed squared
float speed = m.velocity.mag();
float dragMagnitude = c * speed * speed;
// Direction is inverse of velocity
PVector dragForce = m.velocity.get();
dragForce.mult(-1);
// Scale according to magnitude
// dragForce.setMag(dragMagnitude);
dragForce.normalize();
dragForce.mult(dragMagnitude);
return dragForce;
}
PVector buoyancy (Mover m)
{//浮力
PVector l = m. position;
float h=y-l. y;//水面高度和球心高度的差
float R=m. mass*16/2;//半徑
float Sf= (R*R/2) * (2*PI-acos (h/R)) +sqrt (R*R-h*h) *h;//V排
PVector buoyancyForce;
if(-h-R<0)
{
buoyancyForce=new PVector (0, -0.1*Sf, 0);
buoyancyForce.mult(0.001);//由于整個(gè)程序中的力都比較小，浮力也不能太大
println(buoyancyForce);
}
else
{
buoyancyForce=new PVector(0,-density *m.mass,0);
}
return buoyancyForce;
}
void display() {
stroke(#09AAE3);
fill(50);
//rect(x, y, w, h);
float x = 0;
while (x < width) {
line(x, 180 + 50 * noise(x / 100, t), x, height);
x = x + 1;
}
t = t + 0.01;
}
}
主程序類(lèi)：
Mover[] movers = new Mover[9];
// Liquid
Liquid liquid;
void setup() {
size(640, 360);
reset();
// Create liquid object
liquid = new Liquid(0, height/2, width, height/2, 0.1,0.12);
}
void draw() {
background(255);
// Draw water
liquid.display();
for (int i = 0; i < movers.length; i++) {
// Is the Mover in the liquid?
if (liquid.contains(movers[i])) {
// Calculate drag force
PVector dragForce = liquid.drag(movers[i]);
// Apply drag force to Mover
movers[i].applyForce(dragForce);
}
// Gravity is scaled by mass here!
PVector gravity = new PVector(0, 0.1*movers[i].mass);
// Apply gravity
movers[i].applyForce(gravity);
if(liquid.y-movers[i].position.y<movers[i].mass*8)//如果圓接觸液體表面
{
PVector buoyancy=liquid.buoyancy(movers[i]);
movers[i].applyForce(buoyancy);
}
// Update and display
movers[i].update();
movers[i].display();
movers[i].checkEdges();
}
fill(0);
text("click mouse to reset", 10, 30);
}
void mousePressed() {
reset();
}
// Restart all the Mover objects randomly
void reset() {
for (int i = 0; i < movers.length; i++) {
movers[i] = new Mover(random(2, 4), 40+i*70, 0);
}
}
總結(jié)
本文介紹了《代碼本色：用編程模擬自然系統(tǒng)》第2章的主要內(nèi)容泞遗，并在示例的基礎(chǔ)上進(jìn)行了拓展性的創(chuàng)作惰许，創(chuàng)作了流體阻力重力和浮力的模型。
