The deployment of high power radio frequency waves in the ion cyclotron range (ICRF) constitutes an important operational facility in many plasma devices, including ITER. Any charged particle describes a helical motion around a given magnetic field line, the so-called cyclotron motion. ICRF relies on the interaction between charged particles and an injected Radio Frequency (RF) wave, tuned to be at the same frequency as the helical cyclotron motion. It is applied not only for pure heating of the plasma, i.e. Ion Cyclotron Resonant Heating (ICRH), but also for the generation of non-inductive current through Ion Cyclotron Current Drive (ICCD). The numerical code package SCENIC has been developed for self-consistently simulating the effects of ICRH on the resonant ion species within the plasma, the resulting changes in the plasma equilibrium, and finally the back reaction onto the injected wave field. SCENIC is an iterated scheme, which advances the resonant ions' distribution function, the equilibrium and the wave field iteratively until a converged solution, representing a steady state, is reached. The constituents of SCENIC are the MagnetoHydroDynamic (MHD) equilibrium code VMEC,1 the full wave code LE-Man2 and the Hamiltonian guiding centre drift following code VENUS.3 All of these codes are capable of dealing with 3D geometries, and have recently been updated to handle pressure anisotropy, where the energy density parallel and perpendicular to the magnetic field differ. This is important since the RF field resonates mainly with the particle's motion perpendicular to the magnetic field, thus creating pressure anisotropy. After the introduction and description of the different codes and their interfaces, this work verifies the consistency of the numerical results with expected results for simple cases, and a benchmarking effort against the similar code package SELFO is shown. After this validation, SCENIC is applied to different heating scenarios, which are relevant to present (Joint European Torus, JET) and future (ITER) devices. Low power heating simulations with a 1% helium-3 minority in background deuterium plasmas demonstrate that a pressure anisotropy is induced. We show that the hot particle distribution function can be adequately approximated with a particular bi-Maxwellian for the equilibrium and wave field computations. For high power, 3% hydrogen minority heating scenarios, the heating scheme alters the background equilibrium state. This justifies one of the main novelties introduced in this work, namely the inclusion of the equilibrium computation in the self-consistent scheme. Effects due to asymmetric wave injection and different heating locations on the hot particle distribution function, the hot dielectric tensor and the equilibrium will be studied. Here, the emergence of a high energy tail in the minority species distribution function is shown explicitly, and some of its exotic features are observed via the RF driven current, the density and the pressure evolution.