Hi,
I am currently working on a university assignment where I have to calculate the height at which a package should be dropped from a helicopter. The package is dropped and after 3 seconds a parachute opens. It is given that the package must be traveling slower than 20mph (8.94ms^-1) for the package not to be damaged.
The assignment asks me to: estimate the height to drop the package. find time taken to reach the ground, find the displacement from a target the package must be dropped (x-axis displacement),find the resultant velocity that the package contacts the ground with. The results must be plotted.
The way I've approached the question is firstly, calculate the absolute minimum height required for a safe landing (for emergency situations say), then secondly I have calculated a recommended drop height based on when the resultant velocity comes with 0.1ms^1 of the y axis terminal velocity.
The code works exactly how I want it to, there is just a lot of code, can anyone suggest a way to streamline the code or can anyone see any flaws in my logic?
%MODEL VELOCITY AND DISPLACEMENT.
%Constant boundary conditions
syms Vx Vy Xx Xy %Create displacement and velocity symbols.
m=265; %Package mass [kg].
g=9.81; %Gravity [ms^-2].
p=1.225; %Air density at sea level [kgm^-3].
Vxheli=14.3053; %Helicopter velocity [ms^1].
Vsafe=8.9408; %Max safe landing velocity [ms^1].
%Boundary conditions parachute closed.
Axx=0.7; %Area facing x-x air flow.
Ayy=1; %Area facing y-y air flow.
Cdx=0.9; %x-axis drag coefficient.
Cdy=1.3; %y-axis drag coefficient.
%Define time periods for analysis (0.1 second steps).
t1=linspace(0,3,31); %0-3 seconds, parachute closed.
t2=linspace(3,12,91); %3-12 seconds, parachute closed.
%ODE's modelling package displacement and velocity (F=ma).
eqn1 = '-(0.5*p*Axx*Cdx)*Vx^2=m*DVx'; %x-axis velocity.
eqn2 = 'm*g-(0.5*p*Ayy*Cdy)*Vy^2=m*DVy'; %y-axis velocity.
eqn3 = '-(0.5*p*Axx*Cdx)*DXx^2=m*D2Xx'; %x-axis displacement.
eqn4 = 'm*g-(0.5*p*Ayy*Cdy)*DXy^2=m*D2Xy'; %y-axis displacement.
%SOLVE FOR VELOCITY AND DISPLACEMENT
%Solve for x-axis velocity (parachute closed).
inits1 = 'Vx(0)=Vxheli'; %When t=0 seconds, Vx=Vxheli.
Vx1=dsolve(eqn1,inits1,'t1'); %Solve ODE for Vx1(t1).
Vx1=eval(vectorize(Vx1)); %Evaluate Vx1(t1) for all values of t1.
Vxclosed=Vx1(end); %Vxclosed=velocity at 3 seconds.
%Solve for y-axis velocity (parachute closed).
inits2 = 'Vy(0)=0'; %When t=0 seconds, Vy=0ms^1.
Vy1=dsolve(eqn2,inits2,'t1'); %Solve ODE for Vy1(t1).
Vy1=eval(vectorize(Vy1)); %Evaluate Vy1(t1) for all values of t1.
Vyclosed=Vy1(end); %Vyclosed=velocity at 3 seconds.
%Solve for x-axis displacement (parachute closed).
inits3 = 'Xx(0)=0'; %When t=0, Xx=0m.
inits4 = 'DXx(0)=Vxheli'; %When t=3, DXx=Vxheli.
Xx1=dsolve(eqn3,inits3,inits4,'t1'); %Solve ODE for Xx1(t1).
Xx1=eval(vectorize(Xx1)); %Evaluate Xx1(t1) for all values of t1.
Xxclosed=Xx1(end); %Xxclosed=displacement at 3 seconds.
%Solve for y-axis displacement (parachute closed).
inits5 = 'Xy(0)=0'; %When t=0, Xy=0m.
inits6 = 'DXy(3)=Vyclosed'; %When t=3, DXy=Vyclosed.
Xy1=dsolve(eqn4,inits5,inits6,'t1'); %Solve ODE for Xy1(t1).
Xy1=eval(vectorize(Xy1)); %Evaluate Xy1(t1) for all values of t1.
Xyclosed=Xy1(end); %Xyclosed=displacement at 3 seconds.
%Boundary conditions parachute open.
Axx=35; %Area facing x-x air flow.Ayy=50; %Area facing y-y air flow.Cdx=1.25; %x-axis drag coefficient.Cdy=1.75; %y-axis drag coefficient.%Solve for x-axis velocity (parachute open).
inits7 = 'Vx(3)=Vxclosed'; %When t=3 seconds, Vx=Vxclosed.
Vx2=dsolve(eqn1,inits7,'t2'); %Re-solve ODE for Vx2(t2).
Vx2=eval(vectorize(Vx2)); %Evaluate Vx2(t2) for all values of t2.
%Solve for y-axis velocity (parachute open).
inits8 = 'Vy(3)=Vyclosed'; %When t=3 seconds, Vy=Vyclosed.
Vy2=dsolve(eqn2,inits8,'t2'); %Re-solve ODE for Vy2(t2).
Vy2=eval(vectorize(Vy2)); %Evaluate Vy2(t2) for all values of t2.
%Solve for x-axis displacement (parachute open).
inits9 = 'Xx(3)=Xxclosed'; %When t=3, Xxclosed=Xxclosed.
inits10 = 'DXx(3)=Vxclosed'; %When t=3, DXx=Vxclosed.
Xx2=dsolve(eqn3,inits9,inits10,'t2'); %Re-solve ODE for Xx2(t2).
Xx2=eval(vectorize(Xx2)); %Evaluate Xx2(t2) for all values of t2.
%Solve for y-axis displacement (parachute open).
inits11='Xy(3)=Xyclosed'; %When t=3, Xyclosed=0m.
inits12='DXy(3)=Vyclosed'; %When t=3, DXy=Vyclosed.Xy2=dsolve(eqn4,inits11,inits12,'t2'); %Re-solve ODE for Xy2(t2).
Xy2=eval(vectorize(Xy2)); %Evaluate Xy2(t2) for all values of (t2).
%EMERGENCY DROP.
%Calculate emergancy drop distance/height/time.
V1=sqrt((Vx1.^2)+(Vy1.^2)); %Resultant velocity parachute closed.
V2=sqrt((Vx2.^2)+(Vy2.^2)); %Resultant velocity parachute open.
Emergcolumn=find(V2<(Vsafe),1,'first'); %Column where V<Vsafe.
EmergV=V2(1,Emergcolumn); %Find the velocity at time taken.
Emergtime=t1(end)+((Emergcolumn-1)*0.1) %Calculate time taken for V<Vsafe.
EmergXx=Xx2(1,Emergcolumn) %Target distance for emergancy landing.
EmergXy=Xy2(1,Emergcolumn) %Minimum height for emergancy landing.
%Plot y-axis displacement time graph.
figure('Name','Y Axis Displacement of Package','NumberTitle','off') %Plot in figure window.
plot(t1,Xy1,t2,Xy2) %Plot Xy against t.
grid %Display grid lines.
line([Emergtime Emergtime], [0 EmergXy],'LineStyle', '--','Color', 'r') %Mark V at time.
line([0 Emergtime], [EmergXy EmergXy],'LineStyle', '--','Color', 'r' ) %Mark V at time.legend('Parachute Closed','Parachute Open','Emergancy Drop','Location','southeast') %Create legend.
xlabel('Time [s]') %Label x-axis.
ylabel('Y Axis Displacement of Package [m]') %Label y-axis.
title('Package Y Axis Displacement Against Time') %Title graph.
%Plot resultant velocity time graph.
figure('Name','Resultant Velocity of Package','NumberTitle','off') %Plot in figure window.plot(t1,V1,t2,V2,'r') %Plots V against t.
grid %Display grid lines.line([Emergtime Emergtime], [0 EmergV],'LineStyle', '--','Color', 'r') %Mark V at time.line([0 Emergtime], [EmergV EmergV],'LineStyle', '--','Color', 'r' ) %Mark V at time.legend('Parachute Closed','Parachute Open','Emergancy Drop') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('Resultant Velocity of Package [mms^{-1}]') %Label y-axis.title('Resultant Velocity Against Time') %Title graph.%CONTROLLED DROP.
%Calculate controlled drop distance/height/time.
Vyterminal=sqrt((m*g)/(0.5*p*Cdy*Ayy)); %Calc Vyterminal.
Ctrlcolumn=find(V2<(Vyterminal+0.1),1,'first'); %Column where V approx Vyterminal.
CtrlV=V2(1,Ctrlcolumn); %Find the velocity at taken.
Ctrltime=t1(end)+((Ctrlcolumn-1)*0.1) %Calc time taken for V approx Vyterminal.
CtrlXx=Xx2(1,Ctrlcolumn) %Minimum height for safe landing.
CtrlXy=Xy2(1,Ctrlcolumn) %Minimum height for safe landing.%Plot y-axis displacement time graph.figure('Name','Y Axis Displacement of Package','NumberTitle','off') %Plot in figure window.plot(t1,Xy1,t2,Xy2) %Plot Xy against t.grid %Display grid lines.line([Ctrltime Ctrltime], [0 CtrlXy],'LineStyle', '--','Color', 'r') %Mark V at time.line([0 Ctrltime], [CtrlXy CtrlXy],'LineStyle', '--','Color', 'r' ) %Mark V at time.legend('Parachute Closed','Parachute Open','Controlled Drop','Location','eastoutside') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('Y Axis Displacement of Package [m]') %Label y-axis.title('Package Y Axis Displacement Against Time') %Title graph.%Plot resultant velocity time graph.figure('Name','Resultant Velocity of Package','NumberTitle','off') %Plot in figure window.plot(t1,V1,t2,V2,'r') %Plots V against t.grid %Display grid lines.line([Ctrltime Ctrltime], [0 CtrlV],'LineStyle', '--','Color', 'r') %Mark V at time.line([0 Ctrltime], [CtrlV CtrlV],'LineStyle', '--','Color', 'r' ) %Mark V at time.legend('Parachute Closed','Parachute Open','Controlled Drop') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('Resultant Velocity of Package [mms^{-1}]') %Label y-axis.title('Resultant Velocity Against Time') %Title graph.%PLOT VELOCITY AND DISPLACEMENT.
%Plot x-axis velocity time graph.
figure('Name','X Axis Velocity of Package','NumberTitle','off') %Plot in figure window.plot(t1,Vx1,t2,Vx2,'r') %Plot Vx against t.
grid %Display grid lines.legend('Parachute Closed','Parachute Open') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('X Axis Velocity of Package [mms^{-1}]') %Label y-axis.title('Package X Axis Velocity Against Time') %Title graph.%Plot y-axis velocity time graph.
figure('Name','Y Axis Velocity of Package','NumberTitle','off') %Plot in figure window.plot(t1,Vy1,t2,Vy2,'r') %Plot Vy against t.
grid %Display grid lines.legend('Parachute Closed','Parachute Open') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('Y Axis Velocity of Package [mms^{-1}]') %Label y-axis.title('Package Y Axis Velocity Against Time') %Title graph.%Plot x-axis displacement time graph.
figure('Name','X Axis Displacement of Package','NumberTitle','off') %Plot in figure window.plot(t1,Xx1,t2,Xx2) %Plot Xx against t.
grid %Display grid lines.legend('Parachute Closed','Parachute Open','Location','southeast') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('X Axis Displacement of Package [m]') %Label y-axis.title('Package X Axis Displacement Against Time') %Title graph.%Plot y-axis displacement time graph.figure('Name','Y Axis Displacement of Package','NumberTitle','off') %Plot in figure window.plot(t1,Xy1,t2,Xy2) %Plot Xy against t.grid %Display grid lines.legend('Parachute Closed','Parachute Open','Location','southeast') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('Y Axis Displacement of Package [m]') %Label y-axis.title('Package Y Axis Displacement Against Time') %Title graph.%Plot resultant velocity time graph.figure('Name','Resultant Velocity of Package','NumberTitle','off') %Plot in figure window.plot(t1,V1,t2,V2,'r') %Plots V against t.grid %Display grid lines.legend('Parachute Closed','Parachute Open') %Create legend.xlabel('Time [s]') %Label x-axis.ylabel('Resultant Velocity of Package [mms^{-1}]') %Label y-axis.title('Resultant Velocity Against Time') %Title graph.
Best Answer