  We have established an iterative scheme to calculate with 15-digit accuracy
the numerical values of Ambartsumian-Chandrasekhar's H-functions for
anisotropic scattering characterized by the four-term phase function: the
method incorporates some advantageous features of the iterative procedure of
Kawabata (2015) and the double-exponential integration formula~(DE-formula) of
Takahashi and Mori (1974), which proved highly effective in Kawabata (2016).
Actual calculations of the H-functions have been carried out employing 27
selected cases of the phase function, 56 values of the single scattering albedo
$\varpi_0$, and 36 values of an angular variable $\mu(=\cos \theta)$, with
$\theta$ being the zenith angle specifying the direction of incidence and/or
emergence of radiation.
  Partial results obtained for conservative isotropic scattering, Rayleigh
scattering, and anisotropic scattering due to a full four-term phase function
are presented. As a sample application of the isotropic scattering H-function,
an attempt is made in Appendix to simulate by iteratively solving the
Ambartsumian equation the values of the plane and spherical albedos of a
semi-infinite, homogeneous atmosphere calculated by Rogovtsov and Borovik
(2016), who employed their analytical representations for these quantities and
the single-term and two-term Henyey-Greenstein phase functions of appreciably
high degrees of anisotropy, to find that our results are in satisfactory
agreement with theirs.
