The crystal-field splitting and Zeeman splitting of energy levels of Nd 3þ (4f 3) doped into semiconducting GaN (3.2 eV) grown in the hexagonal (huntite) phase by plasma-assisted molecular beam epitaxy have been modeled using a parameterized Hamiltonian defined to operate within the complete 4f 3 electronic configuration of Nd 3þ substituted for Ga 3þ in the lattice. Zeeman splittings were obtained by applying magnetic fields up to 6.6 T with the fields parallel and perpendicular to the crystallographic c-axis. The experimental energy (Stark) levels were obtained from a recent spectroscopic study on the same samples, where the combined excitation emission spectroscopy (CEES) identified the majority of Nd 3þ ions as replacing Ga 3þ in sites of C 3v symmetry. The manifolds of Nd 3þ (4f 3) 2Sþ1 L J modeled for the crystal-field splitting include the ground state, 4 I 9/2 , and excited states 4 I 11/2 , 4 I 13/2 , 4 F 3/2 , 4 F 5/2 , 2 H 9/2 , 4 F 7/2 , 4 S 3/2 , 4 G 5/2 , and 4 G 7/2. The energies of 41 experimental Stark levels from these manifolds were modeled through the use of a Monte Carlo method in which independent crystal-field parameters (CFP) were given random starting values and optimized using standard least-squares fitting between calculated and experimental Stark levels. Irreducible representations (irreps) and crystal field quantum numbers (l) were assigned to the energy level states of the 4 I 9/2 and 4 F 3/2 multiplet manifolds based on an analysis of the Zeeman data. This allows determination of which of the competing local minima should be considered to be the physically significant minimum. Using standard least-squares fitting between calculated and experimental Stark levels for Nd 3þ in C 3v symmetry, we obtain a final standard deviation of 7.01 cm À1 (rms ¼ 5.48 cm À1). V