% We want to demonstrate here the use of numerical analysis with
% Euler for real world problems.
% 
% The problem is to shoot a canon ball with speed v and angle
% a. This can, of course, be computed exactly. But we restrict
% ourselfs to numerical solutions only.
% 
% We fix the speed v=10 m/s and approximate gravity with g=10
% m/s^2.
>v=10; g=10;
% Now we plot a specific ballistic curve for a canon ball fired
% to 45° with the well known formulas for the curve neglecting
% friction.
% 
% Time runs from 0 s to 2 s.
>t=0:0.1:2;
% The movement in x-direction (height) is cos(a)*v*t, but we must
% take care to compute 45° in radians first.
% 
% The movement in y-direction is sin(a)*v*t minus the falling
% of g*t^2/2.
>xplot(cos(45°)*v*t,sin(45°)*v*t-g*t^2/2);
% We notice that we are not interested in the negative y-values
% and clip the plot to the postive quadrant.
% 
% To copy the previous input, one can use the clipboard or recall
% the previous command with Alt-CursorUp.
>setplot(0,10,0,10); xplot(cos(45°)*v*t,sin(45°)*v*t-g*t^2/2);
% We now compute, when the ball has its heighest altitude. There
% is the funciton fmax, which computes the maximum of another
% function or expression in a given range.
>tmax=fmax("sin(45°)*v*x-g*x^2/2",0,2)
     0.707107 
% How high is the maximal height?
>x=tmax; sin(45°)*v*x-g*x^2/2
          2.5 
% To make things a little bit more comfortable, we introduce two
% functions computing the point (sx,sy) at time t.
% 
% We take v and g fromt he global values, but add a as a parameter
% to these functions.
>function sx (t,a)
$global v;
$return cos(a)*v*t;
$endfunction
>function sy (t,a)
$global v,g;
$return sin(a)*v*t-g*t^2/2;
$endfunction
% This makes computing the maximal height more comfortable.
>a=45°; sy(fmax("sy(x,a)",0,2),a)
          2.5 
% We wish to compare several angles. Thus we generate a vector
% of angles from 5° to 85° and generate all curves for all these
% angles simultanously.
% 
% To do so, we need a as a column vektor and t as a row vector.
% Combining both yields a matrix, with a curve in each row for
% each angle.
>a=(5:5:85)°; a=a';
>xplot(sx(t,a),sy(t,a));
>setplot(0,10,0,10); xplot(sx(t,a),sy(t,a));
% This gave a nice picture.
% 
% Lets us now compute the distance of the contact of the ball
% with the ground. We use the simplest of methods, the bijection
% method. Since 0 is already a solution of sy=0, we start a little
% bit off.
% 
% We see the zero of vy. The additional parameter a (after the
% semicolon) is passed from bisect to sy.
% 
% Note: Our function impact will only work for scalar a, not for
% vectors a.
>function distance (a)
$return sx(bisect("sy",0.00001,2;a),a);
$endfunction
% What happens, if we tage 50°?
>distance(50°)
      9.84808 
% We get about 9.85m far.
% 
% Let us plot the impact distance with varying angles.
>a=(1:1:89)°;
% To compute the impact function to all angles, we use "map".
>xplot(deg(a),map("distance",a));
% We can also compute the optimal angle and its distance.
>amax=fmax("distance",1°,89°); deg(amax), distance(amax),
           45 
           10 
% Of course, it is easy to invert the function sx, i.e., compute
% the time t for the value x. We do not need to do that numerically.
% 
% So we write a function that computes the height of the ball
% at distance x, if the angle is a.
>function height (x,a)
$global v;
$return sy(x/(cos(a)*v),a);
$endfunction
% Test ist with 45° from 0 to 10.
>fplot("height",0,10;45°);
% We now want to get as high as possible in 5m distance. Again,
% we use fmax, but the argument x is now the second parameter
% to height, the angle.
>degprint(fmax("height(5,x)",1°,89°))
63°26'5.82''
% We have to shoot about 63° (not 45° as expected).
% 
% Now, we make a function that computes the angles that yields
% the maximal heights for given distances.
>function heighestangle (xx)
$global v;
$return fmax("height(xx,x)",1°,89°);
$endfunction
% Compute various maximal heights and plot them together with
% some shots.
>x=0.2:0.1:10; y=map("heighestangle",x);
>setplot(0,10,0,10); xplot(x,height(x,y));
>a=(5°:5°:85°)'; hold on; xplot(sx(t,a),sy(t,a)); hold off;
% It the curve of maximal heights a parabola?
% 
% Yes it is. To demosntrate, we compute the interpolation polynomial
% of second degree through the know values at -10, 0 and 10.
>d=interp([-10,0,10],[0,5,0]);
% Then we plot the polynomial over the current view in color 13.
% Indeed, we get the same curve.
>hold on; color(13); fplot("interpval([-10,0,10],d,x)",0,10); color(1); hold off;
% Next, we ask how to reach the point x=2, y=2. There seem to
% be two solutions. We get both by searching the angle below and
% the angle above the one with maximal height.
>a=heighestangle(2); degprint(a)
78°41'24.24''
>a1=bisect("height(2,x)-2",0,a); degprint(a1)
51°31'33.49''
>a2=bisect("height(2,x)-2",a,pi/2); degprint(a2)
83°28'26.51''
% Let's plot both solutions to one plot.
>t=0:0.01:2;
>setplot(0,5,0,5); xplot(sx(t,a1),sy(t,a1));
>hold on; xplot(sx(t,a2),sy(t,a2)); hold off;
% Assume we shoot a ball in 80° and another one in 40°. Can we
% hit the first with the secon in the air?
% 
% For a start, we plot the two paths.
>fplot("height(x,80°)",0,4); hold; fplot("height(x,40°)",0,4); hold;
% Now we compute the point, where both paths meet and compute
% the time difference.
>xd=bisect("height(x,80°)-height(x,40°)",0.001,5)
      3.07202 
>xd/(v*cos(80°))-xd/(v*cos(40°))
      1.36808 
% Now we introduce friction. We assume the friction force is proprotional
% to c*v^2, where v is the speed.
% 
% We need to solve a differential equation, using the "runge"
% function. The state x of the ball at any time t is 4-dimensional
% and contains two point coordinates and two speed coordinates.
% We need a function to compute x'=f(t,x).
% 
% We set x=[sx,sy,sx',sy'].
>function f(t,x,c=0.005)
$global g;
$v=sqrt(x[3]*x[3]+x[4]*x[4]); rho=c*v;
$return [x[3],x[4],-rho*x[3],-g-rho*x[4]]
$endfunction
% We need to pass a vector of times t, where x will be computed.
% "runge" returns a matrix s, with one value for s in each column.
% We may plot the path of the ball from the first two rows of
% s.
>t=0:0.01:2; s=runge("f",t,[0,0,cos(45°)*v,sin(45°)*v]);
>setplot(0,10,0,10); xplot(s[1],s[2]);
>hold on; xplot(sx(t,45°),sy(t,45°)); hold off;
% Clearly, we do not get as far as before.
% 
% Now we look at the speed over a larger time. The speed will
% tend to some value, where friction is equal to gravity.
>t=0:0.01:20; s=runge("f",t,[0,0,cos(45°)*v,sin(45°)*v]);
>xplot(t,sqrt(s[3]^2+s[4]^2));
% We wish to compute the angle of maximal width for this problem.
% It will no longer be exactly 45°.
% 
% First we define a function returning the height at time t.
>function height (t,a)
$global v;
$s=runge("f",[0,t],[0,0,cos(a)*v,sin(a)*v],20);
$return s[2,2];
$endfunction
>fplot("height",0,2;45°);
% Then we seek the time when the height is 0 for any angle a,
% and return the width at that time.
>function width (a)
$t=bisect("height",0.001,2;a);
$global v;
$s=runge("f",[0,t],[0,0,cos(a)*v,sin(a)*v],20);
$return s[1,2];
$endfunction
>width(45°)
      9.62578 
% Now we maximize the width over the angles.
% 
% This computation takes some time (approx. 3 seconds on my modern
% Notebook)
% 
% We find that the optimal angle is slightly less than 45°.
>amax=fmax("width",30°,60°); degprint(amax)
44°42'25.31''
>
>