One of the limitations in phase field modeling historically has been a lack of ability to reproduce complex anisotropy in interface energies. Spherical Gaussians are able to produce continuous interface energies by subtracting a divot from a base value to reproduce experimentally or computationally calculated minima. The method is described for more simple 3D anisotropies, as would be seen in dendritic solidification, then is expanded to accurately reproduce 5D grain boundary energies to model polycrystalline grain growth (GG). In the 3D case, a base value is selected for a given system based on atomistic simulations or experimental data. Minima in energy are found with their corresponding directions, then a spherical Gaussian is subtracted out at each minimum direction to produce the desired total interface energy. In 5D quaternions, assigned to individual grains as orientations and as misorientations for grain boundaries, conduct the ongoing mesoscale changes. A 5D space scanning generates meaningful grain boundaries; input into the continuous function developed by Bulatov et al. to calculate grain boundary energy (GBE); which local minima are used in the phase field model. The methodology involves using 2D Gaussian switches, which match the misorientation between grains with misorientations for the GBE minima. Accounting a threshold range for the minimas, the switch activates a Spherical Gaussian to set the GBE to the desired value; creating in combination a full 5D GBE space. Multiphysics Object Oriented Simulation Environment (MOOSE), where reduced order parameters still retain individual grain identification useful for individually assigned quaternions, is used for implementation; with validation performed through bicrystal simulations of known outcomes.