What Projective Angle Makes the Arc-Length of the Trajectory in a Resistive Media Maximum? A Reverse Engineering Approach

Haiduke Sarafian

The Pennsylvania State University, University College, York, USA.

DOI: 10.4236/ajcm.2021.112007                

Abstract

We consider the motion of a massive point-like projectile thrown with initial velocity with respect to horizontal in a two-dimensional vertical plane under the influence of gravity in a viscose media. Two different velocity-dependent resistive media models are considered—linear and quadratic. With an objective to utilizing a Computer Algebra System (CAS), specifically Mathematica [1] numerically we solve the corresponding equations of motions. For a set of compatible parameters characterizing viscose forces graphically we display comparing the trajectories explicitly showing the impact of the models. Utilizing the model-dependent trajectory equations numerically we evaluate their associated arc-lengths. What distinguishes our approach vs. the existing body of work is the notion of the “reverse engineering”. Meaning, utilizing our numeric data we establish their corresponding analytic counter parts. Ultimately, utilizing both outputs numerically and analytically we determine the matching initial projectile angles maximizing their respective arc-lengths.

Keywords

Velocity Dependent Resistive Media, Trajectory Arc-Length, Maximizing Angle, Computer Algebra System, Mathematica

1. Introduction

Motion of a point-like massive projectile thrown in vacuum with initial velocity in a vertical plane under the sole action of gravity is a theme of interest in introductory physics and engineering college courses [2] [3]. Although simple it sets the foundation to building more realistic applicable cases. The upgrade of the situation is a case where the projectile is thrown in a resistive, viscose media. Depending to the size of the projectile and characteristics of the media resistive forces are parameterized either as a velocity or velocity-squared entities. For microscale projectiles such as a bacteria linear velocity-dependent forces are used [4]. For larger scales e.g., a baseball, a parachute velocity-squared are considered [2]. Issues of interest in both cases include physical quantities such as time of flight and geometric quantities such as trajectories, range, maximum height etc. are addressed in the cited references. The next layer of complication includes the impact of the size, spin, buoyant as well as Magnus force—their collective impact reported [5]. References are resourceful yet none with one exception [6] addresses the objective of the current investigation—i.e., the initial projectile angle that maximizes the trajectory’s arc-length. Twenty-one years ago, the author published a report identifying the value of this special initial angle that maximizes the trajectory length in a non-resistive media, that is 56.46˚ [6]. This is a global unique angle; it is evaluated for a point-like projectile thrown in a vacuum. Considering the advances made within the last two decades in Computer Algebra System (CAS) that facilitates addressing non-trivial issues it is natural to revisit the same issues addressed in [6] embodying the resistive media. Therefore, the main objective and interest of our current work is to develop a forum to evaluating issues discussed in [6] with upgrades. In addition to the author’s own work there are ample resources [7] [8] addressing the impact of resistive forces, mostly analytic and in bits. Meaning, no one single article presents the linear and quadratic velocity-dependent viscous resistive forces side-by-side. In order avoiding repeating the existing body of work here we pursue tackling the problem entirely in numeric mode. Nonetheless, utilizing the numeric output we generate its useful analytic counterpart. To fulfill both objectives explicitly we apply Mathematica’s superb numeric and analytic built-in features.

As pointed out analysis of motion of a projectile in vacuum is trivial. Analysis in a velocity-dependent viscose media is somewhat non-challenging as well. The author has pushed the boundaries of the latter revealing the unknown envelope related facts [9]. Cited references addressing the numeric solutions include sections solving the same analytic equation of motion numerically by crafting a fourth order Runge-Kutta numeric algorithm. Crafting such codes during investigation de-focuses the objectives. Noting commonly used Computer Algebra Systems (CAS) embody built in algorithms assisting to solving differential equations; with no need to re-inventing the wheels. Consequently, helping to stay focused on the objectives.

Henceforth, we lay the foundation of our work on numeric approach assisting to stay focus conducive information that comparatively are not reported in the literature. Our work includes Mathematica codes for generating the numeric solutions for equations of motions for three cases: non-resistive, linear-velocity and velocity-squared viscose medias. The solutions are graphically compared assisting quantifying the impact of the modeled resistive forces. Utilizing the solutions, we calculate the arc-lengths of the trajectories. Ultimately empowering to identifying the angles maximizing them. Auxiliary information such as flight times etc. are calculated as well.

One last word before presenting our work—it is gratifying to show that our current numeric-based computation leading to the angle of interest for non-viscose media matches exactly the one reported twenty-one years ago that was derived analytically [6].

Our current work is composed of three sections. In addition to Sect. 1, Introduction, Motion and Goals, Sect. 2 embodies analytic equations of motions and their associated Mathematica codes conducive to their solutions. The trajectories and arc-lengths graphically are depicted yielding to the interested quantities.

2. Equations of Motions, Mathematica Codes

Following the objectives mentioned in Sect 1 Figure 1 depicts the relevant forces acting on a projectile projected in a vertical plane. In addition to gravity the viscous media exerts a velocity dependent resistive force. Since we are interested in two different models, linear and velocity-squared resistive forces, in Figure 1 these two forces are shown as ${\stackrel{\to }{F}}_{\text{resistive}}$ with a vector opposite to the orientation of the instantaneous velocity. According to the models the former is formulates as, ${\stackrel{\to }{F}}_{\text{resistive}}=-mn\stackrel{\to }{v}$, and the latter as, ${\stackrel{\to }{F}}_{\text{resistive}}=-mn{v}^2\hat{v}$. In these equations m is the mass of the projectile, n is the scale factor adjusting the strength of the resistive force, v is the speed of the projectile, and $\hat{v}$ is the unit velocity vector. The minus signs in front of the forces are the indicative of the orientation of the forces opposite to direction of the motion signifies by $\stackrel{\to }{v}$. Applying Newton’s law, ${\stackrel{\to }{F}}_{net}=m\stackrel{\to }{\ddot{r}}$, where $\stackrel{\to }{r}$ is the position vector of the projectile and super double dots indicates its second order derivative with respect to time, is the acceleration. ${\stackrel{\to }{F}}_{net}$ is the net force embodying the active forces acting on the projectile. With this simple formulation we include the impact of either one of the

Figure 1. Direction of motion of the projectile in the vertical plane is shown by velocity, $\stackrel{\to }{v}$. Gravity pulls, $m\stackrel{\to }{g}$ and the resistive force, ${\stackrel{\to }{F}}_{\text{resistive}}$ depicted as well.

resistive forces. We also notice the null value of the scale factor n associates with the motion in a non-resistive, vacuum case.

Because the projectile is projected in a vertical plane, the above-mentioned dynamic equation in the text is to be projected along two axes of a Cartesian coordinate system shown in Figure 1. For the case where the resistive viscous force is a liner-velocity dependent we write,

$\{\begin{array}{l}\ddot{x}\left(t\right)+n\stackrel{\dot{}}{x}\left(t\right)=0\\ \ddot{y}\left(t\right)+n\stackrel{\dot{}}{y}\left(t\right)+g=0\end{array}$ (1)

And for the scenario where the resistive force is velocity-squared it yields,

$\{\begin{array}{l}\ddot{x}\left(t\right)+n\stackrel{\dot{}}{x}\left(t\right)\sqrt{\stackrel{\dot{}}{x}{\left(t\right)}^2+\stackrel{\dot{}}{y}{\left(t\right)}^2}=0\\ \ddot{y}\left(t\right)+n\stackrel{\dot{}}{y}\left(t\right)\sqrt{\stackrel{\dot{}}{x}{\left(t\right)}^2+\stackrel{\dot{}}{y}{\left(t\right)}^2}=0\end{array}$ (2)

In (2) g is the gravity acceleration. Notice also as mentioned setting n = 0 in either one of (1), (2) yields the motion in a vacuum.

Because we would like to compare the trajectories of the projectiles under the same initial conditions, we consider the initial conditions,

$\stackrel{\to }{r}\left(t=0\right)=0\text{and}\stackrel{\to }{\stackrel{\dot{}}{r}}\left(t=0\right)=\left\{v_0\cos \left(\theta \right),v_0\sin \left(\theta \right)\right\}$ (3)

The first term of (3) shows the projectile is projected from the origin, the second term is it is horizontal and vertical components of the initial velocity. The v₀ is the initial speed, angle θ is the initial projectile angle with respect the x-axis. In a nutshell this is the angle we are interested evaluating it leading maximizing its associated arc-length. In the forthcoming sections we have shown how it has been determined.

To make the analysis as general as possible without losing the generalities we set the initial speed and the gravity acceleration unity, $v_0=g=1$. This helps focusing on the impact of the speed scale factor n on our analysis. We apply the same initial conditions for all three cases of interest. A compact Mathematica code is crafted addressing all three cases. However, for the sake of transparency and for the readers not quite fluent with Mathematica language three short codes are prepared.

Case 1. Impact of the linear-velocity dependent resistive force.

code 1

plotTrajectoryV1n={};

values1={v0->1,g->1,n->1.0}

initialx1={x1[0]==0,x1'[0]==v0 Cos[θ]}/.values1;

initialy1={y1[0]==0,y1'[0]==v0 Sin[θ]}/.values1;

eqx1=x1"[t]+n x1'[t]/.values1;

eqy1=y1"[t]+n y1'[t]+g/.values1;

solx1=DSolve[{eqx1==0,initialx1},x1[t],t]//Simplify;

soly1=DSolve[{eqy1==0,initialy1},y1[t],t]//FullSimplify;

plotTrajectoryV1=Table[ParametricPlot[{x1[t]/.solx1,y1[t]/.soly1}/.θ->φ, {t,0,2},PlotRange->{{0,1},{0,0.5}},GridLines->Automatic,AxesLabel->{"x","y"},PlotStyle->{If[Evaluate[n/.values1]==0.0,{Dashing[{0.01}],(*Black*)Hue[0.7φ]},{(*Black*)Hue[0.7 φ]}]}],{φ,π/6,π/2,π/10}];

combinedList1=AppendTo[plotTrajectoryV1n,plotTrajectoryV1];

sV1=Show[Transpose[{combinedList1〚1〛,combinedList1〚2〛}]]

As shown in Figure 2 the impact of the media resistive force being linear-velocity dependent is considerable. Irrespective to the chosen initial projectile angle for a “short” horizontal traversed distance the trajectory in the vacuum and the one in the media do overlap, then they split apart. Meaning, the resistive force at the beginning is less effective than further down the flight. The common overlapping features of these flights may be summarized as: the resistive force shortens the range, lowers the height making the trajectory lump side-ed conducive to non-symmetric trajectories.

Case 2. Impact of the velocity-squared dependent resistive force.

code 2

plotTrajectoryV2n={};

values2={v0->1,g->1,n->1.0};

initialx2={x2[0]==0,x2'[0]==v0 Cos[θ]}/.values2;

initialy2={y2[0]==0,y2'[0]==v0 Sin[θ]}/.values2;

eqx2=x2"[t]+n x2'[t]Sqrt[x2'[t]^2+y2'[t]^2]/.values2;

eqy2=y2"[t]+n y2'[t]Sqrt[x2'[t]^2+y2'[t]^2]+g/.values2;

tablesolxyθ=Table[NDSolve[{eqx2==0,eqy2==0,initialx2,initialy2},{x2[t],y2[t]},{t,0,2}],{θ,π/6,π/2,π/10}];

plotTrajectoryV2=Table[ParametricPlot[{x2[t],y2[t]}/.tablesolxyθ[[i]],{t, 0,2.},PlotRange->{{0,1},{0,0.5}},GridLines->Automatic,AxesLabel->{"x","y"}, PlotStyle->{If[Evaluate[n/.values2]==0.0,{Dashing[{0.01}],(*Black*)Hue[0.7 i]},{(*Black*)Hue[0.7 i]}]}],{i,1,Length[tablesolxyθ]}];

combinedListV2=AppendTo[plotTrajectoryV2n,plotTrajectoryV2];

sV2=Show[Transpose[{combinedListV2〚1〛,combinedListV2〚2〛}]] (Figure 3).

Figure 2. Dashed curves are the trajectories in vacuum, solids curves are the corresponding trajectories in viscus media indexed with scale factor n = 1.0 with the matching colors vs the dashed. Projected angles are: θ = π/6, π/6 + π/10, π/6 + 2π/10, π/6 + 3π/10 radians from shortest to the longest range.

Figure 3. Description is the same as Figure 2. The difference is the viscosity of the media is quadratic speed. Dashed curves are the trajectories in vacuum, solids curves are the corresponding trajectories with matching colors in viscus media with scale factor n = 1.0. Projected angles are:θ = π/6, π/6 + π/10, π/6 + 2π/10, π/6 + 3π/10 radians for the shortest to the longest range.

For the sake of comprehension in Figure 4 we put the graphic information side-by-side. This figure explicitly shows the impact of the velocity-dependent models on the trajectories.

Figure 4 shows the impact of the velocity-dependent resistive forces. The red curve on he left panel corresponds to the resistive force formulated as nmv, the purple on the right panel is the trajectory associated with nmv². The impact of the velocity-square case is hampering the resistivity force allowing the projectile to reaching higher height and a longer range. For the sake of comparison, the scaling factor for both cases is set to unity, e.g., n = 1.0. This figure shows as we commented in the earlier paragraph the trajectories irrelevant of the case of interest for a certain horizontal traversed distance are indistinguishable vs. the corresponding situation for the vacuum. Then at certain point they split. This is the general feature, and we may not make more quantified sentence because the splitting coordinate depends on the value of the scaling factor n.

We note had we applied the standard parametrization of the viscosity media coefficients for each model would have impacted the format of the resistive forces, as such the value of the velocity scale factor would have been different. That in turn would have altered the shape of the trajectories, most likely resulting reordering the shown trajectories. However, the mere objective of our exploratory investigation is to compare the impact of the resistive media for cases preserving the commentaries other than the speed dependencies.

We extend the analysis exploring additional information concerning the arc-length of the trajectories and the angles maximizing their lengths. Two issues need to be rectified. First, the arc-length ought to be formulated based on parametric, t-dependent coordinate of the projectile traversing the trajectory. Customary this is written as, $\mathcal{l}={\displaystyle \int \sqrt{\stackrel{\dot{}}{x}{\left(t\right)}^2+\stackrel{\dot{}}{y}{\left(t\right)}^2}\text{d}t}$. Where the dots are the derivatives of the coordinates with respect to t. Integration is performed over time

Figure 4. The left panel is the velocity-dependent scenarios, the right panel is the quadratic speed viscous cases. For both cases all the other parameters are kept the same.

variable, t. Secondly, the lower and upper limits of integration are zero and $t_{\text{flight}}$, respectively. Here, $t_{\text{flight}}$ is the time-of-flight of the projectile. The reader going over this article would appreciate to learn how Mathematica can evaluate the needed derivatives of the integrand $\mathcal{l}$. Noting, the coordinate of the projectile, {x(t), y(t)} are the numeric solutions of the differential equations of motions given in Code 1 and Code 2. These are the coordinates used plotting the associated trajectories; these are set of implicit tabulated numeric. However, what needed is their slopes, i.e., their derivatives with respect to t; Mathematica is capable evaluating these slopes!—and finally, by evaluating the integrand it evaluates the integration conducive to the arc-length. To evaluate the f_light, we set the ordinate of the trajectory zero forming a numeric equation and search for its roots. With this guided procedure we craft a Mathematica code. This allows the interested reader to run the code confirming our results.

Code 3

x1y1primetθ=Table[D[Evaluate[{x1[t],y1[t]}/.NDSolve[{eqx1==0,eqy1==0,initialx1,initialy1},{x1[t],y1[t]},{t,0,2}]],{t,1}],{θ,0.01π,0.9 π/2,0.005π}];

tableFlightTimeV1=Table[{θ,FindRoot[y1[t]/.NDSolve[{eqx1==0,eqy1== 0,initialx1,initialy1},{x1[t],y1[t]},{t,0,2.5}],{t,1}]},{θ,0.01π,0.9 π/2,0.005π}];

tableArcLengthV1=Table[NIntegrate[Sqrt[x1y1primetθ〚i,1,1〛^2+x1y1 primetθ〚i,1,2〛^2],{t,0,t/.tableFlightTimeV1[[i,2]]}],{i,1,Length[tableFlight TimeV1]}];

Max[tableArcLengthV1];

tableθ=Table[θ,{θ,0.01π,0.9 π/2,0.005π}];

TableForm[Transpose[{tableθ,tableArcLengthV1}]];

insetV1=TableForm[{{1.02101, 0.646601}, {1.03672, 0.646787}, {1.05243, 0.646849}, {1.06814, 0.646789}}];

s1=tableθArcLengthV1=ListPlot[Transpose[{tableθ,tableArcLengthV1}],AxesLabel->{"θ(rad)","ArcLength"},GridLines->Automatic,PlotStyle->Red, Epilog->Rectangle[{0.6,0.3},{1.1,1.0},insetV1],PlotRange->{{0,1.5},{0,1}}, AxesStyle->Directive[Arrowheads[0.05]]]

tableθArcLengthV1=Transpose[{tableθ,tableArcLengthV1}];

MaxtableθArcLengthV1=Select[tableθArcLengthV1,MemberQ[#,Max[tableArcLengthV1]]&];

The output of the given code is converted in a graph as depicted in Figure 5.

For the sake of accurately reading the coordinate of the depicted points in Figure 5, we created an insert table. The first and second columns are the initial projectile angle and the corresponding ArcLength, respectively. As such the angle maximizing the ArcLength is 60.30˚.

As Figure 5 shows the arc-lengths of the trajectories is a gradual increasing function of the initial projectile angles. It eventually passes a maximum value and then it is decreasing. The insert table assists identifying the maximum coordinate with high accuracy. The maximum coordinate is {1.05243, 0.646849}.

Following the same procedure, we craft a Mathematica code like code 3 to handling the second scenario, i.e., for the velocity-squared resistive media. The interested reader may utilize code 3 crafting the needed code. We skip displaying the code however the Figure 6 shows the output.

For comprehensive understanding a briefed version of the graphs is shown in Figure 7.

This figure shows the quadratic-speed dependent media is less resistive vs. the linear-velocity resistance case. As such the less resistive media allows longer trajectories conducive to longer arc-lengths. This is consistent with previous results shown in corresponding figures.

Note also analysis is based on point-like projectiles. If one considers finite-sized projectiles as intuitively expected the trajectories will have reverted orders, i.e., the linear-velocity trajectories will reach to higher heights conducive to longer arc-lengths.

Figure 5. Dots are the values of the arc-lengths, ℓ, for the trajectories in the viscous linear-velocity dependent resistive media. Horizontal axis is the initial projectile angles in radian.

Figure 6. Description of this figure is the same as Figure 5. The only difference is the curve is a generated for quadratic-speed dependent media.

Figure 7. The red and the purple dots are the arc-lengths of the trajectories vs. the initial projected angles. The red is for the linear-velocity and the purple is associated with the quadratic-speed dependent resistive media, respectively.

Summary of the collected data is presented in Figure 8. The dots are the coordinates of the maximum angles vs. the characteristics of the viscous media. The first dot at the left with a zero abscissa corresponds to the vacuum. The other two with velocity-character 1 and 2 correspond to linear and quadratic speed-dependent viscous cases.

Thus far as shown we have accomplished the set goals. In short, for a chosen resistive media applying a CAS numerically we identified the projectile angle that maximizes the arc-length of associated trajectory. The applied methodology is numeric and as such differential equations describing the motions are solved numerically. However, as pointed out in the abstract we wish to apply “reverse engineering” to generate analytic expressions describing trajectory related features. In another words, having the numeric features of the issues at hand aiming to generate their analytic equivalence. Here we show our approach for one case, namely the linear-speed dependent resistive case. The quadratic-speed case easily maybe followed similarly. Data shown in Figure 5 is the arc-length of the trajectory vs. the initial projectile angles. By trial and error, we search for a function best fit the data. Our search identifies a cubic polynomial. The coefficients of the fitted function are given and the plot of the fitted function and its overlap with the data are shown in Figure 9.

fitArcLengthV1=a+b θ+c θ²+d θ³

{a->0.0267161,b->1.64468,c->-1.43487,d->0.406906}

Figure 8. Plot of angles maximizing the arc-lengths of the projectiles in viscous media characterized with the velocity-dependent character 0, 1 and 2. These charters are indicatives of vacuum, linear and quadratic-speed dependent resistive media, respectively.

Figure 9. The arc-length of the projectile vs. the initial angles is shown in red dots, this is the same as Figure 5. The black dashed curve is the cubic fitted analytic function.

As shown the fitted function is just perfect. This is what we claim is “reverse engineering”. That is, we bypassed the regular routine solving the equation of motion analytically, and yet at the end just by utilizing the output of its numeric solution we were able to create analytic function as if we have solved the problem analytically. And at the end to show its usefulness utilizing the fitted function we seek for the angle maximizing its associated arc-length. We set the derivative of the arc-length to zero and search for its root. This results 56.73˚ and is just off by a 3.5˚ vs. its numeric counter version.

3. Discussions and Conclusions

We set out a list of objectives investigating information about the trajectories of projectiles projected in a vertical plane under the action of gravity pull and two different modes of velocity-dependent viscous medias. Globally speaking aside from analyzing the physics of the problems on hand technology wise the main objective is to utilize the ever-growing popularity of Computer Algebra System (CAS), specifically Mathematica generating the needed information. As it is discussed and shown, for the first scenario with a linear-velocity viscous case the equation of motion analytically trivially is solvable. However, for the second case where the viscosity makes the media velocity-squared resistive the equation of motion despite of its look becomes quite challenging. The sited references show how the authors by inventing the wheels have solved the equations numerically applying the popular numeric methods. Codes are useful however de-focuses the objective. Most of the CASs have built in sophisticated numeric methods. Familiarity with these CASs helps avoiding inventing the wheels facilitating a forum to merely focus on the science of the problem. With this objective we have applied the implicit Runge-Kutta numeric method of Mathematica conducive to the solutions of the needed equations of motions. We have provided set of Mathematica codes so that interested readers familiar with this CAS may run the codes reproducing our results. The codes are robust and the control parameters such as velocity scale factor may be adjusted to generate alternate cases of interest. To make the analysis helpful we converted equation of motions to graphics. These figures help to form a visual understanding about the problem on hand. For the second scenario where the viscosity is velocity-squared dependent the arc-length expression needed to be integrated numerically as well. This challenging procedure is accomplished with a one-line code—attesting to the amazing computation power of Mathematica. Meaning the integrand is a list of numeric pairs yet the integration is accomplished without a glitch. Utilizing this information finally we numerically evaluated the initial projectile angle maximizing the arc-length of the associated projectile. For a comprehensive understanding we have supplied all the information graphically, making our analysis comprehensive. Twenty-one years ago [6] the author addressed a similar problem for the case of no resistive media, i.e., a vacuum. The problem was solved analytically conducive to the angle of interest. It is gratifying to justify that our current numeric solution matches our reported angle when the viscous impact is turned off. We also show how the “reverse engineering” was used to generate an analytic expression for the arc-length of the projectile as if the equation of the motion were solved analytically.

As a final comment the author would like to attract the reader’s attention that the entire investigation that is comprised of analytic, numeric and graphics all are done in one single file avoiding to creating separate files for different aspects of the objectives. Interested readers may find [10] [11] resourceful.

Acknowledgements

The author gracefully acknowledges the John T. and Paige S. Smith Professorship funds for completing and publishing this work.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	(2020) Mathematica V12.1.1. http://wolfram.com/
[2]	Halliday, D., Resnick, R. and Walker, J. (2013) Fundamentals of Physics Extended. 10th Edition, John Wiley and Sons, New York.
[3]	Tippler, P. and Mosca, G. (2008) Physics for Scientists and Engineers. 6th Edition, Freeman and Company, New York.
[4]	Kambe, T. (2007) Elementary Fluid Mechanics. World Scientific Publishing Co. Pte. Ltd., New Jersey. https://doi.org/10.1142/5895
[5]	Sarafian, H. (2015) Impact of the Drag Force and the Magnus Effect on the Trajectory of a Baseball. World Journal of Mechanics, 5, 49-58. https://doi.org/10.4236/wjm.2015.54006
[6]	Sarafian, H. (1999) On Projectile Motion. The Physics Teacher, 37, 86. https://doi.org/10.1119/1.880184
[7]	Sketches, N. (2008) Mathematical Physics2. https://www.ph.ed.ac.uk/
[8]	Timmerman, P. and van der Weele, J. (1999) On the Rise and Fall of a Ball with Linear and Quadratic Drag. American Journal of Physics, 67, 538. https://doi.org/10.1119/1.19320
[9]	Sarafian, H. (2020) Envelope of Family of Angled Projectiles and Its Universal Geometric Characteristics. American Journal of Computational Mathematics, 10, 425-430. https://doi.org/10.4236/ajcm.2020.103023
[10]	Wolfram, S. (1996) Mathematica Book. 3rd Edition, Cambridge University Press, Cambridge.
[11]	Sarafian, H. (2019) Mathematica Graphics Examples. 2nd Edition, Scientific Research Publishing, Wuhan.