In fitting the integral flux spectrum J( > E; B; L) of trapped energetic particles above the energy threshold E, it is common practice to use either an exponential function exp( - E/E0) or a power law E-γ where the energy constant E0 or the index γ are adjusted to obtain best fits to the measured flux values Ji( > Ei; B; L). These fit functions have generally been chosen not on firm physical grounds, but because of their mathematical simplicity. They have been used since 1960 in radiation belt modelling. However, a large number of mathematical functions can be selected to fit data with an equally good accuracy. Note that most of them have no physical meaning or relevance. In this work, the discrete AE-8 energy spectra Ji(Ei; B0; L) have been fitted with a sum of two maxwellian functions. The slope of the differential flux Jd(E) is then related to the characteristic temperature of the maxwellian velocity distribution and the normalisation constant is proportional to the number density of the maxwellian population. Besides its greater physical relevance, this new fit function has a non-singular behaviour in the limit E = 0. This is not the case for the power laws mentioned above which fit the data only in limited energy ranges. The fitting parameters have been determined for drift shells parameters ranging between L = 1.2 and L = 10. Their dependence on L is determined and discussed. The analytical fit of the energy spectra is good for energies smaller than 4 MeV, but it is not satisfactory at larger energies for reasons which are discussed. This work indicates that the energy spectra of the AE-8 model can be approximated rather well by two maxwellian distributions for electron energies smaller than 4 MeV, not only near geostationary orbit but also for a wide range of L values.