code: def wcs_add_energy_axis(wcs, energies):
      
      if wcs.naxis != 2:
          raise Exception(
              'wcs_add_energy_axis, input WCS naxis != 2 %i' % wcs.naxis)
      w = WCS(naxis=3)
      w.wcs.crpix[0] = wcs.wcs.crpix[0]
      w.wcs.crpix[1] = wcs.wcs.crpix[1]
      w.wcs.ctype[0] = wcs.wcs.ctype[0]
      w.wcs.ctype[1] = wcs.wcs.ctype[1]
      w.wcs.crval[0] = wcs.wcs.crval[0]
      w.wcs.crval[1] = wcs.wcs.crval[1]
      w.wcs.cdelt[0] = wcs.wcs.cdelt[0]
      w.wcs.cdelt[1] = wcs.wcs.cdelt[1]
      w = WCS(w.to_header())
      w.wcs.crpix[2] = 1
      w.wcs.crval[2] = energies[0]
      w.wcs.cdelt[2] = energies[1] - energies[0]
      w.wcs.ctype[2] = 'Energy'
      return w
docs: Copy a WCS object, and add on the energy axis.
  
      Parameters
      ----------
      wcs : `~astropy.wcs.WCS`
          WCS
      energies : array-like
         Array of energies.