In the present work, the angular distribution of atmospheric cosmic muons is studied using PYTHIA8 simulations, comparing it to experimental data. The standalone code in PYTHIA8 consists of an atmosphere of only nitrogen molecules. Here, the standalone code has been updated adding oxygen molecules with a relative abundance of 78:22 (Nitrogen:Oxygen) to resemble a more realistic atmosphere. The vertical flux intensity ( \(I_\mu (\theta )\) ) of cosmic muons varies with zenith angle as a function of cos \(^n\theta \) , where \(\theta \)  is the zenith angle and n is the exponent term. Simulations are compared to the measurements performed at the surface laboratories of SINP (Kolkata) and JUSL (Jadugoda) at two different latitudes and altitudes. The integrated vertical muon flux ( \(I_{\mu }(0)\) ) and the exponent term (n) are extracted fitting each experimental dataset. The simulated vertical flux intensity ( \(I_\mu (\theta )\) ) as a function of the zenith angle ( \(\theta \) ) of cosmic muons shows a good level of agreement with the spectral shape of our experimental data. This study highlights the relevance of PYHTIA8 in the context of air-shower simulations where incorporation of a more realistic atmospheric model improves the model’s predictive power.