The analysis of the stability and numerical simulation of Costas loop circuits for high-frequency signals is a challenging task. The problem lies in the fact that it is necessary to simultaneously observe very fast time scale of the input signals and slow time scale of phase difference between the input signals. To overcome this difficult situation it is possible, following the approach presented in the classical works of Gardner and Viterbi, to construct a mathematical model of Costas loop, in which only slow time change of signal's phases and frequencies is considered. Such a construction, in turn, requires the computation of phase detector characteristic, depending on the waveforms of the considered signals. While for the stability analysis of the loop near the locked state (local stability) it is usually sufficient to consider the linear approximation of phase detector characteristic near zero phase error, the global analysis (stability in the large) cannot be accomplished using simple linear models. The present paper is devoted to the rigorous construction of nonlinear dynamical model of classical Costas loop, which allows one to apply numerical simulation and analytical methods (various modifications of absolute stability criteria for systems with cylindrical phase space) for the effective analysis of stability in the large. Here a general approach to the analytical computation of phase detector characteristic of classical Costas loop for periodic non-sinusoidal signal waveforms is suggested. The classical ideas of the loop analysis in the signal's phase space are developed and rigorously justified. Effective analytical and numerical approaches for the nonlinear analysis of the mathematical model of classical Costas loop in the signal's phase space are discussed.