(* PROJECTILE *)
M=165/100000kg; (* Mass *)
r=55/20000m; (* Radius *)
BC=0.78lb (in^(-2)); (* Ballistic Coefficient *)
cB=4M/(BC π (r^2)); (* BC Factor *)

(* SETTINGS *)
e0=38kg (m^2)(sec^(-2)); (* Muzzle Energy *)
v0=(Sqrt[2 e0/M]); (* Muzzle Velocity *)
α1=5degree; (* Inclination Angle *)
h = 1m; (* Height *)
sz=400m; (* Scope Zero *)
sh=7cm; (* Scope Height *)
sx=60cm; (* Scope to Barrel Distance, x Axis *)
w=0m/sec; (* Horizontal Wind Speed *)
ts=500sec; (* Simulation Duration *)

(* CONSTANTS *)
cs=3432/10m/sec; (* Speed of Sound *)
G=667384*^-16 (m^3)(kg^(-1))(sec^(-2)); (* Newtons Constant *) 
mE= 597219*^19kg; (* Earth Mass *)
rE=6378000m; (* Earth Radius *)
Ah=π (r^2); Av=Ah; (* Cross Section Area *)
hs=7800m; (* Scale Height *)
ρ0=12041/10000 kg (m^(-3)); (* Density of Air *)
ρ[h_]:=ρ0 (E^(-h/hs)); (* Density by Height *)
g[h_]:=G mE/(((rE+h))^2) ; (* Gravity by Height *)

(* DIMENSIONS *)
kg=1;m=1; cm=m/100; mm=m/1000; sec=1; degree=π/180;  ft=0.3048m; in=0.0254m; lb=0.45359237kg; 

(* DRAG COEFFICIENT *)
p=Interpolation[{0.24,0.23,0.22,0.21,0.2,0.19,0.19,0.21,0.24,0.32,0.5,0.6,0.64,0.67,0.68,0.67,0.656,0.642,0.628,0.614,0.6}]; (* Fitting Curve, Round Nose Diabolo *)
cW[v_]:=-p[10v/cs+1]cB; (* Drag Coefficient Diabolo *)
P=Interpolation[{0.48,0.485,0.49,0.495,0.5,0.51,0.54,0.6,0.67,0.75,0.82,0.9,0.96,1,1.01,1.02,1.02,1.02,1.015,1.01,1.005}]; (* Fitting Curve, Sphere *)
CW[v_]:=-P[10v/cs+1]; (* Drag Coefficient Ball *)

(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];

{"v0"->N[v0]}
τ=Quiet[FindRoot[Evaluate[ry[t]/.sol/.t->time] [[1]] , {time,10}]] (* Time at impact, in seconds *)
d=Quiet[{"distance"-> Evaluate[rx[t]/.sol/.t->τ[[1]][[2]]][[1]] } ](* Impact distance, in meters *)
hMax={"max height" ->Quiet[FindMaximum[Evaluate[ry[time]/.sol], {time,0}]]} (* Maximum heigth *)
H=Quiet[{"impact height"-> Evaluate[ry[t]/.sol/.t->time/.τ] [[1]] }](* Impact height, in meters *)
e=Quiet[{"energy"-> M( ((Evaluate[rx'[t]/.sol/.t->time/.τ] [[1]] )^2)+ ((Evaluate[ry'[t]/.sol/.t->time/.τ] [[1]] )^2))/2}] (* Impact energy, in Joule *)

(* CW DIAGRAM *)
Plot[{-cW[v],-CW[v]}, {v,0,2cs}, Frame->True, ImageSize->540, AspectRatio->2/8, PlotRange->{{0,2cs},{0,1.4}}]

plot=Manipulate[sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
Grid[{
{Show[
ListPlot[Table[{Flatten[Evaluate[{rx[t]/.sol,ry[t]/.sol}]]},{t,0,tm,0.01}], Frame->True, ImageSize->540, AspectRatio->1/2, PlotRange->{{0,d[[1]][[2]]},{0,hMax[[1]][[2]][[1]]*1.2}}, GridLines->{{5,sz, Evaluate[rx[time/.hMax[[1]][[2]][[2]]]/.sol][[1]],Evaluate[rx[tm]/.sol][[1]]},{hMax[[1]][[2]][[1]], Evaluate[ry[tm]/.sol][[1]], h}}, PlotTheme->"Classic"],

Graphics[{Gray, Dashed, InfiniteLine[{{(-sx)Cos[α1]+(h+sh)Sin[α1],(sx)Sin[α1]+(h+sh)Cos[α1]},{sz, Quiet[Evaluate[ry[Quiet[Tsz/.Quiet[FindRoot[Evaluate[rx[Tsz]/.sol]-sz, {Tsz,(v0/sz)}]]]]/.sol]][[1]]}}]}]
]},
{Grid[{
{Quiet["t"->  {Evaluate[tm]}] ,
Quiet["x"->  Evaluate[rx[tm]/.sol]] ,
Quiet["e"-> {M( ((Evaluate[rx'[tm]/.sol] [[1]] )^2)+ ((Evaluate[ry'[tm]/.sol] [[1]] )^2))/2}]},

{Quiet["vt"->{(Sqrt[( ((Evaluate[rx'[tm]/.sol] [[1]] )^2)+ ((Evaluate[ry'[tm]/.sol] [[1]] )^2))])}],
Quiet["dr"->{(Evaluate[ry[tm]/.sol] [[1]])-
(Evaluate[rx[tm]/.sol] [[1]]*((Quiet[Evaluate[ry[time/.Quiet[FindRoot[Evaluate[rx[t]/.sol/.t->time] [[1]]-sz , {time,10}]]]/.sol]][[1]])-((sx)Sin[α1]+(h+sh)Cos[α1]))/((sz)-((-sx)Cos[α1]+(h+sh)Sin[α1])))-((Quiet[Evaluate[ry[time/.Quiet[FindRoot[Evaluate[rx[t]/.sol/.t->time] [[1]]-sz , {time,10}]]]/.sol]][[1]])-(sz)*((Quiet[Evaluate[ry[time/.Quiet[FindRoot[Evaluate[rx[t]/.sol/.t->time] [[1]]-sz , {time,10}]]]/.sol]][[1]])-((sx)Sin[α1]+(h+sh)Cos[α1]))/((sz)-((-sx)Cos[α1]+(h+sh)Sin[α1])))}],
Quiet[ "Δα"->
{-(α1-(α2/.FullSimplify[Solve[((Evaluate[rx'[tm]/.sol] [[1]] +Evaluate[ry'[tm]/.sol] [[1]] )/(Sqrt[Evaluate[((rx'[tm])^2)/.sol] [[1]] +Evaluate[((ry'[tm])^2)/.sol] [[1]]]))==Sin[α2]+Cos[α2]&&α2≤(π/2)&&α2≥-(π/2), α2]][[1]]) )(180/π)}]}}]}}],
{{tm,τ[[1]][[2]],"t"},0,τ[[1]][[2]]}]

Export["ballistik.avi", plot, "FrameRate"->200]

energyslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["ek"-> M( ((Evaluate[rx'[tm]/.sol] [[1]] )^2)+ ((Evaluate[ry'[tm]/.sol] [[1]] )^2))/2] ,
{tm,0,τ[[1]][[2]]}]

velocityslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["vt"-> √( ((Evaluate[rx'[tm]/.sol] [[1]] )^2)+ ((Evaluate[ry'[tm]/.sol] [[1]] )^2))] ,
{tm,0,τ[[1]][[2]]}]

vxslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["vx"->  Evaluate[rx'[tm]/.sol] [[1]] ] ,
{tm,0,τ[[1]][[2]]}]

vyslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["vy"->  Evaluate[ry'[tm]/.sol] [[1]] ] ,
{tm,0,τ[[1]][[2]]}]

dropslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["dr"-> {(Evaluate[ry[tm]/.sol] [[1]])-
(Evaluate[rx[tm]/.sol] [[1]]*((Quiet[Evaluate[ry[time/.Quiet[FindRoot[Evaluate[rx[t]/.sol/.t->time] [[1]]-sz , {time,10}]]]/.sol]][[1]])-((sx)Sin[α1]+(h+sh)Cos[α1]))/((sz)-((-sx)Cos[α1]+(h+sh)Sin[α1])))-((Quiet[Evaluate[ry[time/.Quiet[FindRoot[Evaluate[rx[t]/.sol/.t->time] [[1]]-sz , {time,10}]]]/.sol]][[1]])-(sz)*((Quiet[Evaluate[ry[time/.Quiet[FindRoot[Evaluate[rx[t]/.sol/.t->time] [[1]]-sz , {time,10}]]]/.sol]][[1]])-((sx)Sin[α1]+(h+sh)Cos[α1]))/((sz)-((-sx)Cos[α1]+(h+sh)Sin[α1])))} ] ,
{tm,0,τ[[1]][[2]]}]

xslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["px"->  Evaluate[rx[tm]/.sol] [[1]] ] ,
{tm,0,τ[[1]][[2]]}]

yslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet["py"->  Evaluate[ry[tm]/.sol] [[1]] ] ,
{tm,0,τ[[1]][[2]]}]

alphaslide=Manipulate[
(* DIFFERENTIAL *)
sol=Quiet[NDSolve[{Sign[rx'[t]]*cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Ah ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*(rx'[t]+w))/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×rx''[t],rx'[0]==v0×Cos[α1],rx[0]==0 ,-M×g[ry[t]]+Sign[ry'[t]]×cW[(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)])] Av ρ[ry[t]]×((((((rx'[t]+w))^2)+((ry'[t])^2))*ry'[t])/(Sqrt[(((rx'[t]+w))^2)+((ry'[t])^2)]))==M×ry''[t],ry'[0]==v0×Sin[α1],ry[0]==h},{rx,ry},{t,0,ts}]];
(* PLOT *)
Quiet[ "Δα"->
-(α1-(α2/.FullSimplify[Solve[((Evaluate[rx'[tm]/.sol] [[1]] +Evaluate[ry'[tm]/.sol] [[1]] )/(Sqrt[Evaluate[((rx'[tm])^2)/.sol] [[1]] +Evaluate[((ry'[tm])^2)/.sol] [[1]]]))==Sin[α2]+Cos[α2]&&α2≤(π/4)&&α2≥-(π/4), α2]][[1]]) )(180/π)] ,
{tm,0,τ[[1]][[2]]}]

Export["energyslide.avi", energyslide, "FrameRate"->200]

Export["velocityslide.avi", velocityslide, "FrameRate"->200]

Export["vxslide.avi", vxslide, "FrameRate"->200]

Export["vyslide.avi", vyslide, "FrameRate"->200]

Export["dropslide.avi", dropslide, "FrameRate"->200]

Export["xslide.avi", xslide, "FrameRate"->200]

Export["yslide.avi", yslide, "FrameRate"->200]

Export["alphaslide.avi", alphaslide, "FrameRate"->200]

(* Simon Tyran, Vienna. Yukterez, 11.3.2015 *)

Quit[] (* Ruger Black Hawk Magnum .177 @ 32 Joule *)

(* PROJECTILE *)
M=87/100000kg; (* Mass *)
r=45/20000m; (* Radius *)

(* SETTINGS *)
e0=32kg (m^2)(sec^(-2)); (* Muzzle Energy *)
v0=(Sqrt[2 e0/M]); (* Muzzle Velocity *)
α1=0degree; (* Inclination Angle *)
h = 3m; (* Height *)
ts=500sec; (* Simulation Duration *)

Quit[] (* 100 m Match .177 @ 7.5 Joule *)
(* PROJECTILE *)
M=55/100000kg; (* Mass *)
r=45/20000m; (* Radius *)

(* SETTINGS *)
e0=7.5kg (m^2)(sec^(-2)); (* Muzzle Energy *)
v0=(Sqrt[2 e0/M]); (* Muzzle Velocity *)
α1=1.35degree; (* Inclination Angle *)
h = 1m; (* Height *)
ts=500sec; (* Simulation Duration *)

Quit[] (* CG Test *)

(* PROJECTILE *)
M=165/100000kg; (* Mass *)
r=55/20000m; (* Radius *)

(* SETTINGS *)
e0=38kg (m^2)(sec^(-2)); (* Muzzle Energy *)
v0=(Sqrt[2 e0/M]); (* Muzzle Velocity *)
α1=0degree; (* Inclination Angle *)
h = 2m; (* Height *)
sz=10m; (* Scope Zero *)
sh=5cm; (* Scope Height *)
sx=0cm; (* Scope to Barrel Distance, x Axis *)
w=0m/sec; (* Horizontal Wind Speed *)
ts=500sec; (* Simulation Duration *)

Quit[] (* Benjamin Trail NP XL 1100 .22 @ 38 Joule *)

(* PROJECTILE *)
M=165/100000kg; (* Mass *)
r=55/20000m; (* Radius *)

(* SETTINGS *)
e0=38kg (m^2)(sec^(-2)); (* Muzzle Energy *)
v0=(Sqrt[2 e0/M]); (* Muzzle Velocity *)
α1=5degree; (* Inclination Angle *)
h = 1m; (* Height *)
sz=400m; (* Scope Zero *)
sh=7cm; (* Scope Height *)
sx=60cm; (* Scope to Barrel Distance, x Axis *)
w=0m/sec; (* Horizontal Wind Speed *)
ts=500sec; (* Simulation Duration *)