The Fastest Motion of an Object in a Curved Earth Tunnel — Revisiting the Earth Tunnel Problem Using Python Visualization Toolkit

1. Model Establishment

To solve the problem of curved earth tunnels, we first establish a model : assume that the earth is a uniform sphere with a radius of R, a mass of M, and a density of ρ, then dig a (curved) tunnel of unlimited shape, regard any object in the tunnel as a mass point m, r is the distance between the mass point and the center of the earth, and establish a polar coordinate system with the center of the earth as the origin (as shown in Figure 1). Ignore the centrifugal inertia of the earth's rotation, the Coriolis force, the size of the tunnel, the friction in the tunnel, etc., and discuss the movement of the mass point under the condition of only considering the action of universal gravitation.

2. Python Visualization Toolkit Analysis Modeling

2.1. Theoretical Modeling Process

According to Gauss's theorem, take any point (r,θ) inside the earth, then the gravity of the particle at that point points to the center of the sphere, and the magnitude is

$$F={\displaystyle \frac{4\pi \rho G}{3}}mr$$

Let $k={\displaystyle \frac{4\pi \rho G}{3}}$, assuming that the infinite distance outside the earth is the potential energy zero point, then the gravitational potential energy of a point (r,θ) inside the earth is:

$$V\left(r\right)={\displaystyle \underset{0}{\overset{r}{\int }}kmr\mathrm{d}r}={\displaystyle \frac{1}{2}}km{r}^2$$

Then, the energy conservation equation for the process of moving from a point on the earth's surface to a point inside the earth (r₀ ,θ₀) is listed as follows:

$${\displaystyle \frac{1}{2}}m{v}^2+{\displaystyle \frac{1}{2}}km{r}^2={\displaystyle \frac{1}{2}}km{R}^2$$

It can be solved

$$v=\sqrt{k\left(R^2-r^2\right)}$$

At the same time, according to the representation of velocity in polar coordinate system

$$v=\sqrt{{\left({\displaystyle \frac{dr}{dt}}\right)}^2+{\left(r{\displaystyle \frac{d\theta }{dt}}\right)}^2}=\sqrt{{\left({\displaystyle \frac{dr}{d\theta }}\right)}^2+r^2}{\displaystyle \frac{d\theta }{dt}}$$

the time T required to travel from one point A on the Earth's surface to another point B on the Earth's surface.

$$T={\displaystyle {\int }_A^{B}dt}={\displaystyle {\int }_A^B\frac{1}{v}\sqrt{{\left(\frac{dr}{d\theta }\right)}^2+r^2}d\theta ={\displaystyle {\int }_A^B\frac{\sqrt{{\left(\frac{dr}{d\theta }\right)}^2+r^2}}{\sqrt{k\left(R^2-r^2\right)}}d\theta }}(r\ne R)$$

To find out where the minimum value of T is obtained, let $f\left(r\right)={\displaystyle \frac{\sqrt{{\left(\frac{dr}{d\theta }\right)}^2+r^2}}{\sqrt{R^2-r^2}}}$(where $r\ne R$), that is, to solve the minimum point (extreme point) of f(r). According to the Euler-Lagrange equation, if the function T is required to obtain an extreme value, the function f (r) needs to satisfy the equation

$${\displaystyle \frac{d}{d\theta }}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{dr}{d\theta }\right)}}\right]-{\displaystyle \frac{\partial f}{\partial r}}=0$$

To solve the differential equation more simply, we can use ${\displaystyle \frac{\partial f}{\partial \theta }}=0$ to transform it into

$${\displaystyle \frac{dr}{d\theta }}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{dr}{d\theta }\right)}}\right]-f=C(C\text{ is a constant})$$

Converted to

$${\displaystyle \frac{dr}{d\theta }}=-r\sqrt{{\displaystyle \frac{1}{C^2}}{\displaystyle \frac{r^2}{R^2-r^2}}-1}(\mathrm{in}r\ne R)$$

Assuming that θ₀ is the direction angle of the particle at the initial moment, then solving the transformed differential equation, we can finally obtain the analytical expression of the fastest descent line inside the earth in the polar coordinate system:

① $\theta ={\theta }_0-\arcsin \left(C\sqrt{{\displaystyle \frac{R^2-r^2}{r^2}}}\right)+{\displaystyle \frac{C}{\sqrt{C^2+1}}}\arctan \left[\sqrt{C^2+1}{\displaystyle \frac{\sqrt{\frac{R^2-r^2}{r^2}}}{\sqrt{1-C^2\left(\frac{R^2-r^2}{r^2}\right)}}}\right]$（$r\ne R$）

② $\theta ={\theta }_{0}or\theta ={\theta }_0-\pi \left(1-{\displaystyle \frac{C}{\sqrt{C^2+1}}}\right)$($r=R$, are the endpoint values of the fastest curve )

Due to the domain of the inverse sine function, the premise for the analytical expression ① to be meaningful is $\left|C\sqrt{{\displaystyle \frac{R^2-r^2}{r^2}}}\right|\le 1$

That is

$$r\ge {\displaystyle \frac{CR}{\sqrt{C^2+1}}}$$

2.2. Python Visualization Code

2.2.1. Visual Code Design Ideas

(1) Function thetar

This function is used to calculate a certain θ under a certain r and the angle value θ ₀ of the direction where the particle is located at the initial moment according to the analytical formula of the fastest curve ①. The input parameters are r, R, C and θ₀ are calculated and output in two cases:

①At $C\sqrt{{\displaystyle \frac{R^2-r^2}{r^2}}}\to 1$ that time, $\theta ={\theta }_0-{\displaystyle \frac{\pi }{2}}\left(1-{\displaystyle \frac{C}{\sqrt{C^2+1}}}\right)$

② Otherwise, still according to the calculation

$$\theta ={\theta }_0-\arcsin \left(C\sqrt{{\displaystyle \frac{R^2-r^2}{r^2}}}\right)+{\displaystyle \frac{C}{\sqrt{C^2+1}}}\arctan \left[\sqrt{C^2+1}{\displaystyle \frac{\sqrt{\frac{R^2-r^2}{r^2}}}{\sqrt{1-C^2\left(\frac{R^2-r^2}{r^2}\right)}}}\right]$$

(2) Function plot_q

This function is used to draw the fastest curve when C is a specific value. The thetar function needs to be called to calculate the coordinates of each point required to draw the graph.

(3) Main program

Set R (for the convenience of drawing, set R to 1) and the angle value θ ₀ of the direction where the particle is located at the initial moment, and call the function plot_q to draw the fastest curve under different C values.

For each C value, the plot_q function is called to draw the corresponding curve and set the legend. Finally, plt.show() is used to display a series of images of the fastest curves under different C values.

2.2.2. Visualization Code and Visualization Results

According to the design ideas described in Section 2.2.1, the visualization code is shown above (which is also shown in Figure 2), and the visualization result is shown in Figure 3.

import matplotlib.pyplot as plt

import math

def thetar(r,R,C,theta_0):

 if abs(C*math.sqrt(R**2/r**2-1)-1)<1.0e-15:

  return theta_0-math.pi/2*(1-C/math.sqrt(C**2+1))

 else:

  return theta_0-math.asin(C*math.sqrt(R**2/r**2-1))+C/math.sqrt(C**2+1)*math.atan(math.sqrt(C**2+1)/C*math.tan(math.asin(C*math.sqrt(R**2/r**2-1))))

def plot_q(R,C,theta_0):

 rlist=[];thetalist=[];xlist=[];ylist=[]

 for i in range(1001):

   theta=thetar(R-i*(R-R/math.sqrt(1+1/C**2))/1000,R,C,theta_0)

   rlist.append(R-i*(R-R/math.sqrt(1+1/C**2))/1000)

   thetalist.append(theta)

   xlist.append((R-i*(R-R/math.sqrt(1+1/C**2))/1000)*math.cos(theta))

   ylist.append((R-i*(R-R/math.sqrt(1+1/C**2))/1000)*math.sin(theta))

 for i in range(len(rlist)):

   theta=-thetalist[len(thetalist)-1-i]-math.pi*(1-C/math.sqrt(C**2+1))+2*theta_0

   xlist.append(rlist[len(rlist)-1-i]*math.cos(theta))

   ylist.append(rlist[len(rlist)-1-i]*math.sin(theta))

 plt.gca().set_aspect('equal',adjustable='box');plt.plot(xlist,ylist,label=f'C = {C}')

R=1;theta_0=0

Clist=[0.1,0.2,0.4,0.8,1.6]

for C in Clist:

 plot_q(R,C,theta_0)

 plt.legend(loc='center left',bbox_to_anchor=(1,0.5))

def c(R):

 x=[];y=[]

 for i in range(1001):

   x.append(R*math.cos(i*2*math.pi/1000))

   y.append(R*math.sin(i*2*math.pi/1000))

 return x,y

cx,cy=c(R)

plt.plot(cx,cy)

plt.show()

Figure 2. Python Visualization Code.

Figure 2. Python Visualization Code.

Figure 3. Python Visualization Results.

Figure 3. Python Visualization Results.

3. Python Visualization Result Analysis

As can be seen from Figure 3, the smaller C is, the longer the distance the particle moves inside the earth, and the closer the steepest descent line is to a straight line. When θ₀ is fixed, given an integral constant C, there is only one steepest descent line. Here we can calculate the time required for the particle to move from the initial point to the end point along this steepest descent line:

From Section 2.1, we can see

$$T={\displaystyle {\int }_A^B\frac{\sqrt{{\left(\frac{dr}{d\theta }\right)}^2+r^2}}{\sqrt{k\left(R^2-r^2\right)}}d\theta }$$

Substituting ${\displaystyle \frac{dr}{d\theta }}=-r\sqrt{{\displaystyle \frac{1}{C^2}}{\displaystyle \frac{r^2}{R^2-r^2}}-1}$ into the above formula, we can get

$$T={\displaystyle {\int }_A^B\frac{1}{c\sqrt{k}}\frac{r^2}{R^2-r^2}d\theta }=-{\displaystyle {\int }_A^B\frac{1}{c\sqrt{k}}\frac{r}{R^2-r^2}\frac{dr}{\sqrt{\frac{1}{C^2}\frac{r^2}{R^2-r^2}-1}}}$$

Substituting the upper and lower limits of the integral at points A and B, we get $T={\displaystyle \frac{\pi }{\sqrt{k}}}{\displaystyle \frac{1}{\sqrt{C^2+1}}}$

Using the Python visualization toolkit, write code to plotTagainst C, showing therelationship between the corresponding C values of particles in different tunnel curves and the time T required to travel from point A through the Earth's interior to point B (as shown in Figure 4, the dashed line represents the reference value of 42 minutes in the case of degradation to a straight line). At this point, substitute the real data with a value of 1.24×10^-3 rad/s.

Next, let's consider the actual physical meaning of the integral constant C. Assume that the angular coordinate of the particle leaving the earth tunnel is θ _ex , then from the analytical expression of the fastest curve in 2.1, we know that

$${\theta }_{ex}={\theta }_0-\pi \left(1-{\displaystyle \frac{C}{\sqrt{C^2+1}}}\right)$$

Assume that the difference in angular coordinates of the particle entering and exiting the earth tunnel is θ _ab = θ ₀ - θ _ex ( θ _ab ∈[0, π ]), then the magnitude of the integral constant C is

$$C={\displaystyle \frac{1-\frac{{\theta }_{ab}}{\pi }}{\sqrt{1-{\left(1-\frac{{\theta }_{ab}}{\pi }\right)}^2}}}$$

This formula gives the physical meaning of the integration constant.

Substituting the above formula back $T={\displaystyle \frac{\pi }{\sqrt{k}}}{\displaystyle \frac{1}{\sqrt{C^2+1}}}$ into the equation, we get

$$T={\displaystyle \frac{\pi }{\sqrt{k}}}\sqrt{\left(2-{\displaystyle \frac{{\theta }_{ab}}{\pi }}\right){\displaystyle \frac{{\theta }_{ab}}{\pi }}}\le {\displaystyle \frac{\pi }{\sqrt{k}}}$$

Similarly, toclearly understand the relationship between the angular coordinate difference θ ab and the corresponding dimensionless integration constant C, as well as the required time T, we use Python visualization toolkit,write code to plot T and C against θ_ab, as shown in Figure 5 and Figure 6.

When θ _ab = π, the problem is reduced to the case of digging a straight earth tunnel inside the earth (as shown by the dashed line in Figure 6). From the above equation combined with Figure 4 and Figure 6, it can be seen from the above formula that digging a straight earth tunnel inside the earth does not actually make the time taken by the particle from point A to point B through the earth's interior the shortest, and in general, the straight tunnel is the path that takes the longest time.

4. Conclusion and Expectation

Python's strength lies in its ability to handle large datasets and to integrate with other data analysis tools. However, it requires a certain level of programming knowledge and can be overkill for simple visualization tasks,while students and teachers possibly do not have the knowledge that supports the visualization process in Python.Also,visualization results generated using Python can sometimes appear less smooth when the number of manually controlled sampling points is insufficient. This limitation arises because the quality and smoothness of a visualization are highly dependent on the density and distribution of the data points. With fewer sampling points, the visual representation may not accurately capture the underlying patterns or trends in the data. Additionally, the lack of sufficient sampling points can lead to jagged or discontinuous lines in plots, which can be misleading or difficult to interpret,which might confuse the thinking of students. Furthermore, this issue can also affect the ability to interact with the visualization in real-time. Real-time interaction typically requires a dynamic and responsive system that can update the visualization as the user manipulates the data or explores different aspects of the visualization. When the data is sparse or the sampling points are not optimally distributed, the system may struggle to provide a seamless and responsive user experience. This can result in delays or inaccuracies in the visualization updates, which can hinder the user's ability to explore the data effectively and make informed decisions based on the visual insights.

All in all,in this paper, we explored the problem of the fastest motion of an object in a curved earth tunnel by building a mathematical model and using the Python visualization toolkit. Through in-depth analysis of the model, we not only derived the analytical expression of the brachistocentric curve, but also revealed the physical meaning of the integral constant C, which is related to the difference in the angular coordinates of the particle entering and exiting the earth tunnel. Our results show that, contrary to intuition, a straight tunnel is not always the fastest path, and a curved tunnel can significantly reduce the time it takes for an object to traverse the earth in most cases. In addition, this research method that combines mathematical software and physical models can help students better expand their thinking about physical models. With the assistance of software tools, even complex extended physical problems can become easy to understand and master, which is of great significance for improving teaching effectiveness and students' learning efficiency.