%Mason Averill
%ME-544 Modeling and Simulation of ME Systems: Fall 2020
%Special Project 1: Initial Velocity Required for 600ft Home Run W/
%Drag/Air Resistance
%08_29_2020


%Environmental and Physical Parameters

mass=0.145; %mass of baseball in kg
diameter=0.074;%diameter of baseball in m
air_density=1.225; %density of ambient air in kg/m^3
gravity=9.81; %acceleration due to gravity in m/s^2

%Other Input Parameters
angular_velocity=160;%angular velocity of baseball in rad/s (at 0 no drag will be considered, as function for determining drag coefficients is no longer valid)
wind_speed=2; %wind speed in m/s, + for same direction as baseball travel, - for opposite baseball travel
player_height=1.9; %height of player in m
hit_angle=40; %initial launch angle w.r.t the ground in degrees
hit_distance=600;%desired hit distance in feet

%convert hit_angle to radians for future use

hit_angle=hit_angle*pi()/180;

area=pi()*diameter^2/4;%effective area of interest for drag force calculations



%Initialize lower bound for initial velocity

initial_velocity_ground=30; %initial velocity (relative to the ground) guess in m/s (your CPU will appreciate a close guess at high resolutions)


%%LOOP

sx=0; %Initial position in x
delta_velocity=0.05;

while(sx<hit_distance/3.281)
    
    
    initial_velocity_ground=initial_velocity_ground+delta_velocity;

%Initialize initial conditions

sx=0; %Initial position in x
sy=player_height-.3; %initial position in y

vx=initial_velocity_ground*cos(hit_angle);%Initial velocity in x
vy=initial_velocity_ground*sin(hit_angle);%Initial velocity in y

ax=0;%initial acceleration in x
ay=-1*gravity;%initial acceleration in y


%Set time

t=0;

%Set increments in time desired

delta_t=.05; %in seconds

%Calculate trajectory based on given inputs and conditions

weight_force=-mass*gravity;

spin_ratio_x=0;%create variable spin_ratio_x and assign a value of 0
spin_ratio_y=0;%create variable spin_ratio_y and assign a value of 0
fd_x=0;%create variable Fd_x and assign a value of 0 (drag force)
fd_y=0;%create variable fd_y and assign a value of 0 (drag force)
fl_y=0;%create variable fl_y and assign a value of 0 (lift force)
cd_x=0;%create variable cd_x and assign a value of 0 (drag coefficient)
cd_y=0;%create variable cd_y and assign a value of 0 (drag coefficient)
cl_y=0;%create variable cl_y and assign a value of 0 (drag coefficient)

sx_store=[];
sy_store=[];

vx_store=[];
vy_store=[];

ax_store=[];
ay_store=[];



sx_store=[sx];
sy_store=[sy];

vx_store=[vx];
vy_store=[vy];

ax_store=[ax];
ay_store=[ay];

i=2;

while(sy>0)

    
%Determine Forces acting on projectile in x and y directions

%Drag Forces in x

    %Calculate Spin Ratio

    spin_ratio_x=(angular_velocity*diameter)/(2*(vx-wind_speed));

            %Determine Coefficient of Drag
            
            if(spin_ratio_x<5 && spin_ratio_x>0.05)
            
            cd_x=-0.0001474916*(spin_ratio_x)^10+0.0037232759*(spin_ratio_x)^9+-0.0400229087*(spin_ratio_x)^8+0.2391758192*(spin_ratio_x)^7+-0.8712539421*(spin_ratio_x)^6+1.9971200562*(spin_ratio_x)^5+-2.8586264480*(spin_ratio_x)^4+2.3770168518*(spin_ratio_x)^3+-0.8364561899*(spin_ratio_x)^2+-0.0998227742*spin_ratio_x+0.6286980000;
            
                    %Determine Drag force in x
            end
                    fd_x=-0.5*air_density*(vx-wind_speed)^2*area*cd_x;
%Drag Forces in y

    %Calculate Spin Ratio in y
    spin_ratio_y=abs((angular_velocity*diameter)/(2*vy));
  
            %Determine Coefficient of Drag
            if(spin_ratio_y<5 && spin_ratio_y>0.05)
            
            cd_y=-0.0001474916*(spin_ratio_y)^10+0.0037232759*(spin_ratio_y)^9+-0.0400229087*(spin_ratio_y)^8+0.2391758192*(spin_ratio_y)^7+-0.8712539421*(spin_ratio_y)^6+1.9971200562*(spin_ratio_y)^5+-2.8586264480*(spin_ratio_y)^4+2.3770168518*(spin_ratio_y)^3+-0.8364561899*(spin_ratio_y)^2+-0.0998227742*spin_ratio_y+0.6286980000;
            end 
                %Determine Drag force in y
                
           
                
                fd_y=-0.5*air_density*(vy)^2*area*cd_y;
                
                if(vy<0)
                
                    fd_y=-1*fd_y;
                    
                end

%Lift Forces in y

  %Calculate lift coefficient in y
  if(spin_ratio_x<5 &&spin_ratio_x>0.05)
  cl_y=-0.000319990769*(spin_ratio_x)^10+0.008584023046*(spin_ratio_x)^9+-0.098681399148*(spin_ratio_x)^8+0.633281021323*(spin_ratio_x)^7+-2.472378037439*(spin_ratio_x)^6+5.979932314805*(spin_ratio_x)^5+-8.665063641190*(spin_ratio_x)^4+6.735554629484*(spin_ratio_x)^3+-2.071120367*(spin_ratio_x)^2+0.215694805*spin_ratio_x+0.020034532;
  
  end

  
  %Calculate lift force in y
  
  fl_y=0.5*air_density*(vx-wind_speed)^2*area*cl_y;
  
  %All forces on ball are now known for its current conditions
  
  %Now determine position, velocity, and acceleration in x and y for this
  %increment in time
  
  ax=fd_x/mass;
  ay=(fl_y+fd_y+weight_force)/mass;
  
  %ay=(fd_y+weight_force)/mass;
  
  %ay=(weight_force)/mass;
  vx=vx+ax*delta_t;
  vy=vy+ay*delta_t;
  
  sx=sx+vx*delta_t+0.5*ax*(delta_t)^2;
  sy=sy+vy*delta_t+0.5*ay*(delta_t)^2;
    
 %Position, Velocity, and Acceleration have now been updated for this
 %increment in time
    
    
 %Store the new position, velocity, and acceleration for future plotting
 
sx_store(i)=sx;
sy_store(i)=sy;

vx_store(i)=vx;
vy_store(i)=vy;

ax_store(i)=ax;
ay_store(i)=ay;

i=i+1;

end 


end 

hit_angle=hit_angle*180/pi();
hit_angle=num2str(hit_angle);
hit_angle_string=strcat('Hit Angle=',hit_angle,' degrees');

gravity=num2str(gravity);
gravity_string=strcat('Gravitational Acceleration in Y=', gravity,' m/s^2');

angular_velocity=num2str(angular_velocity);
angular_velocity_string=strcat('Angular Velocity=',angular_velocity, ' rad/s');

wind_speed=wind_speed*2.237;
wind_speed=num2str(wind_speed);
wind_speed_string=strcat('Wind Speed=',wind_speed, ' mph');

initial_velocity_ground=initial_velocity_ground*2.237;
initial_velocity_print=num2str(initial_velocity_ground);
initial_velocity_print=strcat('The required minimum initial velocity= ', initial_velocity_print,' mph');

delta_t_for_print=num2str(delta_t);
delta_t_for_print=strcat('Delta time increments=',delta_t_for_print, ' seconds');

delta_velocity=delta_velocity*2.237;
delta_velocity=num2str(delta_velocity);
delta_velocity=strcat('Delta velocity increments= ',delta_velocity,' mph');

string_for_print=sprintf('%s\n%s\n%s\n%s\n%s\n%s\n%s\n',hit_angle_string,gravity_string,angular_velocity_string,wind_speed_string,initial_velocity_print,delta_t_for_print,delta_velocity);
    

sx_store=sx_store*3.281;
sy_store=sy_store*3.281;

figure 
plot(sx_store, sy_store)
title('Position in Y vs Position in X')
xlabel('Position in X (ft)')
ylabel('Position in Y (ft)')
xlim([0 hit_distance])
ylim([0 400])
annotation('textbox',[0.3 .5 0.4 0.3],'String',string_for_print,'FitBoxToText','on');

fprintf('The minimum required initial velocity is %.2f mph for the following input parameters:\n%s\n',initial_velocity_ground,string_for_print)