In this paper we present a numerical and step by step ray-tracing method on a phenomenological ionospheric electron density (N(e)) model, the TaiWan Ionospheric Model (TWIM), which is constructed from the FormoSat3/Constellation Observing System for Meteorology, Ionosphere and Climate (FS3/COSMIC) ionospheric radio occultation data. With the Earth's magnetic field and horizontal N(e) gradient effects included, efficient methods for calculating ray parameters such as the ground range, reflection height, phase path, and group path are presented. The three-dimensional TWIM consists of vertically fitted alpha-Chapman-type layers, with distinct F2, F1, E, and D layers, for which the layer parameters such as peak density, peak density height, and scale height are represented by surface spherical harmonics. This way the continuity of N(e) and its derivatives is maintained. This is important for practical schemes for providing reliable radio propagation predictions. The methodology is successfully applied to a practical high-frequency transmitter for oblique incidence ray tracing and further evaluated by comparing synthetic vertical ionograms generated by the method with experimental ionosonde observations.