% ch1example4prg1.m
g=9.8;   % 重力加速度
v0=0;    % 初始速度
y0=1;    % 初始位置
m=0.1;   % 小球質量
t0=0;    % 起始時間
K=0.8;   % 彈跳的損耗系數
N=5000;  % 仿真的總步進數
dt=0.005;% 仿真步長
v=v0;    % 初狀態
y=y0;   
vx=0;
vz=0;
sx=0;
sz=0;
for k=1:N 
    if y > 0        % 小球在空中的動力方程計算
        v=v -g*dt;
        y=y +v*dt;
    else            % 碰擊瞬間的計算
        y=-K.*v*dt;  
        v=-K.*v-g*dt;   
    end 
    Fx=randn;       % x水平方向的隨機力,方差為1
    ax=Fx./m;       % Fx導致的x水平方向的加速度
    vx=vx+ax*dt;    % 小球在x水平方向的瞬時速度
    sx=sx+vx*dt;    % 小球在x水平方向的位移
    Fz=randn;       % z水平方向的隨機力,方差為1
    az=Fz./m;       % Fz導致的z水平方向的加速度
    vz=vz+az*dt;    % 小球在z水平方向的瞬時速度
    sz=sz+vz*dt;    % 小球在z水平方向的位移
    plot3(sx,sz,y,'.'); grid on;hold on;
    axis([-2 2 -2 2 0 1]);          % 坐標范圍固定
    set(gcf,'DoubleBuffer','on');   % 雙緩沖避免作圖閃爍
    xlabel('水平方向x');ylabel('水平方向z');zlabel('垂直方向y');
    drawnow;                        % 作圖
end