A simple intuitive method for seeking intersections of hyperbolas for acoustic positioning biotelemetry

We proposed a simple hyperbolic positioning method that does not require solving simultaneous quadratic equations. Moreover, we introduced the mathematical concept of a “pencil” into analytical calculations in the hyperbolic positioning method for a better understanding. In many recent studies using positioning biotelemetry, the specific procedure for intersection calculation of hyperbolas has rarely been described. This might be one of two major obstacles, with the other being clock synchronisation among receivers, for positioning biotelemetry users, including potential users. We focus only on the intersection calculation in this paper. Therefore, we propose a novel method and introduce the mathematical concept into analytical calculations. The computing performances of the novel method, an analytical method applying the concept of a pencil, and an approximating method using the Newton-Raphson method were compared regarding positioning correctness, accuracy, and calculation speed. In the novel method, hyperbolas were represented using the parameter θ, which was treated as a discrete variant. The finer the tick-width of the parameter θ, the more accurate was its positioning, but it took slightly longer to calculate. By setting the tick-width to 0.01°, a simulated trajectory was correctly and accurately localised, as in the analytical method which always correctly returned the accurate solution. The approximating method has a major limitation concerning correctness. It returns a single solution regardless of two intersections of hyperbolas; however, the positioning is accurate when the hyperbolas intersect at a single point. This study approached one major difficulty in positioning biotelemetry and will help biotelemetry users overcome this drawback with a simple and intuitive understanding of hyperbolic positioning.


Introduction
Acoustic biotelemetry techniques have been applied to elucidate the ecology and biology of various free-ranging aquatic animals [1][2][3][4]. Biotelemetry systems comprising ultrasonic tangential line of a target function. It is a powerful approximation algorithm, especially with regard to the convergent speed if an appropriate initial value is set. This method is implemented through several programming languages such as the nleqslv function in R [19], which makes writing a code for it relatively easy. However, the Newton-Raphson method will not converge if an inappropriate initial value is set, and it does not always return correct solutions if there are multiple roots. There are many methods to find intersections of hyperbolas, and all of them can potentially provide the same answer. However, the specific procedure for intersection calculation using acoustic positioning biotelemetry has rarely been described previously. Recent papers have vaguely indicated a calculation method without citation [8,9,[11][12][13][14]20,21]. Calculating the intersection of hyperbolas is a major difficulty that biotelemetry users encounter while obtaining positioning results from acquired receiver data. The analytical (algebraic) way to calculate intersections between two quadric surfaces such as hyperboloids and ellipsoids in a context of the hyperbolic navigation have been described [15,16]. It is mathematically more general than intersection calculation between quadratic curves. However, it is not intuitive because they are not examples of biotelemetry and involves a redundant list of equations. Other papers have illustrated calculation methods in the context of biotelemetry [17,18]. These have been described with only the list of mathematically technical equations, thus it would be essentially the same as the previous set of papers with respect to difficulty associated with intuitive understanding [17]. Furthermore, a calculation method can be used only on a 2-D plane, not in a 3-D space, and has an ambiguous choice of solutions [18]. When solving numerically (approximately) rather than analytically, it may result in the problem of always returning a single answer, even though there are two possible answers, and the calculation method used for this is not clearly stated [8,9,[11][12][13][14]20,21]. Another obstacle in the positioning procedure, is clock synchronisation among the receivers. However, in this study, we focus only on intersection calculation and not on clock synchronisation because the latter has been well described [14,22]. Therefore, to overcome this major obstacle faced during positioning by biotelemetry users, we aimed (1) to propose a simple and intuitive method for obtaining intersections in a hyperbolic positioning method without solving a simultaneous quadratic equations, (2) to introduce a mathematical concept into analytical (algebraic) calculations for better, consistent understanding; and (3) to compare the positioning performance of this novel method, an analytical (algebraic) method using the introduced concept, and a numerical (approximating) calculation method.

Algorithm of proposed positioning method
An outline of the basic concept of a hyperbolic positioning method has been provided previously. Geometrically, an intersection of a set of hyperbolic LOPs is a transmitter position on a 2-D plane, while in an underwater 3-D space, the intersection of a set of hyperboloids, especially one sheet of a circular hyperboloid of two sheets, and a plane defined by a transmitter depth is the transmitter position. Transmitter depth is obtained if a signal includes depth information (a transmitter has a depth sensor); otherwise, depth should be assumed appropriately, in which case, it might be the same depth as the receivers. Note that we assume that all receivers were deployed at the same depth. One can obtain a hyperbola in the 2-D plane when slicing a hyperboloid by a depth plane. Therefore, we calculate the intersection of a set of hyperbolic LOPs.
The hyperbola is defined by a locus of points P(x, y) such that the difference between the distances from two fixed points F 1 (f, 0), F 2 (-f, 0), f > 0 (the foci) is a constant, denoted by 2a, Eq (1) leads to an equation of the hyperbola as follows: where a > b > 0 and f > 0. While considering acoustic positioning in biotelemetry, let foci F 1 and F 2 be acoustic receivers R 1 and R 2 , respectively. Assume that an acoustic signal from a particular transmitter at point P(x, y) is detected by the two receivers at t 1 and t 2 . Then, the time difference of arrival (TDOA) is |t 2 -t 1 |. Let c be the speed of underwater sound. Let 2a be the difference in the distance between the two receivers and the transmitter.
Eq (4) can be used to define the hyperbola. Therefore, the transmitter will be on a hyperbola which satisfies Eq (4). Parameters a and b are defined by transforming Eqs (3) and (4). Let d be the distance between R 1 and R 2 , thus, ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi d 2 A hyperbola can be also represented by introducing a new parameter θ (-π/2 < θ < π/2, π/2 < θ < 3π/2) as follows (Fig 1): Accordingly, a and b represent the coordinates of a hyperbola with θ.
Here, considering the order of arrival times t 1 and t 2 , and the domain of θ, only one branch (R1-or R 2 -side branch) of a hyperbola, that is, a hyperbolic LOP, can be defined. In the case of t 2 -t 1 > 0, in which R 1 detected the signal primarily rather than R 2 , meaning that the transmitter would exist closer to R 1 than to R 2 , so R 1 -side branch should be selected. In the other case of t 2 -t 1 < 0, R 2 -side branch can be defined in the same manner. (Note that no hyperbola exists in the case of t 2 -t 1 = 0, i.e., a vertical bisector of R 1 R 2 .) This fact is equivalent to considering the sign of Eq (7). Assuming that the domain of θ is -π/2 < θ < π/2 such that cosθ > 0, the sign of a prescribes the sign of x. Therefore, by letting the sign of t 2 -t 1 prescribe the sign of a, we redefine a as a signed value as follows: If transmitter depth is known, a hyperbola can be obtained by slicing a circular hyperboloid of two sheets in a plane of transmitter depth (Fig 2). Parameters a' and b' in this case are defined by multiplying a correction coefficient with a and b, respectively. a 0 ¼ a ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where z 0 is the absolute difference between receiver installation depth and transmitter depth. Correction coefficient is derived from the equation of a circular hyperboloid of two sheets. S1 File shows the detailed procedure for derivation. The hyperbola is drawn on a coordinate system in which the receivers R 1 and R 2 are foci on the x-axis (Fig 3A). By rotating and translating the foci to match the actual receiver's position in real coordinates, we can obtain a hyperbola on the real coordinates (Fig 3B and 3C). The translation distances δx and δy are distances from the origin to the midpoint between actual receiver positions of R 1 (x R1 , y R1 ) and R 2 (x R2 , y R2 ).

PLOS ONE
The rotation angle φ is an angle formed by x-axis and R 1 .
Note that φ is easily calculated using arctan2 or atan2, which is a powerful programming function implemented in almost all programming languages that returns an exact and unique value in the range of -π < φ � π. For x and y calculated from Eqs (7) and (8), the coordinates X and Y given a rotation and translation are as follows: Y ¼ x sin φ þ y cos φ þ dy: ð16Þ ( As described above, we can obtain the hyperbolic LOP H 1 from the TDOA between R 1 and R 2 , and those coordinates ( Fig 3C). The same procedure for receivers R 1 and R 3 allows us to draw another hyperbolic LOP H 2 . The intersection of H 1 and H 2 is the position at which the transmitter emitted the signal, that is, the estimated position ( Fig 4A).
We shall present our calculation method for the intersection of H 1 and H 2 . Let H n -be the coordinate system where the hyperbola before rotation and translation are given. Let the R-be the coordinate system where the hyperbola after rotation and translation are given. Let (x n , y n ) and (X n , Y n ) be the coordinates of the hyperbola on H n -and R-coordinate systems, respectively. (X n , Y n ) is a function of (x n , y n ), rotation angle φ n , translation distance δx n and δy n , and parameter θ. Here, we calculate the intersection of H 1 and H 2 on H 1 -coordinate system by rotating and translating H 2 (X 2 , Y 2 ). The new coordinates of H 2 '(x 2 ', y 2 '), which are translated and rotated to fit H 1 -coordinate system, are as follows:

PLOS ONE
Solving Eq (2) for x, we can consider x as a function of y, meaning that we can represent the coordinate x 1 of H 1 as a function of given y on H 1 -coordinate system. Thus, x 1 is defined by y 2 ' that is y coordinate of H 2 on H 1 -coordinate system.
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where a 1 and b 1 are derived from Eqs (9) and (6) (or (10) and (11)), respectively. When H 1 and H 2 are placed on H 1 -coordinate system, there is a unique pair of x 1 and x 2 ' for any y 2 '. Δ is defined as the absolute distance between x 1 and x 2 ': Δ represents x-axis distance between H 1 and H 2 for any y 2 ' on H 1 -coordinate. Therefore, it is clear that H 1 and H 2 must intersect at the point where Δ = 0 ( Fig 4B and 4C). Further, everything from hyperbolic LOP (x n , y n ) to Δ is consistently a function of θ. If, θ is treated as a discrete sequence, there is a consistent correspondence between θ and Δ via index i from θ to Δ ( Fig 5). Hence, seeking an index j that suffices Δ(j) = 0, or Δ(j) � 0, is equal to seeking the intersection of H 1 and H 2 . Thus, the intersection on the R-coordinate system, which we want to find, is (X 2 (j), Y 2 (j)).

Analytical solution
Simultaneous quadratic equations must be solved to analytically determine the intersection of hyperbolas. As it is a troublesome task, we introduce a mathematical concept of a 'pencil' to help understand this calculation procedure. When the pencil is applied to the intersection calculation of two quadratic curves, finding intersections of the two curves is transformed into finding intersections of one of the curves with straight lines, which are generated from the two curves. Let both f(x, y) = 0 and g(x, y) = 0 be quadratic curves, and assuming that those curves intersect at a point of (s, t). Then, f ðs; tÞ ¼ 0^gðs; tÞ ¼ 0 ð21Þ where λ is any real number. This means Eq (22) always passes through (s, t) regardless of the value of λ. Therefore, f(x, y) + λg(x, y) = 0 is a set of quadratic curves including a line passing through (s, t), which is an intersection of f(x, y) = 0 and g(x, y) = 0 ( Fig 6A). This is a brief description of using the concept of pencil.
The specific procedure for intersection calculation of two hyperbolas in positioning biotelemetry is as follows. An acoustic signal emitted by a transmitter is detected by three fixed receivers R a (a x , a y , 0), R b (b x , b y , 0), and R c (c x , c y , 0) at time t a , t b and t c , respectively (Fig 7A). Note that the receivers are assumed to be deployed at the same depth. The TDOA between R a and R b , and R a and R c are T ab = t a -t b , and T ac = t a -t c , respectively. Let c be the underwater speed of the sound. The position of a transmitter (x, y, z) suffices following equations: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À a x Þ 2 þ ðy À a y Þ 2 þ z 2 q À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À a x Þ 2 þ ðy À a y Þ 2 þ z 2 q À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where, These equations provide the definition of a hyperbola. Transposing the first term of the left-hand side to the right-hand side and squaring the two sides, we obtain 2ðb x À a x Þx þ 2ðb y À a y Þy À 2R ab ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À a x Þ 2 þ ðy À a y Þ 2 þ z 2 2ðc x À a x Þx þ 2ðc y À a y Þy À 2R ac ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðx À a x Þ 2 þ ðy À a y Þ 2 þ z 2

PLOS ONE
Substituting this λ into f(x, y) + λg(x, y) = 0 and solving it, we obtain the following equation for a straight line. where, 2ðR ac ðb y À a y Þ À R ab ðc y À a y ÞÞ : The intersection of Eqs (25) and (26) is equivalent to the intersection of Eqs (25) and (28) (Fig 7B). Substituting Eqs (28) into (25) and solving it, we obtain If the depth of transmitter is known as z 0 , that is, the transmitter has a depth sensor, substituting z = z 0 with Eq (29); otherwise, use an appropriate value, which might be zero in most cases, to z 0 . If the installation depth of the receivers is not zero, which is a natural situation, z 0 is the absolute difference between the installation depth and the transmitter depth. Applying the quadratic formula to Eq (29) to find x, we obtain ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi : 30 ð Þ Then, substituting x into Eq (28) to find y, two candidates for the intersection are obtained. If there is only one intersection, one of the two candidates should be selected after considering the condition. For example, substituting x and y in Eq (23), the pair of x and y that holds equality is the coordinate of the intersection. Where there are two intersections, both pair of x and y pairs hold equality. In the other case where there is no intersection, x and y become imaginary numbers.

Comparison of positioning results
We compared the positioning results of a simulated trajectory using the proposed method, analytical method, and approximation method. A simulated trajectory comprising 200 points was constructed in a 3-D space as it passed through an array of three receivers. The array formed an equilateral triangle of 50 m on each side. Installation depths of the three receivers were set at the same depth of 15 m. It was assumed that the signal, including depth information from a transmitter on the trajectory propagated at the underwater sound speed of 1500 m/s, and that it was detected by all three receivers. There were no misdetections or multipath effects; thus, all data would be correctly calculated. The comparison items are as follows: (1) the number of false answers, i.e., whether the number of solutions is correctly returned, (2) positioning accuracy, and (3) computation time. For (1), the analytical method is treated as it always returns true solutions. For (2), accuracy was the average distance from the true position and the estimated position only when the hyperbolas intersected at one point. For (3), it was expressed as a ratio to the average of 200 calculations using the analytical method, to clearly show differences between the methods, as actual calculation period may highly depend upon the machine power of a computer used. The proposed method was set to three types, with θ tick width of 1˚, 0.1˚, and 0.01˚. The analytical method was established by applying the concept of a pencil. An approximation method was constructed using the Newton-Raphson method. The initial coordinates were set at the centroid of the receivers. Iterations were terminated when the difference between the x-and y-coordinates in the update was <0.001, or when the number of iterations was >5. All the procedures of the three methods were coded in R ver. 4.0.5 [23]. The source codes of the three methods and coordinates of the trajectory are shown in S2 and S3 Files. Computing was performed on a laptop computer of MacBook Pro (13-inch, 2017, Four Thunderbolt 3 Ports) with a processor of 3.1 GHz and 8 GB random-access memory.

Results
The five computing results of the three methods are graphically presented in Fig 8. The analytical method had no false answer with positioning accuracy of zero ( Fig 8E and Table 1). Similar to the analytical method, the proposed method with θ tick width of 0.01˚also had no false answer with approximately zero accuracy ( Fig 8A and Table 1). The proposed method with 0.1˚had only one false answer ( Fig 8B and Table 1), and that with 1˚θ tick width had ten false answers including three returns with no answer where there was a single intersection (Fig 8C  and Table 1). Positioning accuracy of the proposed method got lower as θ tick width got larger ( Table 1). The approximating method had 29 false answers, meaning this method always returned single answer even though there were two intersections (Fig 8E). However, positioning accuracy of this method was zero in the case that there was only one intersection (Table 1). Computing period of the analytical method was 3.81 × 10 −4 sec in average (n = 200). In comparison, the average computing period of other methods ranged for 0.89-9.05 times ( Table 1).

Discussion
We proposed a simple, intuitive method and introduced a mathematical concept for intersection calculation in hyperbolic positioning. The calculation procedure in acoustic positioning biotelemetry would have two difficulties in obtaining the positions of target animals: clock synchronisation among receivers and intersection calculation between hyperbolas. A brief explanation of the clock synchronisation technique has been described [14,22], so we focused on intersection calculation in this study. The clock can be synchronised by linear regression using detection data recorded in receivers [14,22]. However, it is slightly more difficult to obtain positions by intersection calculation because of the need to solve simultaneous quadratic equations. This task for finding intersections could be a great obstacle for calculating positioning by biotelemetry users, including potential users. Therefore, we focused on intersection calculation methods.

PLOS ONE
From the computing results for the positioning of the simulated trajectory, it was found that the proposed method with θ tick width of 0.01˚was able to localise the trajectory with an accuracy of zero m without any false answers, including the analytical method (Fig 8A-8E). The proposed method with θ tick width of 1˚had several false answers and noticeable errors (Fig 8C). The method with θ tick width of 0.1˚had a higher accuracy, but one false answer was given (Fig 8D). When the returned answer was incorrect, it was because the shape of a hyperbola could only be roughly drawn owing to the resolution of θ. If there were two intersections and the method returned only one solution, the two intersections were too close together to be distinguished. If there was a unique solution but no solution was returned, that possible solution was eliminated because a local minimum of Δ was larger than the threshold that was set to omit false intersections (Fig 9B). There was no false answer with θ tick width of 0.01˚; however, computing period increased considerably as the tick width became finer, owing to the larger length of the variables handled. Although the computing period should be short, it took an average of nine times longer than 3.81 × 10 −4 s, which was on the order of 10 −3 s. Although the period depends largely on machine specifications of the computer used for the calculations, the computer used in this study was a laptop, which is commonly available on the market. Therefore, the computing period is unlikely to be a major bottleneck. For the approximating method, there was a serious limitation regarding false answers, that is, the method always returned a single solution even if there were two intersections of hyperbolas. This is a weakness of approximating methods, including the Newton-Raphson method. The approximating method is relatively easy to code using existing functions (e.g., the nleqslv function in R [19]); therefore, special care should be taken when using it. Note that the positioning error is almost zero within <5 iterations (Fig 8B), indicating that the convergence of this method is very fast. This is a strong point of the approximation method using the Newton-Raphson method.
Our proposed method is a little complicated because it employs parameter θ to express coordinates of a hyperbola. However, although the equations may seem slightly complicated as they involve trigonometric functions, the operation is essentially simple since it only requires rotations and translations. We drew enough figures to make it intuitive and easy to understand. The proposed method can simply find intersections during position calculation without

PLOS ONE
the requirement to solve quadratic equations. We employed two ideas for the calculation: parameter θ, signed a. First, parameter θ enables us to cover the possible range of a hyperbola. We can also obtain the x coordinate of a hyperbola using Eq (19) in the H n -coordinate system. However, considering the possible range of a hyperbola by Eq (19), y must range from -inf to inf, which cannot be represented using programming languages. (Note that x 1 is calculated directly from y 2 ' in Eq (19) for calculating Δ, but this does not matter here because y 2 ' has already been finitely determined by parameter θ.) A larger range of y would simply require a larger computational complexity if the tick value is consistent. Although it might be possible to apply an optimal range of y for each case, this method would not be simple. Parameter θ is convenient because we do not have to consider it. Second, the signed a enables us to omit the case classification when defining a hyperbolic LOP. Generally, a is larger than 0 because a is the distance from the origin to the vertex of a hyperbola. In the case of a > 0, we must choose one branch of a hyperbola (e.g., R 1 -or R 2 -side branch) based on the order of arrival times. However, integrating the order of arrival times into a plus-minus sign would simply and uniquely define which branch is a hyperbolic LOP.
The proposed method can detect an index when there are two intersections of hyperbolic LOPs if the tick width of θ is adequate (Fig 8). To find an index j that suffices Δ(j) = 0, we can explore the raw Δ value without using a technique such as differentiating Eq (20) by x 2 ' because Δ comprises an array of discretised values. However, due to the discretisation, there is not always an index j where Δ(j) = 0. In addition, there are theoretically one or two intersection points. For these two reasons, it is convenient to consider the local minimum value instead of the minimum value or zero. Although there can be a case where three indices of local minimal value exist, we can omit the third index that does not detect an intersection by some threshold or the order of local minimal values, because intersections must theoretically be up to two points (Fig 9). Practically, the existence of two intersections means that it is impossible to provide a unique position unless a unique position is calculated from another combination of three receivers in the array.
We introduced the mathematical concept of a pencil to help with understanding the calculation procedure intuitively and graphically for hyperbolic positioning. The listed calculation procedure using the equations is approximately similar to when simultaneous equations of two hyperbolic LOPs are solved [15]. However, rather than simply listing the preceding equations, we have tried to present them in a way that makes what is being calculated here, explicit, by employing a pencil. We have also attempted to introduce a pencil so that it can be applied in the future. The concept of a pencil is a powerful tool for finding a sound source for acoustic positioning biotelemetry. Although it is assumed that installation depth of the three receivers is the same in the intersection calculation described in this study, this technique can be applied to different positioning cases. Such cases include, 3-D positioning using four receivers without depth information of a transmitter. The calculation procedure in each case is highly redundant, therefore all of it has not been discussed here; it has been described in S4 File. In positioning biotelemetry, fine-scale positioning of <1 m order will be one trend to observe inter-individual or inter-specific interactions between animals (e.g. [12,24]). Positioning results are geometrically dependent upon coordinates of the receivers, including installation depth, and are particularly sensitive in the vicinity of the receiver because of the shape of a hyperbola or hyperboloid. Mathematically rigorous methods of intersection calculations under various conditions are important in this area.

Conclusions
This study aims to achieve two goals: proposing a novel positioning method which does not require solving simultaneous quadratic equations, and introducing a mathematical concept of a 'pencil' to an analytical (algebraic) positioning method. First, it was found that our proposed method has potential to become a new method in positioning biotelemetry, so far as the tick width remains adequately fine. From the results of positioning simulation, our proposed positioning method could accurately obtain the intersection between hyperbolic LOPs with no false return similar to the analytical method, which is mathematically correct, if tick width of parameter θ was �0.01˚. In contrast, Newton-Raphson method, and our proposed method without fine tick width found some false answers that come from approximating methods, especially in Newton-Raphson method. Second, it was also proved through positioning simulation that an analytical method could be constructed by introducing the concept of a 'pencil,' and it could correctly calculate an intersection in positioning biotelemetry. We proved that the analytical method constructed using the concept of a pencil is highly applicable, and that it can be applied to the more complex case of positioning biotelemetry such as calculating 3-D positioning using four receivers as presented in S4 File. However, approximately identical results will be obtained by all methods i.e., the analytical method, proposed method, and approximating method used in this study, although this scenario is only in the case that there is single intersection (i.e., excluding false answers caused by approximating way), because all methods are constructed based on mathematics. Nevertheless, seeking the intersection of hyperbolas may be a major obstacle for calculating positioning by biotelemetry users, including potential users. We tried to overcome this obstacle by enabling the users to solve this issue through a simple and intuitive method that we proposed and introduced in this paper, with a better understanding, rather than through troublesome tasks.