In this work, the angular distribution of atmospheric cosmic muons is studied using PYTHIA8 simulations, comparing it to the experimental data. The standalone code in PYTHIA8 consists of an atmosphere of only nitrogen molecules. Here, the standalone code has been augmented adding oxygen molecules with a relative abundance of 78:22 (Nitrogen:Oxygen) to resemble a more realistic atmosphere. The vertical flux intensity ([Formula: see text]) of cosmic muons varies with zenith angle as a function of [Formula: see text], where [Formula: see text] is the zenith angle and n is the exponent term. The vertical flux intensity of cosmic muons has been characterized as a function of zenith angle ([Formula: see text]), momentum (p) and energy (E). Simulations are compared to the experimental results found in the literature and to the measurements performed at the surface laboratories of SINP (Kolkata) and JUSL (Jadugoda) at two different latitudes and altitudes in May 2024. The integrated vertical muon flux ([Formula: see text]) and the exponent term (n) are extracted by fitting each experimental dataset. The simulated vertical flux intensity ([Formula: see text]) as a function of the zenith angle ([Formula: see text]) of cosmic muons shows a good level of agreement with the spectral shape of our experimental data as well as with that from literature. This study highlights the relevance of PYHTIA8 in the context of air-shower simulations, improving its predictive power by incorporating a more realistic atmospheric model.