% DEMO of using Matlab ODE solvers % % This script m-file solves for the initial velocity needed for a pojectile % to hit a target a specified distance from the lauch point. The program % includes air drag and height dependent gravity. % % There are many ways to solve this problem and we will use the Matlab ODE % solvers with event detection (the event in this case is the z-value % hitting zero). % global targdist g0 dgdh rho0 scaleht mass xarea cd g0 = -9.806; % m/sec^2 rho0 = 1.29; % kg/m^3 %----- Get the inout that we need ----- prompt = {'Distance (km)', 'Error (m)', 'Launch Angle (deg)', 'Mass (kg)', ... 'CrossSection Area (m^2)', 'Drag Coefficient', ... 'AirScale height (m)', 'dg/ht (m/sec^2)/m'}; title = 'Ballistic Trajectory Parameters'; nlines = 1; defaults = {'100','1','45','500','0.01','0.1','7500.e6', '3.0786e-6'}; % Get user inputs and save in the appropriate valuyes result = inputdlg(prompt,title,nlines, defaults,'on'); % targdist = eval(result{1})*1000.0; terr = eval(result{2}); ltheta = eval(result{3})*pi/180; mass = eval(result{4}) ; xarea = eval(result{5}) ; cd = eval(result{6}) ; scaleht = eval(result{7}) ; dgdh = eval(result{8}) ; % % Based on no-drag, no dght get initial velocity vx = sqrt(-targdist*g0/(2*tan(ltheta))); vz = vx*tan(ltheta); vtot = vx/cos(ltheta); dvxdD = sqrt(-g0/(2*tan(ltheta)*targdist))/2; dvdD = dvxdD/cos(ltheta); tmax = targdist/vx*2*(1+cd); % Set the ODE options options = odeset('AbsTol',[terr 1 1 1],'Events','Lec18_hit'); %%% Loop until target distance is hit done = 0; prevDD = 0; iter = 0; while done == 0 % Set the initial conditions for the ODE iter = iter + 1; y0 = [0.0; 0.0; vx; vz]; [t,y,te,ye,ie] = ode23(@Lec18_bacc,[0:1:tmax],y0,options); % Check the time if te/tmax > 0.5 , tmax = 2*te; end % See how close we got to target dd = targdist - ye(1); fprintf('Result te %f ye %f dd %f \n',te, ye(1), dd); if abs(dd) < terr done = 1; % We are finished else % Adjust the initial velocity. Initially we use the analytic value % for the directivative but this updated once we start if prevDD == 0 vtot = vtot + dvdD*dd; prevDD = dd; else dvdD = dvdD*(1+dd/prevDD); vtot = vtot + dvdD*dd; end vx = vtot*cos(ltheta); vz = vtot*sin(ltheta); end if vtot > 3.e8 fprintf('\nVelocity faster than speed of light!\n') done = 1; end end % % Output the final results fprintf('TARGET DISTANCE: %8.3f km, within %6.1f m, Launch Angle %6.3f deg\n', ... targdist/1000, terr, ltheta*180/pi) fprintf('MASS %7.1f kg, XArea %7.4f m**2, Drag Coefficient %7.3f\n', ... mass, xarea, cd) fprintf('GRAVITY %8.5f m/s**2 %11.4e (m/sec**2)/m\n', ... g0, dgdh) fprintf('RHO AIR %8.5f kg/m**3, Scale Height %11.5e m\n', ... rho0, scaleht) fprintf('Initial Velocity needed %8.3f n/s\n',vtot) fprintf('Calculation took %d iterations with final error %8.2f m\n',iter,dd) % Now animate (make as function) Lec18_animate(t,y); fprintf('To show animation: animate(t,y)\n'); %