function [initial_v]=long_jump_full(density,desired_distance)

%Long_jump_full.m will run through a set of initial velocities for a given
%density until the length of the modeled jump is approximately equal to
%distance

p=density;                                  %set p=density for subfunction purpose
v0_guess=11;                                %Guess for vo is 10
distance=7;                                 %Approximate value for distance with vo=10
theta0=pi/8; xo=0;yo=0;                     %Initial values
y0=[theta0;v0_guess;xo;yo];                                 %Initial value vector
opts=odeset('events',@long_jump_landing,'reltol',10^-13)    %set up event finder
difference=abs(distance-desired_distance);

while difference > .0001                            %Quit when x(tf)~9.8
    y0=[theta0;v0_guess;xo;yo];                     %Initial values for iteration
    [t,y,te]=ode45(@long_jump,[0,2],y0,opts,p);    %numerical solution
    num_of_elements=length(t);                      %how many elements are in t
    dis=y(:,3);                                     %x(tf)
    distance=dis(num_of_elements);                  %distance
    difference=abs(distance-desired_distance);       %difference between dis and 9.8
    v0_guess=v0_guess+0.0001;                        %increment velocityed
    
end

initial_v=v0_guess-0.0001;                           %initial velocity that gives x(tf)=8.90


%%subfunctions for call to ode and event finder

function long_jump=f(t,y,p)

%This function will attempt to find the trajectory of a long jumper using
%the differential equations provided in the online textbook problem #7.16

g=9.8;      %acceleration due to gravity
c=0.72;     %drag coefficient
m=80;       %mass of long jumper in kg
s=0.5;      %effective surface area of jumper

long_jump=[-g*cos(y(1))/y(2);-(c*p*s/(2*m))*(y(2))^2*((cos(y(1)))^2+sin(y(1))^2)-g*sin(y(1));y(2)*cos(y(1));y(2)*sin(y(1))];

function [gstop,isterminal,direction]=long_jump_landing(t,y,p)

gstop=y(4);     %Measure when y=0--jumper has landed
isterminal=1;   %stop the numerical solving
direction=-1;   %jumper is falling