function Fd = drag_3DOF(vx,vy,vz,h)

%	This function gives a shape to the drag curve
%	This function is in the FIRST ORDER stage of development.

%	VARIABLES and UNITS
%
%	conventions:	Large arrays are in all caps
%						Vectors and other variables are in mixed case, as appropriate
%
%	h					:Height above sea level, in m.
%	vx, vy, vz		:Component velocities of the vehicle, in m/s.
%	v					:Velocity vecotr for the vehicle, in m/s.
%	cs_area			:Frontal cross-sectional area of the vehicle, in m2.
%	Cd					:Coefficient of drag, unitless.
%	M					:Mach number the vehicle is at, unitless.
%	gamma				:Coefficient of specific heats, unitless.
%	R					:Gas constant for air, in J/kg/K.



%	CONSTANTS
cs_area = .026;			% m2
gamma = 1.4005;		% unitless
R = 287.10;				% J/kg/K



%	MAIN PROGRAM

v = [vx vy vz];

%a = sqrt(gamma * R * temperature(h));
a = sqrt(gamma * R * 273);
M = norm(v) / a;


if (norm(v) <= 100)
   Cd = .15;
   
elseif (norm(v) <= 200)
   Cd = .35;
   
elseif (norm(v) <= 250)
   Cd = 1.10;
   
elseif (norm(v) <= 300)
   Cd = 1.50;
   
elseif (norm(v) <= 500)
   Cd = 1.15;
   
elseif (norm(v) <= 800)
   Cd = .95;
   
else 
   Cd = .80;
   
end


%Fdx = sign(vx) * Cd * 1/2 * density(h) * vx^2 * cs_area;
%Fdy = sign(vy) * Cd * 1/2 * density(h) * vy^2 * cs_area;
%Fdz = sign(vz) * Cd * 1/2 * density(h) * vz^2 * cs_area;
Fdx = sign(vx) * Cd * 1/2 * 3 * vx^2 * cs_area;
Fdy = sign(vy) * Cd * 1/2 * 3 * vy^2 * cs_area;
Fdz = sign(vz) * Cd * 1/2 * 3 * vz^2 * cs_area;


Fd = [Fdx Fdy Fdz];