Chapter 2 
Realistic Projectile Motion

The effect of air resistance

Sample code 2.1.1

P=400.0 (* in unit of Watt *);
v0=4.0 (* in unit of m/s *);
Deltat=0.1 (*in unit of second *);
NStep = (200/Deltat); (* Number of total time step *)
m=70.0;(* in unit of 70 kg *)
(*A quantity must be declared as real number if it is not intrisically integer *)

v[0]=v0;
For [i=0,i≤NStep,i++,
v[i+1]=v[i] + P*Deltat/(m*v[i]);
];(* end for *)
tvf=Table[{i*Deltat,v[i]},{i,0,NStep}];
Print[tvf];
(*ListPlot[tvf]*)

Sample code 2.1.2: Add some feature to your ListPlot. 
You can modify the ListPlot by adding in the option: Joinned --> True; AxesLabel -> {"time (s)", "velocity (m/s)'} via: 
ListPlot[tv, Joined -> True, AxesLabel -> {velocity (m/s), time (s)}]

ListPlot[tvf,Joined->True,AxesLabel->{"time (s)","velocity (m/s)"}]

Sample code 2.1.3:

P=400.0 (* in unit of Watt *);
v0=4.0 (* in unit of m/s *);
Deltat=0.1 (*in unit of second *);
NStep = 200/Deltat; (* Number of total time step *)
m=70.0;(* in unit of 70 kg *)
A =0.33 (* area, in unit of m^2 *);
coef=1(* drag coefficient *);
rho=1.29 (*in unit of kg/m^3*);
v[0]=v0;
For [i=0,i≤NStep,i++,
v[i+1]=v[i] + P*Deltat/(m*v[i])-Deltat*coef*rho*A*v[i]^2/(2m);
];(* end for *)
tvr=Table[{i*Deltat,v[i]},{i,0,NStep}];
ListPlot[{tvr},Joined->True,AxesLabel->{"time (s)","velocity (m/s)"},PlotRange->All]

Sample code 2.1.4:

g=9.81 (*in unit of m/s^2*);
v0=700.0 (* in unit of m/s, magnitude of initial velocity *);
x0=0 (* in unit of m*);y0=0 (* in unit of m*);
Deltat=0.1 (*in unit of second *);
NStep = IntegerPart[(40000/Deltat)+1]; (* Number of total time step *)
B2overm=4*10^(-5);(* B2/m, in unit of m^(-1) *)
v[0]=v0;x[0]=x0;y[0]=y0;

angle0=30*(Pi/180);(* in unit of radian *)
vx[0]=v0*Cos[angle0];vy[0]=v0*Sin[angle0];

For [i=0,i≤NStep,i++,
	v[i]=Sqrt[vx[i]^2+vy[i]^2];
	x[i+1]=x[i]+vx[i]*Deltat;
	vx[i+1]=vx[i]-B2overm*v[i]*vx[i]*Deltat;
	y[i+1]=y[i]+vy[i]*Deltat;
	vy[i+1]=vy[i]-g*Deltat-B2overm*v[i]*vy[i]*Deltat;

	If[y[i+1]<0,
	ilast=i+1;
	Print["y[i+1]<0 ",y[i+1]<0];Break[]
	](*end if *)

];(* end for *)
txy=Table[{x[i]/1000,y[i]/1000},{i,0,ilast}];
ListPlot[{txy},Joined->True,AxesLabel->{"y (km)","x (km)"},PlotRange->All]

Sample code 2.1.5: Modify code 2.1.4 to generate a familiy of curve on the same graph, with each curve associated with a different initial angle.

g=9.81 (*in unit of m/s^2*);
v0=700.0 (* in unit of m/s, magnitude of initial velocity *);
x0=0 (* in unit of m*);y0=0 (* in unit of m*);
Deltat=0.1 (*in unit of second *);
NStep = IntegerPart[(40000/Deltat)+1]; (* Number of total time step *)
B2overm=4*10^(-5);(* B2/m, in unit of m^(-1) *)

For[j=6,j≤11,j++,

v[0]=v0;x[0]=x0;y[0]=y0;
angle0=( j*5)*(Pi/180);(* in unit of radian *)
vx[0]=v0*Cos[angle0];vy[0]=v0*Sin[angle0];
For [i=0,i≤NStep,i++,
v[i]=Sqrt[vx[i]^2+vy[i]^2];
x[i+1]=x[i]+vx[i]*Deltat;
vx[i+1]=vx[i]-B2overm*v[i]*vx[i]*Deltat;
y[i+1]=y[i]+vy[i]*Deltat;
vy[i+1]=vy[i]-g*Deltat-B2overm*v[i]*vy[i]*Deltat;

If[y[i+1]<0,
ilast=i+1;
Print["y[i+1]<0 ",y[i+1]<0];Break[];
](*end if *)
];(* end for i*)

txyj[j]=Table[{x[i]/1000,y[i]/1000},{i,0,ilast}];
];(*end for j*)
txy=Table[txyj[j],{j,6,11}];
ListPlot[txy,Joined->True,AxesLabel->{"y (km)","x (km)"},PlotRange->All]
(*ListPlot[{txyj[1],txyj[2]},Joined->True,AxesLabel->{"y (km)","x (km)"},PlotRange->All]*)

Sample code 2.1.6: Curve with Isothermal approximation

Y0=10^4;(* in unit of m*)
g=9.81 (*in unit of m/s^2*);
v0=700.0 (* in unit of m/s, magnitude of initial velocity *);
x0=0 (* in unit of m*);y0=0 (* in unit of m*);
Deltat=0.1 (*in unit of second *);
NStep = IntegerPart[(40000/Deltat)+1]; (* Number of total time step *)
B2overm=4*10^(-5);(* B2/m, in unit of m^(-1) *)

For[j=0,j≤1,j++,

v[0]=v0;x[0]=x0;y[0]=y0;
angle0=(35 + j*10)*(Pi/180);(* in unit of radian *)
vx[0]=v0*Cos[angle0];vy[0]=v0*Sin[angle0];
For [i=0,i≤NStep,i++,
v[i]=Sqrt[vx[i]^2+vy[i]^2];
x[i+1]=x[i]+vx[i]*Deltat;
vx[i+1]=vx[i]-(Exp[-y[i]/Y0])*B2overm*v[i]*vx[i]*Deltat;
y[i+1]=y[i]+vy[i]*Deltat;
vy[i+1]=vy[i]-g*Deltat-(Exp[-y[i]/Y0])*B2overm*v[i]*vy[i]*Deltat;

If[y[i+1]<0,
ilast=i+1;
Break[];
](*end if *)
];(* end for i*)

txyjthermal[j]=Table[{x[i]/1000,y[i]/1000},{i,0,ilast}];
];(*end for j*)
txythermal=Table[txyjthermal[j],{j,0,1}];
ListPlot[txythermal,Joined->True,AxesLabel->{"y (km)","x (km)"},PlotRange->All]
(*ListPlot[{txyj[1],txyj[2]},Joined->True,AxesLabel->{"y (km)","x (km)"},PlotRange->All]*)

Sample code 2.1.7: Plotting simulataneously curves with and without thermal approximation for curves with initial angle 35 and 45 degree.

ListPlot[{txyj[7],txyj[9],txyjthermal[0],txyjthermal[1]},PlotStyle->{Black,Blue,Dashed,Dashed},Joined->True,AxesLabel->{"y (km)","x (km)"}]