% Script file ball.m
%
% Purpose:
%    This program calculates the distance traveled by a ball
%  thrown at a specified angle "theta" and a specified velocity
%  "vo" form a point,ignoring air friction .It calculates the angle
%  yeilding maximum range,and also plots selected trajectories.
%
% Define variables:
%   conv      degrees to radians conv factor
%   grav      The gravity accel
%   ii,jj     Loop index
%   index     The maximum range in array
%   maxangle  The angle that gives the maximum range
%   maxrange  Maximum range 
%   range     range for a specified angle
%   time      Time
%   theta     Initial angle
%   fly_time  the total trajectory time
%   vo        The initial velocity
%   vxo       x-component of the initial velocity
%   vyo       y-component of the initial velocity
%   x         x-position of ball
%   y         y-position of ball

%  定義常數數值
conv=pi/180;
grav=-9.82;
vo=input('Enter the initial velocity:');

range=zeros(1,91);

% 計算最大的水平距離
for ii=1:91
     theta=ii-1;
     vxo=vo*cos(theta*conv);
     vyo=vo*sin(theta*conv);
     max_time=-2*vyo/grav;
     range(ii)=vxo*max_time;
end

%顯示計算水平距離的列表
fprintf('Range versus angle theta:\n');
for ii=1:5:91
    theta=ii-1;
    fprintf('%2d %8.4f\n',theta,range(ii));
end

% 計算最大的角度和水平距離
[maxrange index]=max(range);
maxangle=index-1;
fprintf('\n Max range is %8.4f at %2d degrees.\n',maxrange,maxangle);

% 繪制軌跡圖形
for ii=5:10:80
     theta=ii;
     vxo=vo*cos(theta*conv);
     vyo=vo*sin(theta*conv);
     max_time=-2*vyo/grav;
     %  計算小球軌跡的x，y坐標數值   
     x=zeros(1,21);
     y=zeros(1,21);
     for jj=1:21
         time=(jj-1)*max_time/20;
         x(jj)=vxo*time;
         y(jj)=vyo*time+0.5*grav*time^2;
     end
     plot(x,y,'g');
     if ii==5
         hold on;
     end 
end
% 添加圖形的標題和坐標軸名稱
title('\bf Trajectory of Ball vs Initial Angle\theta');
xlabel('\bf\itx \rm\bf(meters)');
ylabel('\bf\ity \rm\bf(meters)');
axis([0 max(range)+5 0 -vo^2/2/grav]);
grid on;

% 繪制最大水平的軌跡圖形
vxo=vo*cos(maxangle*conv);
vyo=vo*sin(maxangle*conv);
max_time=-2*vyo/grav;
     %  Calculate the (x,y) position   
     x=zeros(1,21);
     y=zeros(1,21);
     for jj=1:21
         time=(jj-1)*max_time/20;
         x(jj)=vxo*time;
         y(jj)=vyo*time+0.5*grav*time^2;
     end
     plot(x,y,'r','Linewidth',2);
     hold off