## 1 Introduction
To date, the most credit on this subject is given to researchers who focused on the scale of crustal seismology and who were mainly active before the 1990s. However, with the flourishment of the surface-wave methods for near surface characterization [1; 2; 3; 4; 5; 6; 7], this subject has regained new interest in recent years. For example, [PERSON] recently applied the reduced delta matrix method along with the implicit function theory to achieve the analytical computation of PDs of Rayleigh-wave ellipticity [35] as well as the necessary phase velocity [36]. In order to increase the computational efficiency, he formulated the PDs of Rayleigh-wave phase velocity in a vector-transferring manner instead of the matrix multiplication. Nevertheless, this computational manner makes it intractable to further calculate the PDs of group velocity analytically since more complex second-order partial differentiations must be dealt with. In addition, he only considered the Rayleigh waves in the study. Due to the complexity of this analytical method, it is not widely used in current surface-wave linearized inversion strategies and many researchers have resorted to the global optimization methods as previously mentioned, bypassing the calculation of PDs. However, these methods require a large amount of calculation, particularly when thousands of DCs need to be inverted, as in the context of oil seismic exploration. Based on the above-mentioned reasons, here we rederive a set of new analytical expressions in the manner of matrix multiplication and present the relevant algorithms to compute the PDs of surface-wave DCs. This work can be taken as another option or a complement to the previous studies. We derive the relevant formulas and algorithms in such another way that is more convenient and more straightforward for considering all of the cases (Love and Rayleigh waves; phase and group velocity). Furthermore, we develop an open-source MATLAB package called SWPD (Surface Waves Partial Derivatives) accompanying this paper, and it has the following merits: (i) it is able to deal with both Rayleigh and Love waves; (ii) it can calculate the PDs of surface-wave phase-velocity DCs analytically or accomplish the calculations for Rayleigh-wave group-velocity DCs by a hemi-analytical method with high precision; (iii) all the relevant expressions and algorithms are derived from a widely accepted frame (i.e., the reduced delta matrix method), which makes the package easy to use and to be extended. It is worth mentioning that the hemi-analytical method is a compromise: it calculates all the relevant first-order PDs analytically and approximates the only second-order mixed PD term with a central finite-difference scheme. By doing so, we still can obtain satisfying results comparable to the analytical approach [32], meanwhile avoiding the intricate analytical derivations for the second-order term. We first outline the reduced delta matrix theory and its application in modeling the surface-wave phase-velocity DCs. We then present the algorithms and related formulas to calculate analytically the PDs of the phase velocity (Section 3). Next, both the fully analytical and the hemi-analytical approaches for calculating the PDs of the surface-wave group velocity are elaborated (Section 4). After that, we also describe the accompanying MATLAB code with focus on some key scripts or functions (Section 5). Last, four models are utilized for verification of our method and of the package (Section 6).