## 1 Introduction  
Investigating highly energetic material behavior is critical for various applications, ranging from designing safer explosive devices to understanding high-energy astrophysical events. Central to this investigation is the accurate modeling of the energetic material detonation reaction products, where the Jones-Wilkins-Lee (JWL) equation of state (EoS) is a widely adopted model used to describe their behavior [1, 2, 3, 4]. Upon detonation, an energetic material produces highly compressed and high-temperature gases and condensed products, with the density at the detonation wave front exceeding the original density of the explosive at the Chapman-Jouguet (CJ) state. This product mixture performs work on its surroundings under the hypothesis of an adiabatic process where no heat transfer occurs between the system and its environment [5]. The JWL EoS, hence, offers a semi-empirical description of the pressure-volume relationship of detonation products under isentropic conditions, valued for its relative simplicity and effectiveness. The detonation energy can be evaluated using the JWL model, which assumes that the detonation products' thermodynamics are entirely describable in a pressure-volume space. The CJ model of the detonation process assumes an instantaneous compression from ambient conditions to the CJ point, followed by expansion along the isentrope. Integrating the JWL equation yields the energy of the detonation products at any given relative volume. The coefficients in the JWL EoS are crucial for accurate detonation modeling. These coefficients are typically derived from cylinder test data by matching wall velocity-time histories using finite element hydrocodes [6, 7] and are subsequently refined through the nonlinear fitting of pressure-volume data [8]. Alternatively, theoretical evaluation involves thermochemical calculations to generate pressure-volume data along the expansion isentrope, fitting these data to the JWL equation to obtain the necessary coefficients. Despite its broad application, heuristic (such as particle swarm algorithms) or empirical methods [9] are usually applied to find the JWL parameters; the literature lacks a comprehensive and explicit description of a robust numerical least squares fitting algorithm specifically applied to find the coefficients of the JWL EoS. A detailed methodology for fitting the JWL EoS parameters using a nonlinear least squares (NLLS) approach [10] is presented. We provide a complete algorithmic framework for enforcing constraints at the CJ point, ensuring adherence to essential thermodynamic conditions. Furthermore, we delve into the distinct contributions of various terms in the JWL EoS, providing a granular analysis of their roles and impact on the overall equation. In this study, a Gauss-Newton method [11, 12] combined with a trust-region strategy [13] is employed to achieve rapid and robust convergence. The proposed algorithm constructs a sequence of feasible iterates, incorporating a trust-region subproblem that effectively balances the trade-off between fast convergence and stability. By applying singular value decomposition (SVD) [14], our method ensures numerical robustness, particularly in overdetermined systems, which is crucial for accurate parameter estimation. In such a way, even if the rank of the coefficient matrix is not full-rank, the solution of the nonlinear least squares problem can be found. The presented algorithm was implemented in HEMSim, a toolbox aimed at simulating highly energetic materials [15, 16]. This article is structured as follows. We begin with an in-depth explanation of our fitting algorithm; this is followed by a theoretical overview of the JWL EoS, highlighting the significance of precise parameter fitting: the imposition of the CJ constraint and the treatment of different terms within the JWL expression are explained.