Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["GUIKit`"];

xzm[a_,g_,h_,V10_,α_]:=Module[{so,Vh,tfin,de},Vh=If[α==2/19,(10^(-α)) V10 ((z[t])^α),1.42 V10 (((z[t]/366))^α)];
de={
x'[t]     == Vmh[t],
z'[t]     == - Vmv[t],
Vmh'[t]== a (Vh-Vmh[t])(Sqrt[((( Vh-Vmh[t]))^2)+((Vmv[t])^2)]),
Vmv'[t] == g- a Vmv[t](Sqrt[(((Vh-Vmh[t]))^2)+((Vmv[t])^2)]),
x[0]==0,z[0]==h, Vmh[0]==0, Vmv[0]==0
};
so=Quiet[NDSolve[de,{x,z,Vmh,Vmv},{t,0,Infinity},
Method->{EventLocator,Event:>Re[z[t]]}]];
tfin = InterpolatingFunctionDomain[First[z /. so]][[1,-1]];
{so[[1]],tfin}]

goldenSearch[f_ ,x_,aa_ ,bb_,eps_]:=Module[{h=(bb-aa)/2,a=aa,b=bb,c,d,fc,fd,fcd=2,xm},c=b-(Sqrt[5.]-1)h;fc=f/.x-> c;While[h>eps, If[c>a+h,d=c;fd=fc;c=a+b-d;fc=f/.x-> c,d=a+b-c;fd=f/.x-> d];If[fc<fd,b=d,a=c;c=d;fc=fd];h=(b-a)/2;fcd=fcd+1];xm=a+h;{xm,h,f/.x-> xm,fcd}]

am={0.0021,0.0026,0.0042,0.0057,0.0885,0.176,0.079};

alpha={2/19,1/7};

missileVelocities[missileType_,terrain_,V10_,h_]:=Module[{α,a,g=9.81,tif,xz,vh,vt,vfh,vfv,vft,af,thm,ttm,hm,vhm,ht,vth,vtv,vtm,gsh,gst},
α=alpha[[terrain]];a=am[[missileType]];
{xz,tif}=xzm[a,g,h,V10,α];
{vfh,vfv}=Re[{Vmh[tif],Vmv[tif]}]/.xz;
af= Round[ArcTan[vfv/vfh]/Degree];vft=Sqrt[vfh^2+vfv^2];
vh=Re[Vmh[t]]/.xz;
gsh=goldenSearch[-vh,t,0 ,tif,.001];thm=gsh[[1]];
{hm,vhm}=Re[{z[thm],Vmh[thm]}/.xz];
 vt =Sqrt[Re[Vmh[t]]^2+Re[Vmv[t]]^2]/.xz;
  gst=goldenSearch[-vt,t,0 ,tif,.001];ttm=gst[[1]];
  {ht,vth,vtv}=Re[{z[ttm],Vmh[ttm],Vmv[ttm]}]/.xz;
 vtm=Sqrt[vth^2+vtv^2];
Print["missileType = ",missileType,",  terrain = ",terrain , ",  V10 =  ", V10, " m/s,  h = ", h," m"];
Print["terminal horizontal velocity = ",NumberForm[vfh,{4,1}], " m/s"];
Print["terminal total velocity = ",NumberForm[vft,{4,1}], " m/s"];
Print["angle of implact with ground = ", af,  " degrees"];
Print["maximum horizontal velocity = ",NumberForm[vhm,{4,1}], " m/s", " attained at a height of ", NumberForm[hm ,{4,2}], " m "];
Print["maximum total velocity = ",NumberForm[vtm,{4,1}], " m/s", " attained at a height of  ", NumberForm[ht ,{4,2}], " m "]; {vfh,{vft,af},{vhm,hm},{vtm,ht}};
]

missileVelocitiesA[a_,terrain_,V10_,h_]:=Module[{α,g=9.81,tif,xz,xh,xht,xv,xvt,vz,vzt,vh,vht,vt,vtt,vfh,vfv,vft,at,att,af,thm,ttm,hm,vhm,ht,vth,vtv,vtm,gsh,gst},
α=alpha[[terrain]];
{xz,tif}=xzm[a,g,h,V10,α];
{vfh,vfv}=Re[{Vmh[tif],Vmv[tif]}]/.xz;
af=ArcTan[vfv/vfh]/Degree;vft=Sqrt[vfh^2+vfv^2];
vh=Re[Vmh[t]]/.xz;
vht[s_]:=vh/.t-> s;
xh=Re[x[t]]/.xz;
xht[s_]:=xh/.t-> s;
xv=Re[z[t]]/.xz;
xvt[s_]:=xv/.t-> s;
vz=Re[Vmv[t]]/.xz;
vzt[s_]:=vz/.t-> s;
gsh=goldenSearch[-vh,t,0 ,tif, .1^12];thm=gsh[[1]];
{hm,vhm}=Re[{z[thm],Vmh[thm]}/.xz];
 vt =Sqrt[Re[Vmh[t]]^2+Re[Vmv[t]]^2]/.xz;
vtt[s_]:=vt/.t-> s;
 at =ArcTan[(Re[Vmh[t]]/Re[Vmv[t]])/.xz]/Degree;
 att[s_]:=at/.t-> s;
  gst=goldenSearch[-vt,t,0 ,tif, .1^12];ttm=gst[[1]];
  {ht,vth,vtv}=Re[{z[ttm],Vmh[ttm],Vmv[ttm]}]/.xz;
 vtm=Sqrt[vth^2+vtv^2];
Print["a = ",a,",  terrain = ",terrain , ",  V10 =  ", V10, " m/s,  h = ", h," m"];
Print["terminal horizontal velocity = ",vfh, " m/s"];
Print["terminal total velocity = ",vft, " m/s"];
Print["angle of implact with ground = ", af,  " degrees"];
Print["maximum horizontal velocity = ",vhm, " m/s", " attained at a height of ", hm , " m "];
Print["maximum total velocity = ",vtm, " m/s", " attained at a height of  ", ht, " m "]; {{vfh,tif},{vft,af},{vhm,hm,thm},{vtm,ht,ttm},{xht,xvt,vht,vzt,vtt,att}}
]