問題：想知道炮彈打到那裏
物理量：軌跡，即位置隨時間變化的函數 x(t), y(t)
物理定律（數學模型）：f = m a ，在此 f 純粹來自於重力 （粗體字代表是向量）
粒子在地表附近所受的重力是 g （g = GmM/R），g = 9.8 牛頓/公斤 or 米/(秒^2)
f = ma = m ( ax ex + ay ey ) ，其中 ax = dvx/dt 、 ay = dvy/dt
只有重力時，f = fy ey = g m ey ，也就是說 fy = g m
演算法：
利用 Euler 演算法：
位置 xn+1 = xn + vx,n Δt ； yn+1 = yn + vy,nΔt
速度 vx,n+1 = vx,n + ax,n Δt ； vy,n+1 = vy,n + ay,nΔt
進階補充：對計算穩定性及精確度，有更好的選擇
Euler-Richarson 演算法， 先求出
ax,n = fx(rn,vn,t)/m
vx,mid = vx,n + (1/2) ax,nΔt
xmid = xn + (1/2) vx,nΔt
ax,mid = fx( xmid, vx,mid, t+(1/2)Δt ) / m
然後才做
vx,n+1 = vx,n + ax,midΔt
xn+1 = xn + vx,midΔt
註：為什麼 Euler-Richarson 演算法 比 Euler 的好，可參考 Gould 課本上的推導： p1、p2
程式流程：
(0) 開始，宣告變數
(1) 問質量、初速度（含初速率與角度）
(2) 建立位置 x 與 y、速度 vx 與 vy 的初始值
(3) 用演算法求下一個小時段 Δt 後的位置與速度
(4) 畫出粒子位置並檢查高度是否己達地面，未達地面則回到 (3)
(5) 問要不要再射一次，若要則回到 (1)
(6) 結束
你等不及了嗎？先偷玩一下老師己經做好的可執行檔 cannon.x
你做不出來嗎？先偷看一下老師己經寫好的範圍程式 cannon.f cannon.f.txt （A simpler version for starter cannon_simple.f）
（非不得已請勿先看）
動動腦、動動手，自己再改造（或回家做）：
1. 畫個大砲把炮管的角度表現出來
2. 以火藥包裝填數表現炮彈動能，轉換成初速率
3. 物體被炸破片的彈出（配合亂數產生器 ran3.f，使用方法）
4. 試寫一個煙火程式，會多段開花變色的
分析與討論：
想預測炮彈打到那裏，主要是基於那一條物理定律？
什麼樣的先決（初始）條件決定了炮彈的軌跡（和落點）？
同樣的火藥包充填數，怎麼打炮彈才飛得最遠？
問題：在有空氣阻力的情況下，炮彈軌跡變得如何？
物理量：一樣是軌跡，即位置隨隨時間的變化 x(t), y(t)
物理定律（數學模型）：f = m a，受力 f 同時包含了重力及空氣阻力
空氣阻力非常複雜（因此飛機或高速火車的阻力設計需要靠風洞實驗或電腦模擬來進行），基本上它是與速率有關，速率越高阻力越大，靜止的物體則沒有空氣阻力。可見它是一個與速率有關的量，常見的簡化公式，有隨速率一次方呈正比的，也有隨兩次方呈正比的，方向是速度反向（負號）的
Fd = - k1 v
Fd' = - k2 v2
由於有 x 與 y 兩個分量，假設速度與水平線的夾角為 θ，即 θ = tan-1(vy/vx) [ Fortran 程式寫成 theta = atan(vy/vx) ] ，
Fd,x = Fd cosθ
Fd,y = Fd sinθ
（或者是，直接用 vx/√(vx2+vy2) 及 vy/√(vx2+vy2) 會比先求 θ = atan(vx/vy) 取角度才去求 cos(θ) 及 sin(θ) 更好）
因此，在有關力之部分的公式，
fx = Fd cosθ
fy = - m g + Fd sinθ
其他的部分，與上一個議題相同
新議題：如何求得室氣阻力係數 k？
利用終端速度，即垂直落下過程中重力與阻力達平衡，使合力為零不再加速，而呈等速運動時的速度。
合力 = 0 = - m g + Fd = - m g + k1 vy2
即 k1 = (m g) / vy2，此終端速度可藉由實驗量得
一個粉筆頭大小的石子其終端速度大約是每秒三十公尺。另外，一個重 0.254 克 、半徑 2.54 公分的保利龍球有以下的實測數據（引用 Physics Teacher 24, 153 (1986) )
| 時刻（秒） |
位置（公尺） |
| -0.132 |
0.0 |
| 0.0 |
0.075 |
| 0.1 |
0.260 |
| 0.2 |
0.525 |
| 0.3 |
0.870 |
| 0.4 |
1.27 |
| 0.5 |
1.73 |
| 0.6 |
2.23 |
| 0.7 |
2.77 |
| 0.8 |
3.35 |
想想看，我們要怎麼從上列的數據求得空氣阻力係數 k？
老師寫的範例程式：styrofoam_k.f.txt styrofoam_k.f styrofoam_k.x
上面表格數據的資料檔：styrofoam.txt
（儘量自己想，不會再參考）
演算法：
同樣採用 Euler 或 Euler-Richardson 演算法
程式設計：
與上一個議題相同
參考範例程式
cannon_drift.f cannon_drift.f.txt cannon_drift.x
簡化一些不必要功能，邊想邊做三部曲：
fly_free.f
只有自由落體，抵達地面停止並列印高度供檢查，沒有繪圖功能。
fly_free_plot.f
加上 繪圖功能
fly.f
最後，再加上空氣阻力
質量與空氣阻力係數對拋射體軌跡的影響
mass = 1.0, k = 2.0
mass = 0.5, k = 1.0
mass = 2.0, k = 2.0
分析與討論：
以同樣的力道向上拋一個小石子，在有空氣阻力與無室空氣阻力的兩種情況比較下，那一種會在空中停留較久？
那一種軌跡像羽毛球？那一種像鉛球？
在有空氣阻力的情況下，拋射物的角度要如何設定其距離才會最遠？
挑戰：
試著做出你自己的憤怒鳥遊戲 (技術上的細節歡迎與教授討論 ) 。
挑戰：
試作沒有手操控的「打磚塊」遊戲。
彈珠與光滑釘子的交互作用模式
N = Nx ex + Ny ey
| N | = 1
V dot N = Vx * Nx + Vy * Ny
碰撞後 V 平行於 N 上的分量因反彈而方向相反，與 N 垂直的分量則繼續維持，因此碰撞後
V┴ = (-1)* (V dot N) N
V// = V - (V dot N) N = V + V┴
新的速度則為 V// + V┴
老師提供的的範例程式 pinball.f pinball.f.txt pinball.x ，輸入 範例
以 ran3.f 內含的 ran3(iseed) 函式為例：
iseed = 7
x = ran3(iseed)
y = ran3(iseed)
準備任意一個整數值放在某個整數變數 iseed 的值內，依上述的方法每次使用 ran3( ) ，則 x 與 y 都會獲得一個 0.0 到 1.0 之間的數值，每次的數值雖然不一樣，但這些數值在 0 到 1 之間出現的機率是均等的。另外，直接把整數填入 ran3( ) 卻並不支援，如 xx = ran3(7)，會導致 "區段錯誤 (Sgementation Fault)" 。
試試看，寫一個丟骰子程式。例範：dice.f ，記得與 ran3.f 一起編譯，如 gfortran diec.f ran3.f <Enter>。
另外 gfortran 支援內建副程式 random_number(實數變數)，就不必另外用外部的亂數函式了，可參考另一骰子程式範例 dice_gfortran.f。重設亂數種子的副程式是 random_seed 細節請見： http://www.nsc.liu.se/~boein/f77to90/a5.html#section21c
