The Legendre polynomials in this form are orthogonal in the [-1, 1] window. Therefore, the scattering domain was first mapped into the interval [-1, 1] to calculate the coefficients, \(C_k\), and then the coefficients were converted to fit the correction in the scattering domain. 
Since for every set of parameter choice, there is a unique error function, a single general model can be obtained by fitting a multiple regression model to the coefficients of each Legendre polynomial term as a function of the physical parameters of the aerosols. As shown in Eq. (\ref{eq:P22-correctionMR}) coefficients can be a function of six physical parameters of the aerosol:  the number of monomers, monomer size parameter, refractive index of the monomers, and the average optical depths observed by the monomers within the aggregate. A general error function for each three scattering regimes can be written as;
Coefficients are fitted by multiple regressions which are optimized to give the maximum predicted and adjusted \(R^2\) with the simplest possible models while also maximizing the normality of the residuals of error. The optimization process involves both manually adding and removing parameters based on our expectations on each parameter and automatic statistical algorithms. More details about how the coefficients are fitted and validated are discussed in the Appendix. Table A.1 summarizes the multiple regression statistics.
The angular distribution of a correction function is highly sensitive to the accuracy of its coefficient terms, particularly \(C_1\)and \(C_2\). While this is a limitation of this method, it also provides a way to validate the reliability of the model by extrapolating to the outside of the tested range. Namely, overfitted or too simple models for the coefficients quickly result in unreasonable phase functions when extrapolated. In summary, each of the regression models for the coefficients is selected based on the standard regression criteria, expert judgment, predicted coefficient of validation method and consistency of the model for extrapolations. 
After the empirical corrections model error stays within a 10% limit for the run shown in Fig (\ref{434001}) while the error for the old parametrization model reaches up to 40%. The new model usually performs better than the case 
Coefficients are fitted by polynomials via a combination of stepwise regressions and manually adjusting the t  to add or remove polynomial terms based on procedure. In order to eliminate the overfitting or 
which are optimized to have the minimum possible residual errors and the maximum normality of the distribution of the residual errors.  The optimization procedure involves stepwise regression is made in a combination of manual elimination of some polynomial terms and an automated algorithm that uses Gaussian elimination process. 
For spherical particles, \(P_{22}=P_{11}\) and \(P_{44}=P_{33}\) while these terms are not equal for aggregates due to the depolarization caused by cross-polarization within aggregates  \citep{Tazaki_2016}.  Therefore, \(P_{11}\) is modeled by adding empirically calculated depolarization, \(depol\left(\theta\right)\), to \(P_{22}\) (Eq. (\ref{eq:P11})). This method produces \(P_{11}\) which fits well with the exact model (See Fig. (x)).  \(P_{44}\) is calculated in a similar fashion except that the depolarization is scaled to fit at \(\theta=0\degree\) and \(\theta=180\degree\) because of lack of agreement for the angular dependency of \(P_{44}\) with the exact calculations.