The following three m-files are a little more elegant than the ones in the "My Big Saab Jump" article. As provided, they will run on Octave version 4.2.1 and FreeMat version 4.0. They can be easily modified for other vehicles or other projectiles. To use these files, copy all of the text between dashed lines into your favorite text editor and save each to a directory of your choice using the name provided in the first dashed line. The first file, the main script, can be saved under any name of your choosing as long as it retains the m extension. However, if you change the name of either of the other two files, the first script must be changed to reference the new name. File #1, main script ----------------------------------------------- MySaab96Jump.m ----------------------------------------------- %% MySaab96Jump.m, Ver. 5, 3/31/2018 %% %% Copyright 2018, George (Ted) Yurkon %% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To %% view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative Commons, %% PO Box 1866, Mountain View, CA 94042, USA. %% %% In short, this work may be shared provided terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 %% International License are adhered to. %% % Based on projectile theory by Colin Ratcliffe at... % http://ratcliffe.site/USNA/EM375/labs/07Project/ProjectileTheory.pdf % Assumed units: meters, seconds, Newtons, kg, radians % This program runs successfully on FreeMat v4.0 and Octave 4.2.1 % It may run on Matlab as well, but will need modifications % for Scilab. Scilab provides a Matlab to Scilab converter which may help. %% % Some enhancements to the plot figure are available but Octave and % FreeMat are not compatibe. Seach for "---FreeMat---" or "---Octave---" % to find these options and follow the instructions therein. %% % Modified 2/28/2018 to include aerodynamic lift % Ver. 5 changes... % Modified 3/31/2018 to be compatible with both Octave and Freemat, and % ...to utilize the odeset Maxstep and Event options. %% clear;clc global Cd Cl0 Clslope rho A S m Pa Pr g format long CaseName=''; while isempty(CaseName) CaseName = input('Enter name:','s'); end % V0mph=0; while (V0mph <= 0) V0str=''; while isempty(V0str) V0str = input('Enter V0 in mph:','s'); % Workaround for bug in Freemat V4 end V0mph=eval(V0str); % entry speed mph end V0=V0mph*0.44704; % entry speed mph m/s % theta0deg=0; while (theta0deg < 5 || theta0deg > 90) theta0str=''; while isempty(theta0str) theta0str = input('Enter launch degrees (5 to 90):','s'); % Workaround for bug in Freemat V4 end theta0deg=eval(theta0str); % initial launch angle in degrees end theta0=theta0deg*pi/180; % radians % PaDeg=-99; while (PaDeg < -15 || PaDeg > 15) PaStr=''; while isempty(PaStr) PaStr = input('Enter attack degrees (-15 to 15):','s'); % Workaround for bug in Freemat V4 end PaDeg=eval(PaStr); % Initial attack angle degrees (different than launch angle) end Pa=PaDeg*pi/180; % PrDeg=-99; while (PrDeg < -15 || PrDeg > 15) PrStr=''; while isempty(PrStr) PrStr = input('Enter pitch rate degrees/sec (-15 to 15):','s'); % Workaround for bug in Freemat V4 end PrDeg=eval(PrStr); % attack angle rate degrees/sec end Pr=PrDeg*pi/180; % Y0=2; % Initial height (m) m=880; % mass of Saab 96, kg (car+passenger 1775lb+165lb) A=1.795+0.039; % front silhouette area plus extra for hanging wheels, m^2 S=4.146; % surface area of floor, m^2 (approx 44.625 sq ft) Cd=0.32; % most commonly cited figure Cl0=0.3; % coefficient of lift at attack angle=0 Clslope=(1.6-Cl0)/(15*pi/180); % coefficient of lift increase per radian (1.6 at 15 degrees) g=9.81; % m/s^2 rho=1.2041; % density of air, kg/m^3 %% calculate launch velocity based on ramp energy loss of m*g*h E0=0.5*m*V0^2; % entry energy Eloss1=m*g*Y0; % ramp height energy loss %% Further reduce launch energy for drag force acting through length of slope %% This is a slight overestimate of energy loss because it assumes constant speed Lslope=Y0/sin(theta0); % Length of slope assuming ramp angle is launch angle Eloss2=(0.5*Cd*A*rho*V0^2)*Lslope; % energy loss due to drag while on slope EL=E0-Eloss1-Eloss2; % adjusted launch energy VL=sqrt(2*EL/m); % launch velocity mps based on launch energy VLmph=VL/0.44704; % launch velocity mph %% perform projectile calcs tmax=5.0; % do calculations for 5 seconds of flight time tspan=[0 tmax]; %% % initial conditions as [x0, vx0, y0, vy0] IC = [0; VL*cos(theta0); Y0; VL*sin(theta0)]; options = odeset('MaxStep',0.01,'Events',@MySaab96JumpEvent); [t,oput,tf,tfVALS] = ode45(@MySaab96JumpODE, tspan, IC,options); %% x= oput(:,1)*3.281; % extract x-position from 1st column vx= oput(:,2)/0.44704; % extract x-velocity from 2nd column y= oput(:,3)*3.281; % extract y-position from 3rd column vy= oput(:,4)/0.44704; % extract y-velocity from 4th column %% %%%% Trajectory Calculations Done - Time for Data and Plot Output %%%% %% % tf is time at y-->0 xf= tfVALS(1)*3.281; % final x-position at tf vxf= tfVALS(2)/0.44704; % final x-velocity at tf yf= tfVALS(3)*3.281; % final y-position at tf vyf= tfVALS(4)/0.44704; % final y-velocity at tf %% Paf=(Pa+Pr*tf)*180/pi % landing attack angle deg vf = (vxf^2+vyf^2)^0.5 % final velocity Ymax = max(y) % height of jump at nearest time point thetaf=atan(vyf/vxf)*180/pi % landing trajectory angle deg Dy=(sin(30*pi/180)*4+0.833+0.5)*0.3048 % Occupant/driver vertical stopping distance Gocc=0.5*(vyf*0.44704)^2/Dy/g % (Vyf mph=>m/s) Dy=(sin(30*pi/180)*4+0.833)*0.3048 % Car vertical stopping distance Gcar=0.5*(vyf*0.44704)^2/Dy/g % (Vyf mph=>m/s) Dy=(sin(theta0)*4+0.833+0.5)*0.3048 % Launch perp-to-ramp stopping distance Vyl=V0*sin(theta0); % velocity perpendicular to ramp m/s Glaunch=0.5*(Vyl)^2/Dy/g % Launch g on driver Vylmph=Vyl/0.44704 % velocity perpendicular to ramp mph %% figure(1);clf; % default figure options for ---FreeMat--- and ---Octave--- %% To use the following option for ---Octave---, remove the %% and insert %% in the above figure statement... %% figure(1,'Position',[2,128,1354,320]);clf; % ---Octave--- optional plot resolution [startx,starty,width,height] %% To use the following option for ---FreeMat---, simply remove the %%... %% sizefig(1354,320) % ---FreeMat--- optional plot resolution width and height str=sprintf('%s, V0%6.2f mph, VL%6.2f mph, Angle %4.1f, Attack %4.1f at%4.1f/s, Y0%5.2f ft',CaseName,V0mph,VLmph,theta0deg,PaDeg,PrDeg,Y0*3.281); if ((PaDeg == 0) && (PrDeg == 0)) plot(x,y,'b-'); % blue plot if no angle of attack elseif (PrDeg == 0) plot(x,y,'m-'); % magenta plot if constant angle of attack else plot(x,y,'r-'); % red plot if variable angle of attack end % NOTE: The following axis statement is optional and is of the form axis([xmin xmax ymin ymax zmin zmax]) % NOTE: Use it only if you want all plots to be of exactly the same scale and axis limits axis([0 400 0 60 0 1]); % The following axis equal statement makes sure x and y are spaced equally to give an undistorted plot axis equal; grid on; title(str) xlabel('Span (ft)') ylabel('Ht. (ft)') PlotPng=sprintf('%s.png',CaseName); print(PlotPng) FileName=sprintf('%s.txt',CaseName); str=sprintf('V0:\t\t%6.2f\nVL:\t\t%6.2f\nTheta0:\t\t%6.2f\nPaDeg:\t\t%6.2f\nPrDeg:\t\t%6.2f\ntf:\t\t%6.3f\nxf:\t\t%6.2f\nYmax:\t\t%6.2f\nvxf:\t\t%6.2f\nvyf:\t\t%6.2f\nGlaunch:\t\t%6.2f\nGocc:\t\t%6.2f\nGcar:\t\t%6.2f\nvf:\t\t%6.2f\nPaFinal:\t\t%6.2f\nThetaFinal:\t\t%6.2f\nY0:\t\t%6.2f\n',V0mph,VLmph,theta0deg,PaDeg,PrDeg,tf,xf,Ymax,vxf,vyf,Glaunch,Gocc,Gcar,vf,Paf,thetaf,Y0*3.281); fp = fopen(FileName,'w'); fprintf(fp,'%s',str); fclose(fp); -------------------------------------------------------------------------------------------------------------- File #2, function called by ode45 to update first and second derivatives (velocity and acceleration) ----------------------------------------------- MySaab96JumpODE.m -------------------------------------------- function [ Fout ] = MySaab96JumpODE( t, p ) %% For use with MySaab96Jump.m, Ver. 5, 3/31/2018 %% %% Copyright 2018, George (Ted) Yurkon %% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To %% view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative Commons, %% PO Box 1866, Mountain View, CA 94042, USA. %% %% In short, this work may be shared provided terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 %% International License are adhered to. %% % Simultaneous second order differentials for projectile % motion with air resistance and aerodynamic lift % modified 2/28/2018 to include aerodynamic lift % output vector Fout has the four differential outputs x', x'', y', y'' % assumed units: meters, seconds, Newtons, kg, radians global Cd Cl0 Clslope rho A S m Pa Pr g % these are defined globally so % they are available to the main script and all functions. % % Cdrag=Cd*Afront*rho/2/m drag force constant (Afront varies with attack angle) % Clift=Cl*Afloor*rho/2/m lift force constant (Cl varies with attack angle) % Pat=Pa+Pr*t; % current attack angle Afront=A+S*abs(sin(Pat)); % current frontal area includes floor projection due to attack angle Cdrag=Cd*Afront*rho/2/m; Cl=Cl0+Clslope*Pat; % current coefficient of lift Clift=Cl*S*rho/2/m; Fout = zeros(4,1); % initialize space Fout(1) = p(2); Fout(2) = -Cdrag*sqrt(p(2)^2 + p(4)^2)* p(2) - Clift*sqrt(p(2)^2 + p(4)^2)* p(4); Fout(3) = p(4); Fout(4) = -g -Cdrag*sqrt(p(2)^2 + p(4)^2)* p(4) + Clift*sqrt(p(2)^2 + p(4)^2)* p(2); end -------------------------------------------------------------------------------------------------------------- File #3, an event function which stops ode45 integration when the trajectory crosses the x axis ----------------------------------------------- MySaab96JumpEvent.m ------------------------------------------ function [value,isterminal,direction] = MySaab96JumpEvent(t,p)% when value is equal to zero, an event is triggered. %% For use with MySaab96Jump.m, Ver. 5, 3/31/2018 %% %% Copyright 2018, George (Ted) Yurkon %% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To %% view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative Commons, %% PO Box 1866, Mountain View, CA 94042, USA. %% %% In short, this work may be shared provided terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 %% International License are adhered to. %% % Setting isterminal to 1 stops the solver at the first event. % Setting isterminal to 0 allows the solver to compute past all events. % Set direction=0 if all zeros are to be computed (the default) % Set direction=+1 for only zeros where the event function is increasing. % Set direction=-1 for only zeros where the event function is decreasing. value = p(3); % when value = 0 (y=0), an event is triggered isterminal = 1; % terminate after the first event direction = -1; % only the zeros when decreasing end --------------------------------------------------------------------------------------------------------------