In this thesis I present particle in cell (PIC) simulations of ion cyclotron emission (ICE). ICE comprises suprathermal radiation in the ion cyclotron frequency range, whose spectrum peaks at successive local cyclotron harmonics of the emitting energetic ion population. ICE has previously been observed in all large toroidal MCF plasmas [Dendy and McClements, 2015; McClements et al., 2015]. ICE is caused by a collective instability, which in its linear phase corresponds to the magnetoacoustic cyclotron instability (MCI). The passive, non-invasive character of ICE measurements, suggest ICE is an attractive way forward for future energetic ion measurements in ITER. In the simulations in this thesis we use the EPOCH particle-in-cell code to solve the self-consistent Maxwell-Lorentz system of equations for fully kinetic electrons and thermal background ions, together with the minority energetic ion distribution that drives the primary ICE. We first perform a detailed quantitative comparison between fusion born proton driven chirping ICE observed during KSTAR ELM crashes and fully nonlinear direct numerical simulations of the MCI. We find good quantitative agreement between the simulated and observed spectra, to the extent that the simulations can be used to infer fast (-μs) time scale dynamics of the local electron number density in the emitting region. We then extend this study to determine the origin of a faint, time delayed proton chirping feature observed in one of the KSTAR plasma pulses. We do this using bicoherence analysis of both experimental and simulation data. We then run MCI PIC simulations of the pre ELM crash \steady state" ICE observed on KSTAR, which is believed to be driven by neutral beam injected (NBI) deuterons. PIC simulations of MCI excited ICE in JET and ASDEX Upgrade (AUG) plasmas are then discussed, and we show that AUG observations of the fundamental ICE harmonic can only be explained in terms of the MCI if nonlinear wave-wave interactions between higher harmonics are taken into account. Motivated by recent observations of ICE in the core region of several tokamaks, including AUG and DIIID, we then compare MCI simulations using two types of energetic ion distribution function, a spherical shell of varying thickness, and a ring beam of varying width. It is found that both distribution functions lead to MCI excited waves, and their nonlinear properties are discussed.