PRO SKY_IN_BANDS, bandlist=bandlist, banddirnames=banddirnames, bandnames=bandnames, bandunits=bandunits, polarised=polarised, double=double, $
                  group_gal=group_gal, group_strong_ps=group_strong_ps, group_faint_ps=group_faint_ps

  COMMON psm, param
  COMMON obs, obser

  IF NOT DEFINED(polarised) THEN polarised = SKY_POLARISED()
  IF NOT DEFINED(double) THEN double = PRECISION_ISDOUBLE()

  ;; Find existing components, select those to be observed
  ;;-------------------------------------------------------
  possible_comp = PSM_ALLCOMP()
  npc = N_ELEMENTS(possible_comp)
  existing_comp = EXISTING_COMPONENT_LIST(ncomp=nec)
  complist = existing_comp
  ncomp = N_ELEMENTS(complist)
  IF nec EQ 0 THEN BEGIN
     PSM_INFO_MESSAGE, 'No existing component in sky model'
     RETURN
  ENDIF

  tobe_observed = BYTARR(nec) & tobe_observed[*] = 1B
  PSM_DEBUG_MESSAGE, 'Fix list of observed and coadded components', routine='SKY_IN_BANDS', infolimit=3, stoplimit=3

  ;; Check whether band integration has already been done
  ;;------------------------------------------------------  

  nbands = N_ELEMENTS(bandlist)
  tobedone = BYTARR(nbands,ncomp)+1B
  
  IF param.tasks.sky_task EQ 'restore' THEN BEGIN

     tobedone = BYTARR(nbands,ncomp)

     tocheck = complist
     IF KEYWORD_SET(group_gal) THEN BEGIN
        FOR i=0, ncomp-1 DO BEGIN
           IF IS_GAL_COMP(complist[i]) THEN tocheck[i]='galaxy'
        ENDFOR
     ENDIF
     IF KEYWORD_SET(group_faint_ps) THEN BEGIN
        FOR i=0, ncomp-1 DO BEGIN
           IF IS_FAINTPS_COMP(complist[i]) THEN tocheck[i]='faintps'
        ENDFOR
     ENDIF
     IF KEYWORD_SET(group_strong_ps) THEN BEGIN
        FOR i=0, ncomp-1 DO BEGIN
           IF IS_STRONGPS_COMP(complist[i]) THEN tocheck[i]='strongps'
        ENDFOR
     ENDIF
     
     FOR i=0, N_ELEMENTS(bandlist)-1 DO BEGIN

        newband = *(bandlist[i])
        subdir = GETPSMDIR('skyinbands')+banddirnames[i]
        
        ;; Select polarisation options
        skypol = SKY_POLARISED()
        polar = polarised AND skypol
        
        ;; Units to use in the coaddition
        bandunit = bandunits[i]

        newinfo = {polar:polar,bandunit:bandunit,complist:complist,double:double,group_gal:group_gal,group_strong_ps:group_strong_ps,group_faint_ps:group_faint_ps}
        
        ;; Test whether the band directory exists, the band is the
        ;; same, and the component has already been 'bandintegrated'
        ok = FILE_TEST(FIX_SEPARATOR(subdir)+'band.sav')
        IF NOT ok THEN tobedone [i,*] = 1B ELSE BEGIN
           RESTORE, file=FIX_SEPARATOR(subdir)+'band.sav'
           IF NOT BANDS_EQUAL(newband,band) THEN tobedone[i,*] = 1B ELSE BEGIN
              ok = FILE_TEST(FIX_SEPARATOR(subdir)+'band.sav')
              IF NOT ok THEN tobedone[i,*] = 1B ELSE BEGIN
                 RESTORE, file=FIX_SEPARATOR(subdir)+'info.sav'
                 IF NOT STRUCT_EQUAL(newinfo, info) THEN BEGIN
                    tobedone[i,*] = 1B 
                    ;; empty the directory
                    FILE_DELETE, subdir, /quiet, /recursive
                 ENDIF ELSE BEGIN
                    ;; Loop on components
                    FOR j=0, ncomp-1 DO BEGIN
                       compname = complist[j]
                       neglect = 0B
                       PSM_DEBUG_MESSAGE, 'Revise neglecting of components', routine='SKY_IN_BANDS', infolimit=2
                       IF PSM_CAREFULNESS() LT 1 THEN neglect = COMPONENT_NEGLIGIBLE(compname,band)
                       IF NOT neglect THEN BEGIN
                          obscompmapname = FIX_SEPARATOR(subdir) + tocheck[j] + '_map_' + bandnames[i] + '.fits'
                          ok = FILE_TEST(obscompmapname)
                          IF IS_STRONGPS_COMP(compname) THEN BEGIN
                             IF NOT ok THEN IF (obser.what.strong_sources_to_map EQ 'yes') THEN tobedone[i,j] = 1B $
                             ELSE BEGIN
                                PSM_DEBUG_MESSAGE, 'Check that the component is in the file', routine='SKY_IN_BANDS', infolimit=2, stoplimit=2
                                ;; Check that the component is in the
                                ;; map (using CMIX in header)
                             ENDELSE
                             obscompcatname = FIX_SEPARATOR(subdir) + tocheck[j] + '_cat_' + bandnames[i] + '.sav'
                             ok = FILE_TEST(obscompcatname)
                             IF NOT ok THEN IF (obser.what.strong_sources_to_cat EQ 'yes') THEN tobedone[i,j] = 1B $
                             ELSE BEGIN
                                ;; Check that the component is in the
                                ;; catalogue (using CMIX in header)
                             ENDELSE
                          ENDIF ELSE BEGIN
                             IF NOT ok THEN tobedone[i,j] = 1B $ 
                             ELSE BEGIN
                                ;; Check that the component is in the
                                ;; map (using CMIX in header)
                             ENDELSE
                          ENDELSE
                       ENDIF
                    ENDFOR
                 ENDELSE
              ENDELSE
           ENDELSE
           DESTROY_DATA, band
        ENDELSE
     ENDFOR
           
  ENDIF

  IF MAX(tobedone) EQ 0 THEN GOTO, closing

  ;;-----------------------------------------------------------------------------------------------------------
  ;; Create beams on sky for strong sources
  ;;-----------------------------------------------------------------------------------------------------------
  
  ;; This needs to be done only if maps of strong sources are produced
  ;;-------------------------------------------------------------------
  
  CASE obser.what.strong_sources_to_map OF
     'yes': do_strong_sources = 1B
     'no': do_strong_sources = 0B
  ENDCASE

  IF do_strong_sources THEN BEGIN
        
     skyfwhm = SKY_RESOLUTION()
     
     ;; Infrared
     ;;----------
     strongirps_pointing_filename = PSM_TMPFILENAME(exten='.sav')
     whir = WHERE(complist EQ 'strongirps', nwhir)
     IF nwhir NE 0 THEN BEGIN
        RESTORE_COMP, comp, compname='strongirps'
        polar = polarised AND (comp.polarised)
        GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar
        SAVE, file=strongirps_pointing_filename, listpix, shape
        DESTROY_DATA, listpix
        DESTROY_DATA, shape
        DESTROY_DATA, comp
     ENDIF
     
     ;; Radio
     ;;-------
     strongradiops_pointing_filename = PSM_TMPFILENAME(exten='.sav')
     whrad = WHERE(complist EQ 'strongradiops', nwhrad)
     IF nwhrad NE 0 THEN BEGIN
        RESTORE_COMP, comp, compname='strongradiops'
        polar = polarised AND (comp.polarised)
        GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar
        SAVE, file=strongradiops_pointing_filename, listpix, shape
        DESTROY_DATA, listpix
        DESTROY_DATA, shape
        DESTROY_DATA, comp
     ENDIF
     
     ;; uchii
     ;;-------
     stronguchii_pointing_filename = PSM_TMPFILENAME(exten='.sav')
     whhii = WHERE(complist EQ 'stronguchii', nwhhii)
     IF nwhhii NE 0 THEN BEGIN
        RESTORE_COMP, comp, compname='stronguchii'
        polar = polarised AND (comp.polarised)
        GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar
        SAVE, file=stronguchii_pointing_filename, listpix, shape
        DESTROY_DATA, listpix
        DESTROY_DATA, shape
        DESTROY_DATA, comp
     ENDIF
     
     ;; wmap
     ;;-------
     strongwmapps_pointing_filename = PSM_TMPFILENAME(exten='.sav')
     whwmap = WHERE(complist EQ 'strongwmapps', nwhwmap)
     IF nwhwmap NE 0 THEN BEGIN
        RESTORE_COMP, comp, compname='strongwmapps'
        polar = polarised AND (comp.polarised)
        GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar
        SAVE, file=strongwmapps_pointing_filename, listpix, shape
        DESTROY_DATA, listpix
        DESTROY_DATA, shape
        DESTROY_DATA, comp
     ENDIF
     
     ;; ercsc
     ;;-------
     strongercscps_pointing_filename = PSM_TMPFILENAME(exten='.sav')
     whercsc = WHERE(complist EQ 'strongercscps', nwhercsc)
     IF nwhercsc NE 0 THEN BEGIN
        RESTORE_COMP, comp, compname='strongercscps'
        polar = polarised AND (comp.polarised)
        elliptic = (comp.elliptic)
        IF elliptic EQ 0 THEN BEGIN
           GAUSSIAN_BEAM_ON_SKY, skyfwhm, !dpi/2d0-(*comp.lat.value), (*comp.lon.value), listpix, shape, nside=SKY_NSIDE(), polar=polar
        ENDIF ELSE BEGIN
           ELLIPTIC_GAUSSIAN_BEAM_ON_SKY, ((*comp.el_fwhm_maj)+skyfwhm), (*comp.el_fwhm_min)+skyfwhm, $
                                          !dpi/2d0-(*comp.lat.value), (*comp.lon.value), (*comp.el_theta), listpix, shape, nside=SKY_NSIDE(), polar=polar
        ENDELSE
        SAVE, file=strongercscps_pointing_filename, listpix, shape
        DESTROY_DATA, listpix
        DESTROY_DATA, shape
        DESTROY_DATA, comp
     ENDIF
  ENDIF
  
  ;;-----------------------------------------------------------------------------------------------------------
  ;; Loop on bands
  ;;-----------------------------------------------------------------------------------------------------------

  FOR i=0, N_ELEMENTS(bandlist)-1 DO BEGIN

     IF MAX(tobedone[i,*]) EQ 0 THEN GOTO, skipthisband

     band = *(bandlist[i])
     
     ;; Create directory if it does not exist
     subdir = GETPSMDIR('skyinbands')+banddirnames[i]
     exists = FILE_TEST(subdir,/directory)
     IF NOT exists THEN FILE_MKDIR, subdir     

     ;; Save the band
     IF PSM_VERBOSITY() GE 1 THEN SAVE, file=FIX_SEPARATOR(subdir)+'band.sav', band, /verbose $
     ELSE SAVE, file=FIX_SEPARATOR(subdir)+'band.sav', band

     ;; Select polarisation options
     skypol = SKY_POLARISED()
     polar = polarised AND skypol

     ;; Units to use in the coaddition
     bandunit = bandunits[i]

     info = {polar:polar,bandunit:bandunit,complist:complist,double:double,group_gal:group_gal,group_strong_ps:group_strong_ps,group_faint_ps:group_faint_ps}
     IF PSM_VERBOSITY() GE 1 THEN SAVE, file=FIX_SEPARATOR(subdir)+'info.sav', info, /verbose $
     ELSE SAVE, file=FIX_SEPARATOR(subdir)+'info.sav', info

     ;;---------------------
     ;; Loop on components |
     ;;---------------------

     FOR j=0, ncomp-1 DO BEGIN

        IF tobedone[i,j] EQ 0B THEN GOTO, skipthiscomp

        compname = complist[j]
        neglect = 0B
        PSM_DEBUG_MESSAGE, 'Revise neglecting of components', routine='SKY_IN_BANDS'
        IF PSM_CAREFULNESS() LT 1 THEN neglect = COMPONENT_NEGLIGIBLE(compname,band)

        IF NOT neglect THEN BEGIN
           
           IF COMP_IS_CAT(compname) THEN BEGIN
              compol = ISPOLARISED(compname)
              ispol = (compol AND polarised)
              RESTORE_COMP, comp, compname=compname
              npol = 1+2*ispol
              channelbeam = SKY_RESOLUTION()          
              CASE STRMID(compname,0,5) OF
                 
                 ;; strong sources ------------------------------------------------------------------
                 ;; draw them directly on the map, or make a catalogue, or both

                 'stron': BEGIN
                    CASE obser.what.strong_sources_to_map OF
                       'yes': dopsmap = 1B
                       'no': dopsmap = 0B
                    END
                    CASE obser.what.strong_sources_to_cat OF
                       'yes': dopscat = 1B
                       'no': dopscat = 0B
                    END

                    ;; Catalogue of strong sources
                    ;;------------------------------

                    IF dopscat THEN BEGIN
                       catfilename = FIX_SEPARATOR(subdir) + compname + '_cat_' + bandnames[i] + '.sav'
                       IF ispol THEN datacat = MODEL_TO_DATA_CAT(comp, band, /polar) ELSE datacat = MODEL_TO_DATA_CAT(comp, band, /nopolar)
                       IF ispol THEN BEGIN
                          result = [[*datacat.obs1.iflux],[*datacat.obs1.qflux],[*datacat.obs1.uflux]]
                       ENDIF ELSE BEGIN
                          result = [*datacat.obs1.iflux]
                       ENDELSE
                       WRITE_CATALOG, datacat, catfilename, /sav, /name
                       DESTROY_DATA, datacat
                    ENDIF ELSE BEGIN
                       IF ispol THEN OBSERVE_MODEL_CAT, comp, band, result ELSE OBSERVE_MODEL_CAT, comp, band, result, /nopolar
                    ENDELSE

                    ;; Map of strong sources
                    ;;-----------------------
                    IF dopsmap THEN BEGIN
                       ok = EXECUTE('RESTORE, file='+compname+'_pointing_filename, /verbose')
                       nwhstrong = N_ELEMENTS(shape)
                       IF PRECISION_ISDOUBLE() THEN strongpsmap = DBLARR(SKY_NSIDE()^2*12L,npol) ELSE strongpsmap = FLTARR(SKY_NSIDE()^2*12L,npol)
                       FOR isource=0L, nwhstrong-1 DO BEGIN
                          lis = *(listpix[isource])
                          shpi = result[isource] * (*(shape[isource])).amplitude
if finite((*(shape[isource])).amplitude[0]) eq 0 then stop 
                          strongpsmap[lis,0] = strongpsmap[lis,0] + shpi  ;; check dim of result                    
                          IF ispol THEN BEGIN
                             ang = (*(shape[isource])).rotation + (*comp.pol_angle)[isource]
                             shpi = shpi*(*comp.pol_fraction)[isource]
                             shpq = shpi*COS(2*ang)
                             shpu = shpi*SIN(2*ang)
                             strongpsmap[lis,1] = strongpsmap[lis,1] + shpq                      
                             strongpsmap[lis,2] = strongpsmap[lis,2] + shpu                      
                          ENDIF
                       ENDFOR
                       DESTROY_DATA, listpix
                       DESTROY_DATA, shape
                       ;; convert map to channel units: map in mJy/sr, convert to band units
                       fac = CONVERSION_FACTOR('mJy/sr',bandunit, band, double=double)
                       strongpsmap = TEMPORARY(strongpsmap)*fac
                       compdata = PTR_NEW(strongpsmap,/no_copy)
                    ENDIF

                 END

                 ;; faint sources -------------------------------------------------------------------
                 ;; put them on the map and do harmonic space convolution
                 'faint': BEGIN   
                    IF ispol THEN OBSERVE_MODEL_CAT, comp, band, result, /approx ELSE OBSERVE_MODEL_CAT, comp, band, result, /approx, /nopolar
                    nsources = N_ELEMENTS(result[*,0])                    
                    IF PRECISION_ISDOUBLE() THEN faintpsmap = DBLARR(SKY_NSIDE()^2*12L,npol) ELSE faintpsmap = FLTARR(SKY_NSIDE()^2*12L,npol)
                    ANG2PIX_RING, SKY_NSIDE(), !dpi/2d0 - *comp.lat.value, *comp.lon.value, ipix
                    IF KEYWORD_SET(double) NE 0 THEN fourpi = 4d0*!dpi ELSE fourpi=4.*!pi 
                    omegapix = fourpi / NSIDE2NPIX(SKY_NSIDE())
                    FOR is=0L, nsources-1L DO faintpsmap[ipix[is],*] = faintpsmap[ipix[is],*] + result[is,*]
                    fac = CONVERSION_FACTOR('mJy/sr',bandunit,band, double=double)
                    fac = fac/omegapix
                    faintpsmap = TEMPORARY(faintpsmap)*fac
                    CHANGE_MAP_RESOLUTION, TEMPORARY(faintpsmap), newfaintpsmap, beam_in=GET_BEAM_STRUCT('GAUSSIAN',beamsize=0.,sizeunits='arcmin'), $
                                           beam_out=GET_BEAM_STRUCT('GAUSSIAN',beamsize=SKY_RESOLUTION(),sizeunits='arcmin'), /nohdr, lmax=SKY_LMAX()
                    compdata = PTR_NEW(newfaintpsmap,/no_copy)

                 END
              ENDCASE
                            
              DESTROY_DATA, comp

           ENDIF ELSE BEGIN
              ispol = ISPOLARISED(compname)
              COMP_IN_BAND, compname, band, compdata, units=bandunit, dtshp='array', dtfm='map', $
                            /pointer, polarised=ispol, hpx_nside=nside, lmax=lmax, double=double
           ENDELSE
              
           IF PTR_VALID(compdata) THEN BEGIN
              dimsout = SIZE(*compdata, /dimensions)
              ndimout = SIZE(*compdata, /n_dimensions)
              
              ;; Write individual component observation
              ;;----------------------------------------
              IF tobe_observed[j] EQ 1 THEN BEGIN
                 block1 = PSM_BASE_HDR(datatype='OBS',dataform='MAP',/setval)
                 block2 = PSM_BAND_HDR(band)
                 block3 = PSM_BEAM_HDR(SKY_RESOLUTION())
                 block4 = PSM_MAP_HDR(pixtype='HEALPIX',lmax=SKY_LMAX())                    
                 block5 = PSM_OBS_HDR(compname, obsid=obsid)
                                ;block6 = PSM_SKYMIX_HDR(compname)
                 MK_PSMHDR, xhdr, block1, block2, block3, block4, block5 ;, block6
                 ;; write the map in appropriate format
                 obscompmapname = FIX_SEPARATOR(subdir) + compname + '_map_' + bandnames[i] + '.fits'
                 skypix = SKY_PIX()
                 
                 IF ispol THEN WRITE_HPXMAP, obscompmapname, tqustokes=compdata, xhdr=xhdr, units=bandunit, coordsys=skypix.coordsys $
                 ELSE WRITE_HPXMAP, obscompmapname, tstokes=compdata, xhdr=xhdr, units=bandunit, coordsys=skypix.coordsys
              ENDIF                    
              
              ;; Free the pointer
              PTR_FREE, compdata
           ENDIF
     
        ENDIF
        skipthiscomp:
     ENDFOR                     ; (end loop on components)

     IF param.tasks.sky_task NE 'restore' THEN BEGIN
        IF KEYWORD_SET(group_gal) NE 0 THEN GROUP_GAL_INBAND, subdir;, /no_delete
        IF KEYWORD_SET(group_strong_ps) NE 0 THEN GROUP_STRONG_PS_INBAND, subdir;, /no_delete
        IF KEYWORD_SET(group_faint_ps) NE 0 THEN GROUP_FAINT_PS_INBAND, subdir;, /no_delete
     ENDIF ELSE BEGIN
        IF KEYWORD_SET(group_gal) NE 0 THEN BEGIN
           ok = FILE_TEST(FIX_SEPARATOR(subdir)+'galaxy_map_'+bandnames[i]+'.fits')
           IF NOT ok THEN GROUP_GAL_INBAND, subdir
        ENDIF
        IF KEYWORD_SET(group_strong_ps) NE 0 THEN BEGIN
           ok = FILE_TEST(FIX_SEPARATOR(subdir)+'strongps_map_'+bandnames[i]+'.fits')
           IF NOT ok THEN GROUP_STRONG_PS_INBAND, subdir   
        ENDIF                                   
        IF KEYWORD_SET(group_faint_ps) NE 0 THEN BEGIN
           ok = FILE_TEST(FIX_SEPARATOR(subdir)+'faintps_map_'+bandnames[i]+'.fits')
           IF NOT ok THEN GROUP_FAINT_PS_INBAND, subdir 
        ENDIF
     ENDELSE

     skipthisband:

     ;; Clean memory leaks
     ;;--------------------
     PSM_DEBUG_MESSAGE, 'Find memory leaks', routine='SKY_IN_BANDS'
     HEAP_GC

     ;; Make figures
     ;;--------------

  ENDFOR                        ; (end loop on bands)

closing:

;; Clean-up temporary pointing files for point sources
;;-----------------------------------------------------
 IF DEFINED(strongradiops_pointing_filename) THEN FILE_DELETE, strongradiops_pointing_filename, /quiet, /allow_nonexistent
 IF DEFINED(strongirps_pointing_filename) THEN FILE_DELETE, strongirps_pointing_filename, /quiet, /allow_nonexistent
 IF DEFINED(stronguchii_pointing_filename) THEN FILE_DELETE, stronguchii_pointing_filename, /quiet, /allow_nonexistent
 IF DEFINED(strongwmapps_pointing_filename) THEN FILE_DELETE, strongwmapps_pointing_filename, /quiet, /allow_nonexistent
 IF DEFINED(strongercscps_pointing_filename) THEN FILE_DELETE, strongercscps_pointing_filename, /quiet, /allow_nonexistent
END