I ran into an error in the calculation of H3 at the bottom of the loop. Looking into it, it was suspicious to me that H2, H3 and H4 are overwritten there having been calculated previously, and I suspected that the H2 and H3 which are used there to calculate H3 and H4, respectively, should actually refer to the previously calculated H2 and H3. Changing those bottom H2, H3 and H4 to other names, I was able to generate a plot with 3 distinct trajectories (no idea if they are correct). Here is the code with my changes:
clear all
clc
H0 = 8;
theta_0=26;
g=9.81;
dAC=16;
eR=.94 ;
H3=0;
v0=30;
err = 5;
counter=1;
while err>0.008
%Iteration Control
v0 = v0 - .01;
counter=counter+1;
%initial velocities in both directions
v0x=v0*cosd(theta_0); %Trig
v0y=v0*sind(theta_0); %Trig
%Point B = Maximmum height along trajectory path 1
H1=H0+v0y^2/(2*g);
%trajectory path 1
H2=H0+dAC*tand(theta_0)-g/(2*v0x^2)*dAC^2;
%velocities in x and y direction at max height
v1x=v0x; %Zero Acceleration
v1y=(sqrt(v0y^2-2*g*(H1-H0)));
%velocity just prior to hitting wall
v2x=v1x;
v2y=sqrt(v1y^2+2*g*((H0+H1)-(H0+H2)));
v2=sqrt(v2x^2+v2y^2);
%NEW NEW
%H2 = H1 - v2y^2/(2*g);
%velocity just after hitting the wall
v3x=v2x*eR;
v3y=v2y;
v3=sqrt(v3x^2+v3y^2);
theta_3_ref=atand(v3y/v3x); % Trig: Reference Triangle
theta_3_abs=theta_3_ref+180; % Absolute Angle (Quadrant "3")
%velocity just prior to hitting ground
v4x=v3x;
v4y=(sqrt(v3y^2-2*g*(H3-H2)));
v4=sqrt(v4x^2+v4y^2);
theta_4_ref=atand(v4x/v4y);
theta_4_abs=theta_4_ref+180;
%Distance between wall and ground
A=(-g/(2*v3x^2));
B=tand(theta_3_abs);
C=(1/(2*g)*(v4^2-v3^2));
%Quadratic formula for distance between wall and ground
%%based on energy Conservation 2-3 022& Trajectory Eq. 2-3 equations
dCD=(-B+sqrt(B^2-4*A*C))/(2*A);
%trajectory path 2
H3=H2+dCD*tand(theta_3_abs)-g/(2*v3x^2)*dCD^2;
%velocity just after hitting the ground
v5x=v4x;
v5y=v4y*eR;
v5=sqrt(v5x^2+v5y^2);
theta_5_ref=atand(v5x/v5y);
theta_5_abs=theta_4_ref+90;
%distance between ground and return height
dDA=-dAC-dCD;
%trajectory path 3
H4=H3+dDA*tand(theta_5_abs)-g/(2*v5x^2)*dDA^2;
%Flight time between point A & B
time_AB=(H1-v0)/-g;
%Flight time between point A & C
time_AC=dAC/v0;
%Flight time between point C & D
time_CD=abs(dCD/v4);
%Flight time between point D & return height
time_D_return=abs(dDA/v5);
%Update Error
err=abs(H4-H0);
x=[0:.01:16];
H22=H0+x*tand(theta_0)-g/(2*v0x^2)*x.^2;
x1=[16:-.01:9.05];
H33=H2+x1*tand(theta_3_abs)-g/(2*v3x^2)*x1.^2;
x2=[9.05:-.01:0];
H44=H3+x2*tand(theta_5_abs)-g/(2*v5x^2)*x2.^2;
plot(x,H22,x1,H33,x2,H44)
end
v0
v0y
v0x
v1y
v1x
v2
v2y
v4
dCD
dDA
H1
H2
H3
H4
time_AB
time_AC
time_CD
time_D_return
counter
