Rational function-based models have proved to be very efficient for accurate frequency-dependent modeling of power system components. These models are able to characterize the components terminal behaviours (analysing the admittance matrix) for nodal analysis. This provides a fast convergence and inherent stability to the solution routine of the model. This work presents a general framework for interfacing the dynamic phasor method to the rational models. That would be promising for the electromagnetic transient analysis (under harmonic distortion), in the frequency domain. Therefore, Y-element rational pole-residue models (employing the vector fitting method) are developed. Moreover, the pole-residue model is converted into the state-space representation. Next, the dynamic harmonic approach (in the frequency domain) is employed for harmonic analysis. It is shown that the order of state-space system can become a major concern for frequency-dependent networks analysis. Therefore, to generate a reduced-order model, the balanced realization theory is applied. Moreover, (for the sake of simplicity and efficiency) the trapezoidal integration rule is employed to discretise the state-space equations. For validation of the modelling, it is applied on three test case studies and results of these studies are compared with their time-domain analysis results.