Effectiveness of the Chebyshev Approximation in Magnetic Field Line Tracking

: The tracking of magnetic ﬁeld lines can be very expensive, in terms of computational burden, when the ﬁeld sources are numerous and have complex geometries, especially when accuracy is a priority, because an evaluation of the ﬁeld is required in many situations. In some important applications, the computational cost can be signiﬁcantly reduced by using a suitable approximation of the ﬁeld in the integrated regions. This paper shows how Chebyshev polynomials are well-suited for ﬁeld interpolation in magnetic ﬁeld-line tracking, then discusses the conditions in which they are most appropriate, and quantiﬁes the effectiveness of parallel computing in the approximation procedures.


Introduction
The tracking of magnetic field lines (FL) is required in many scientific and industrial applications. Line-tracking can be carried out by integrating an ordinary differential equations (ODEs) system. A thorough survey of the methods for the numerical integration of a differential equation that preserves geometric properties, such as the volume for solenoidal fields, can be found in [1]. Such methods, including that of Runge-Kutta, can be used for tracing the FLs of the magnetic field [2,3], and can find employment in several scientific applications, including fusion devices based on magnetic confinement [3][4][5][6].
Unfortunately, the FL tracking procedures can be very computationally expensive if the sources are complex since the continuous recalculation of the field along the line trajectory is required. The computational cost also increases with accuracy, which implies more detailed source modeling and a reduced integration step. The effectiveness of parallel computation is rather limited in the tracking formulations since the points where the field has to be evaluated cannot be known a priori because they depend on the particular trajectory under construction.
In this class of applications, the use of interpolation techniques [7] based on the projection of the field in a suitable finite-dimensional functional space can be very beneficial. The best choice for the dimension of the space and for the class of interpolating functions depends on the desired accuracy and on the size of the interpolation domain, which also impacts on the number of direct computations of the "actual" field required to finalize the interpolating functions.
In general, the approximation is convenient if, for a given accuracy, the interpolating procedure is faster than the direct field calculation and, in addition, the number of the points where the field is required is larger than the number of direct calculations required by the interpolation. The aim of the paper is focused on the methodological aspects, but to assess points where the field is required is larger than the number of direct calculations r by the interpolation. The aim of the paper is focused on the methodological asp to assess the proposed approach, the paper also includes a number of demon applications, particularly in the field of fusion technology. This paper aims to show, from the methodological point of view, how Ch polynomials [7] (CP) are well-suited for field interpolation in FL tracking proc discuss the conditions under which they are convenient, and to show how parallel computing can be in the approximation procedures.
The paper is organized as follows. Section 2 gives a brief description of the m system in a Tokamak device. Section 3 discusses the analytical and numerical as the approach. Finally, in Section 4, the accuracy of the approach and its perform FL tracking are evaluated, and some of its applications are examined.

The DTT Tokamak
The applications considered in this paper are related to the analysis of certain aspects of the magnetic field inside a Tokamak device [8]. The Tokamak is t promising scheme for future commercial reactors devoted to the production o electricity. It includes a very complex magnetic-field system that is designed to and stabilize the plasma in a doughnut-shaped container. Due to its toroidal cylindrical coordinate system [R, φ, Z] is commonly used, with ∈ 0, 2 .
Here, the DTT (Divertor Tokamak test) facility [9] is considered. The DTT Tokamak device under construction on the ENEA site in Frascati, Italy, designed power and particle exhaust and the possible divertor concepts. The main parame major radius 2.19 m, minor radius 0.70 m, pulse length 100 s, plasma current toroidal field 6 T, and heating power 45 MW. Figure 1 shows a rendering of the D assembly, and testing area.   In the DTT device, both the PF and CS systems are up-down symmetrical. Therefore, Table 1 shows only the main parameters of the upper-side coils, including the coordinates of the coil centers, the dimensions of the cross-sections, and the number of turns, in both radial and vertical directions.   In the DTT device, both the PF and CS systems are up-down symmetrical. Therefore, Table 1 shows only the main parameters of the upper-side coils, including the coordinates of the coil centers, the dimensions of the cross-sections, and the number of turns, in both radial and vertical directions.  In the DTT device, both the PF and CS systems are up-down symmetrical. Therefore, Table 1 shows only the main parameters of the upper-side coils, including the coordinates of the coil centers, the dimensions of the cross-sections, and the number of turns, in both radial and vertical directions. The finite number of TF coils breaks the axial symmetry of the device, introducing a periodical deformation of the magnetic field (ripple) in the toroidal direction that can be deleterious to the plasma. Different solutions to counteract this effect are proposed, such as ferromagnetic inserts [10]. More details about the model used to represent the field of the TF coil system are given in Section 4.

FL Tracking
By definition, the FL is the line that is tangential to the magnetic field everywhere. Therefore, line tracking can be carried out by integrating the following ODEs system: where r is the field point and t is an FL abscissa.
Except for a few trivial cases, no analytical solutions are available; therefore, the previous ODEs system must be solved numerically. Many approaches, based on explicit or implicit schemes and characterized by different accuracy levels, computational costs, and performance, are available in the scientific literature. Here, the Runge-Kutta-4 (RK4) method is adopted, based on an explicit fixed-step scheme. For each step, the method consists of four consecutive field calculations that are used to estimate the field component derivatives and accumulates an error in the order of O h RK 4 , where h RK is the step size [11].
It is worth noting that the direct integration of (1) is not enough to force the magnetic flux density to be divergence-free. Alternative approaches are able to preserve the volume in a strong form [3,5,12,13], but, unfortunately, this entails an increasing number of field computations.

The Field Approximation
The field approximation approaches can provide an effective contribution to this class of problems, reducing the required computational effort while preserving the target accuracy. The general philosophy is to select, in a functional space, S I , a suitable interpolation function F I (r) defined in an interpolation domain, V I, and then to calibrate F I (r) by superimposing the fitting with the actual function F A (r) in a suitable number N samp of sampling points in V I .
Several methodologies have been suggested, each being characterized by the space S I adopted for the approximation, the most suitable geometry of the domain V I , and, finally, the number of points N samp needed to achieve the required accuracy in V I [7].
For the applications described herein, the space of CP with a brick-shaped interpolation domain V I is recommended [7]. In a Cartesian coordinates system [x, y, z], a single reference parallelepiped that is centered in the coordinate system's origin, with edges that are parallel to the coordinate axis, and bounded in The general form of CP used to approximate a vectorial field function F A (x, y, z) in the reference brick is: where N px , N py , and N pz are the polynomials' degrees in the x, y, and z directions, respectively,x,ŷ, andẑ are the unit vectors in the coordinate directions, α ijk , β ijk , and γ ijk are the coefficients of each polynomial, and, finally, T i (x), T j (y), and T k (z) are the first CPs of the i-th, j-th, and k-th degrees, respectively, defined as: T n (ξ) = cos(narccosξ) with |ξ| < 1 and n 0 where n is the degree of the n-th polynomial basis function and ξ is one of the coordinates x, y, and z. In order to apply this approach to a general brick with any size and orientation, a bijective map projecting the coordinates of an actual domain point (x b , y b , z b ) into the reference domain point (x r , y r , z r ) can be introduced: To be able to describe the rigid displacements as translations and rotations, as well as homothetic deformations as compression and expansion, the map has to be an affine transformation of R 3 into R 3 .
Here, for sake of clarity, the CP-based approximation of the flux density is reported only for a homothetic transformation of the reference domain: where the metric coefficients c kh are easily derived by (4), and U i−1 , U j−1 , and U k−1 ; the second-kind CPs of (i − 1)-th, (j − 1)-th and (k − 1)-th degrees, respectively, are defined as: It is worth recalling that each second-kind CP is related to the first-kind CP by a simple derivative relationship: In case the transformation also includes a rigid displacement, the generalized form of (5) can be derived using (4), with simple mathematical operations.
Actually, to impose its flux-preserving property, instead of the flux density B, a magnetic vector potential A (with a Coulomb gauge) is used for the approximation, while the symbols in (2) are used for its formulation. Then, the flux density approximation can be derived as B = rot A.
To calculate the coefficients α ijk , β ijk , γ ijk in (5), the Chebyshev approximation is constrained, in a number N samp of sampling points, to fit the field calculated by a direct procedure. Actually, together with (5), an additional constraint for the gauge on the vector potential is also imposed. In such a way, at each sampling point, three linear equations are imposed, one for each of the three components of the field, with a fourth equation for the gauge on the magnetic vector potential. The number of overall unknown coefficients: α ijk , β ijk , γ ijk is N unk = 3 × (N px + 1) × (N py + 1) × (N pz + 1), and their values can be found by solving the system.
The degrees of the polynomials N px , N py , and N pz , and the number of points N samp must be carefully chosen, balancing the need for accuracy and the computing burden, and taking into account the regularity of the field map. To look for a stable solution and, in addition, to counteract the impact of uncertainties, an overdetermined system is recommended. Therefore, the number 4 × N samp of equations should be suitably larger then N unk .
In the applications described herein, the sampling points are chosen to sample the interpolation domain uniformly, and the linear overdetermined system is addressed with a standard SVD technique.

Results
The main aim of this paper is to propose the use of the CP approximation for the magnetic field computation. The analysis of specific applications is beyond the scope of this paper.
However, in order to assess the performance of the approach, some realistic problems related to the magnetic system of the DTT Tokamak are considered. The main aspects discussed in this section deal with the accuracy of the methodology (Section 4.1), the performance in terms of computing burden (Section 4.2) and FL tracking (Section 4.3), and, finally, the capability to take advantage of general periodicity (Section 4.4).
In this paragraph, depending on the specific application, both the cylindrical [R, ϕ, Z] and Cartesian [x, y, z] coordinate systems are adopted. The two coordinate systems have the same center, and, at the ϕ = 0 section, the correspondence R ≡ x and Z ≡ z stand.
Most of the numerical results discussed herein are carried out with the MISTIC code [14][15][16][17][18], based on the discretization of the sources in interconnected elementary sources, and are able to be run with parallel computing systems [19].

Accuracy Assessment: Magnetic Field Reconstruction Inside a Domain
The accuracy of the CP approximation can be very high. In order to assess the actual level, the field produced by the DTT TF coil system has been considered. The system includes 18 fully superconducting coils that are designed to guarantee 6 T at R = 2.11 m. Each coil includes 80 series-connected filamentary currents, each carrying 44 kA.
In Figure 4, the filamentary current distribution in the cross-section of each coil is sketched, as described by two local Cartesian coordinates, χ and η.
Energies 2022, 15, x FOR PEER REVIEW 6 of 14 discussed in this section deal with the accuracy of the methodology (Section 4.1), the performance in terms of computing burden (Section 4.2) and FL tracking (Section 4.3), and, finally, the capability to take advantage of general periodicity (Section 4.4).
In this paragraph, depending on the specific application, both the cylindrical , , and Cartesian , , coordinate systems are adopted. The two coordinate systems have the same center, and, at the = 0 section, the correspondence ≡ and ≡ stand.
Most of the numerical results discussed herein are carried out with the MISTIC code [14][15][16][17][18], based on the discretization of the sources in interconnected elementary sources, and are able to be run with parallel computing systems [19].

Accuracy Assessment: Magnetic Field Reconstruction Inside a Domain
The accuracy of the CP approximation can be very high. In order to assess the actual level, the field produced by the DTT TF coil system has been considered. The system includes 18 fully superconducting coils that are designed to guarantee 6 T at = 2.11 m. Each coil includes 80 series-connected filamentary currents, each carrying 44 kA.
In Figure 4, the filamentary current distribution in the cross-section of each coil is sketched, as described by two local Cartesian coordinates, χ and η. The magnetic field can be evaluated using the MISTIC code, parallelized on a specific GPU architecture and based on the modeling of each filamentary current in a set of connected current segments. The number of segments used to discretize a single filamentary current (3000) has been chosen in such a way as to guarantee an accuracy of within 1 μT.
In order to quantify the approximation error induced by the Chebyshev interpolation, a simple cubic domain with a side length of 5 cm, centered in [ , , ] = [1.80, 0, 0] m, and with normal edges parallel to the respective coordinate's axis, has been assumed as the testbed.
Each cube-side is sampled in , uniformly distributed, points to be used for CP fitting, so that = . Consequently, the entire domain could be regarded as consisting of (  The magnetic field can be evaluated using the MISTIC code, parallelized on a specific GPU architecture and based on the modeling of each filamentary current in a set of connected current segments. The number of segments used to discretize a single filamentary current (3000) has been chosen in such a way as to guarantee an accuracy of within 1 µT.
In order to quantify the approximation error induced by the Chebyshev interpolation, a simple cubic domain with a side length of 5 cm, centered in [x, y, z] = [1.80, 0, 0] m, and with normal edges parallel to the respective coordinate's axis, has been assumed as the testbed.
Each cube-side is sampled in N side , uniformly distributed, points to be used for CP fitting, so that N samp = N side 3 . Consequently, the entire domain could be regarded as consisting of (N side − 1) 3 elementary domains, m i , with i = 1, 2, . . . (N side − 1) 3 . The same CPs degree, defined in the three coordinate directions for all the CPs, has been chosen (N px = N py = N pz = N pol ). Therefore, the system to be solved consists of 4 × N side 3 linear equations with 3 × N pol + 1 3 unknowns. To estimate the accuracy of the approximation, the relative error has been defined as where (x b i , y b i , z b i ), i= 1, 2, . . . , (N side − 1) 3 , are the barycenter's coordinates of m i , B is the approximated field, B ref is the field obtained by direct evaluation, and 2 is the Euclidean norm.
A parametric analysis has been carried out in terms of the N side and N pol parameters. To guarantee the overdetermination of the system to be solved, N side ≥ N pol + 1 must be imposed. Table 2 shows that the accuracy strongly depends on N pol and slowly decreases with N side . Then, N side = N pol + 1 can be recommended to save the computing burden for any level of accuracy. The approximation error (see Table 2) strongly depends on the sources, the brick's shape, position, and orientation, but in any case, degree levels to the order of a few units are able to provide high accuracy. It should be noticed that the CP approximation method does not depend on the tool used to evaluate the "direct" field. Therefore, the CP interpolation can be applied for any field-calculation procedure [20].

Performance of the CP Approach in Magnetic FL Tracking
The CP approximation can effectively be used for FL tracking, but the efficiency depends on several elements, including the required accuracy and the characteristics of the magnetic field.
The RK4 method, here adopted for the FL tracking, calls for 4 field evaluations for each of the NP line integration steps.
If the CP approach is adopted, the number of direct field evaluations is limited to the total number of sampling points. Therefore, assuming that N brick bricks are required to track the line and that the same number N samp of samples is used in each brick, the total number of direct calculations is N brick × N samp . As a consequence, a gain g CP provided by the CP approach can be defined as the ratio between the direct field calculations, both without and with the CP approximation: Here, an analysis has been made to evaluate the performance of the CP approach in terms of the computational burden. The TF coil system source has also been considered. Starting from the initial point, where x = 2.8 m and y = z = 0, the line is tracked for 20 degrees in a toroidal direction, covering about 1 m of length, and with an integration step h RK = 1 mm, so that NP line is about 1000. For the CP approach, 11 cubic bricks of 10 cm side length are used, while the same number of sample points N side for each brick side is chosen as well.
Without the use of parallelization, the N samp = N side 3 points are serially evaluated in each brick, one brick at a time.
Assuming a CP degree of N pol = 2 (then, N side = 3) the gain is larger than 13 (see Table 3). The CP gain strongly depends on N side because of its impact on (8) in Table 3. It is worth noting that the use of the CP approach may not be convenient for high values of N side (N side > 7 in the case at hand). Table 3 also reports the computing time (using a processor with 2 × AMD Epyc 32-Core 7452 2.35 GHz architecture and 1 terabyte of RAM) required to track the FL with the CP approximation. The ratio between the computing times without (values not reported in Table 3) and with the CP approximation is very close to the CP gain, g CP . This demonstrates that in the CP approach, the impact of the calculation of the polynomial coefficients (by pseudo-inversion) and the calculation of the approximated field along the line path are fairly negligible. As a consequence, with the reduction in the integration step size (useful for increasing the tracking accuracy) or the tracking of multiple FLs in the same bricks, there is a negligible impact on the computing burden. This further increases the appeal of the CP approximation in terms of practical applications.
The CP approach is even more efficient for FL tracking when a parallel computing system is available. As a matter of fact, in a tracking procedure, the field calculations need to be carried out in strictly sequential order because each point where the field is required is only revealed by the last integration step. Therefore, in a standard procedure, direct field calculations can only be made sequentially. On the other hand, in the CP approach, the direct field calculations are only required in the sampling points of the brick, all of them being already known. Therefore, a procedure based on parallel computing is optimal because it can treat all the sampling points together.
The last column of Table 3 also reports the times required to track the same FL by using an HPC system based on a GPU with an architecture of 2 × NVIDIA Tesla V100S 32 GB HBM2. The speeding-up due to the GPU, depending on the specific case, can achieve two orders of magnitude. It is worth noticing that the flat-time demand required in a number of different cases is due to the capability of taking advantage of the full performance of the parallel system.

Applicative Examples of Magnetic Field Lines Tracking
This section is devoted to assessing, by means of several examples, the accuracy of the FL tracking in both 2D (Two-Dimensional) and 3D (Three-Dimensional) when a CP approximation is applied. The same RK4 approach, with an integration step of h RK = 1 mm, is adopted. Standard bricks with sides of 10 × 5 × 5 cm in length, with 8 × 5 × 5 sampling points, are used for the CP approximation.
The first example analyses with the fully 3D CP approximation used a 2D axisymmetric plasma single null equilibrium configuration [21]. The set of sources includes both the CS and PF coils system, a suitable current able to generate an axisymmetric toroidal field of 6 T at [R, Z] = [2.11,0], and the plasma, modeled as a set of filamentary currents. For the sake of simplicity, the full set of plasma filaments is not reported, but only some synthetic plasma data are here reported: poloidal β p = 0.65, internal inductance l i = 0.8, plasma current I p = 5.5 MA. In Table 4, the total currents feeding the CS/PF coils system are also reported.
When a fully 2D source system is present, it is possible to define the poloidal flux ψ(R, Z) [Wb] as the magnetic flux linked with the axisymmetric circumference, crossing any poloidal section at a point with the same [R, Z] coordinates [8]. The plasma boundary is conventionally defined as the last closed magnetic surface (LCMS), characterized by a well-defined poloidal flux ψ b .
In Figure 5, the poloidal cross-section of the plasma boundary is reported (ψ b = 2.35 Wb). The LCMS cannot be crossed by any FL because it is a magnetic surface. As a consequence, an FL starting from a point of the plasma core is fully included within it. internal inductance = 0.8, plasma current = 5.5 MA. In Table 4, the total currents feeding the CS/PF coils system are also reported. When a fully 2D source system is present, it is possible to define the poloidal flux ( , ) [Wb] as the magnetic flux linked with the axisymmetric circumference, crossing any poloidal section at a point with the same , coordinates [8]. The plasma boundary is conventionally defined as the last closed magnetic surface (LCMS), characterized by a well-defined poloidal flux .
In Figure 5, the poloidal cross-section of the plasma boundary is reported ( = 2.35 Wb). The LCMS cannot be crossed by any FL because it is a magnetic surface. As a consequence, an FL starting from a point of the plasma core is fully included within it.  FL tracking with CP approximation has been applied to calculate the 3D capability. Therefore, due to integration errors, the calculated FL can leave the plasma surface. Then, its coordinates [R, ϕ, Z] (with ϕ being the cumulative toroidal angle) can refer to points with a poloidal flux value that is different from ψ b . An effective metric by which to estimate the integration accuracy could be the normalized flux discrepancy: (9) Figure 6 shows that, in the examined case, ∆ψ is less than 0.1%, corresponding to a displacement smaller than 1 mm, within a ϕ FL path of ϕ = 2000 degrees, corresponding to five-and-a-half toroidal turns.
FL tracking with CP approximation has been applied to calculate the 3D capability. Therefore, due to integration errors, the calculated FL can leave the plasma surface. Then, its coordinates , , (with being the cumulative toroidal angle) can refer to points with a poloidal flux value that is different from . An effective metric by which to estimate the integration accuracy could be the normalized flux discrepancy: (9) Figure 6 shows that, in the examined case, Δ is less than 0.1%, corresponding to a displacement smaller than 1 mm, within a FL path of = 2000 degrees, corresponding to five-and-a-half toroidal turns. The procedure is able to fully analyze 3D magnetic maps. To assess such capabilities, the effect on the FL trajectory (described above), due to the rigid displacement of one PF coil in the direction, has been analyzed. This is a typical analysis that is required to evaluate the impact of construction and assembly uncertainties on a field map [22][23][24].
Due to the 3D characteristics of the perturbation, the FL leaves the 2D plasma boundary. The best way to appreciate the entity of the FL deformation is the Poincaré map, which collects the crossing points of the FL with an assigned poloidal section. In Figure 7, the footprints on the = 0 poloidal section of the FL passing through the point , are reported for the different displacements of the PF6 coil. As expected, the greater the coil displacement, the larger the distance of the FL from the plasma boundary. The procedure is able to fully analyze 3D magnetic maps. To assess such capabilities, the effect on the FL trajectory (described above), due to the rigid displacement of one PF coil in the x direction, has been analyzed. This is a typical analysis that is required to evaluate the impact of construction and assembly uncertainties on a field map [22][23][24].
Due to the 3D characteristics of the perturbation, the FL leaves the 2D plasma boundary. The best way to appreciate the entity of the FL deformation is the Poincaré map, which collects the crossing points of the FL with an assigned poloidal section. In Figure 7, the footprints on the ϕ = 0 poloidal section of the FL passing through the point P 0 , are reported for the different displacements of the PF6 coil. As expected, the greater the coil displacement, the larger the distance of the FL from the plasma boundary. Due to its efficiency, this procedure can be used effectively to follow the FL path for long distances. To assess such a capability, a number of FLs have been tracked at t = 0 s in the pre-magnetization phase of a plasma scenario, when the plasma current is zero.
Only the upper side of the CS and PF coil currents are reported in Table 5 because of Due to its efficiency, this procedure can be used effectively to follow the FL path for long distances. To assess such a capability, a number of FLs have been tracked at t = 0 s in the pre-magnetization phase of a plasma scenario, when the plasma current is zero.
Only the upper side of the CS and PF coil currents are reported in Table 5 because of the up-down symmetry, and a toroidal axisymmetric field able to provide a field of 6 T at R = 2.11 m has been used. In Figure 8a, for each FL, the corresponding Poincaré diagram is reported, together with the total number of toroidal turns performed before leaving the vessel. Moreover, in Figure 8b, the 3D map is also shown. The computing burden, when using the abovedescribed parallel architecture, is in the order of a few tens of seconds for each toroidal turn.

Treatment of Periodicity
The CP approach can take advantage of a possible geometrical periodicity. In fact, the evaluation of the field in the sampling points of the interpolation bricks, which is the most demanding part of the process in terms of computing time, can be limited to only one period and the results can be exploited everywhere.
To assess such an opportunity, the TF coils system of the DTT has been analyzed by only considering a toroidal 20° slice. A set of 100 FLs has been tracked for a full toroidal turn, in order to test the actual periodicity of the integrated FL. The RK4 integration step considered is ℎ = 1 mm, as in the previous examples.
It is not straightforward to provide a picture of the lines. Therefore, only the maximum displacement of the FLs has been considered. The contour map reported in

Treatment of Periodicity
The CP approach can take advantage of a possible geometrical periodicity. In fact, the evaluation of the field in the sampling points of the interpolation bricks, which is the most demanding part of the process in terms of computing time, can be limited to only one period and the results can be exploited everywhere.
To assess such an opportunity, the TF coils system of the DTT has been analyzed by only considering a toroidal 20 • slice. A set of 100 FLs has been tracked for a full toroidal turn, in order to test the actual periodicity of the integrated FL. The RK4 integration step considered is h RK = 1 mm, as in the previous examples.
It is not straightforward to provide a picture of the lines. Therefore, only the maximum displacement of the FLs has been considered. The contour map reported in Figure 9 shows how the maximum deformation affects only those lines closer to the outboard of the TF coils.

Treatment of Periodicity
The CP approach can take advantage of a possible geometrical periodicity. In fact, the evaluation of the field in the sampling points of the interpolation bricks, which is the most demanding part of the process in terms of computing time, can be limited to only one period and the results can be exploited everywhere.
To assess such an opportunity, the TF coils system of the DTT has been analyzed by only considering a toroidal 20° slice. A set of 100 FLs has been tracked for a full toroidal turn, in order to test the actual periodicity of the integrated FL. The RK4 integration step considered is ℎ = 1 mm, as in the previous examples.
It is not straightforward to provide a picture of the lines. Therefore, only the maximum displacement of the FLs has been considered. The contour map reported in Figure 9 shows how the maximum deformation affects only those lines closer to the outboard of the TF coils.

Discussion
The accuracy provided by field approximation techniques depends on the regularity of the field and the characteristics of the approximating functions. The cost of calculating the approximate values depends on the size of the domain, the required accuracy, and the number of field points necessary for numerical integration.
Chebyshev polynomials are well-suited for describing magnetic fields as they are able to provide good accuracy when using limited-degree polynomials because of their capabilities to fit the magnetic field features.

Conclusions
FL tracking is a typical application requiring a high level of field calculations. CP interpolation is very effective in this class of applications, allowing for a strong reduction of the computing burden while guaranteeing high accuracy standards. The effectiveness of CP interpolation can be increased by using parallel architecture since the required direct field calculations can be made together because the sampling points in the Chebyshev domain are known on an a priori basis.
Funding: This paper is partially supported by the Italian MUR, PRIN grant #20177BZMAH, in progress among Italian Universities in cooperation with international labs.