The observation of secondary γ-rays provides an alternative method of measuring cross sections that populate excited final states in nuclear reactions. The angular distributions of these γ-rays also provide information on the underlying reaction mechanism. Despite the large number of data of this type in the literature, publicly available R-matrix codes do not have the ability to calculate these types of angular distributions. In this paper, the mathematical formalism derived by C. R. Brune and R. J. deBoer [Phys. Rev. C 102, 024628 (2020)2469-998510.1103/PhysRevC.102.024628] is implemented in the R-matrix code azure2 and calculations are compared with previous data from the literature for the N15(p,α1γ)C∗12 reaction. In addition, new measurements, made at the University of Notre Dame Nuclear Science Laboratory using the Hybrid Array of Gamma Ray Detectors (HAGRiD), are reported that span the energy range from Ep=0.88 MeV to Ep=4.0MeV. Excellent agreement between the data and the phenomenological fit is obtained up to the limit of the previous fit at Ep=2.0MeV and the R-matrix fit is extended from Ex≈13.5 MeV up to Ex≈15.3 MeV, where N15+p and C12+α reactions are fit simultaneously for the first time. An excellent reproduction of the N15(p,α1γ)C∗12 and C12(α,α)C12 data is achieved, but inconsistencies and difficulty in fitting other data are encountered and discussed.