PRO GAUSSIAN_CMB, cmb

  COMMON psm, param
  COMMON sky, model

  cmbdir = GETPSMDIR('cmb')

  ;; Filenames for the theory alm and map 
  ;;--------------------------------------
  cmb_unlensedalm_file = cmbdir + param.names.cmb_unlensed_alm_file
  cmb_unlensedmap_file = cmbdir + param.names.cmb_unlensed_map_file
  cmb_lensedalm_file = cmbdir + param.names.cmb_lensed_alm_file
  cmb_lensedmap_file = cmbdir + param.names.cmb_lensed_map_file
  potential_alm_file = cmbdir + param.names.potential_alm_file 
  potential_map_file = cmbdir + param.names.potential_map_file

  ;; Read unlensed and lensed Cl
  ;;-----------------------------
  cmb_lensedcl_file = cmbdir + 'cmb' + '_lensed_cl.fits' 
  cmb_unlensedcl_file = cmbdir + 'cmb' + '_unlensed_cl.fits' 
  PSM_DEBUG_MESSAGE, 'Fix problem with the format of the fits file', routine='GAUSSIAN_CMB'
  ;stop
  READ_CL, cmb_lensedcl_file, lensedcl, hdr=hdrlensed
  READ_CL, cmb_unlensedcl_file, unlensedcl, hdr=hdrunlensed

  ;; Create cmb component structure
  ;;-----------------------------------
  cmb_id = STRTRIM(SXPAR(hdrunlensed,'PSM_CPID'),2)
  cmb = COMPONENT_STRUCT('cmb',name='cmb',id=cmb_id)
  cmb.polarised = (model.global.sky_fields NE 'T')
  cmb.e1.ampl.unit = model.units.cmb
  *cmb.e1.numin.value = 0.
  *cmb.e1.numax.value = FLOAT('inf')

  polar = FIX(cmb.polarised)

  ;; Select the cl to be used to generate the CMB map and the name of
  ;; the output cmb alm 
  ;;------------------------------------------------------------------

  CASE model.cmb.cmb_lensing OF
     'cl': BEGIN
        cl = lensedcl
        alm_writefile = cmb_lensedalm_file
     END
     ELSE: BEGIN
        cl = unlensedcl
        alm_writefile = cmb_unlensedalm_file
     END
  ENDCASE

  ;;--------------------------------------------------------------------------------------------------------------------;;
  ;; Convolve cl with beam and pixel window function
  ;;--------------------------------------------------------------------------------------------------------------------;;
  ;; Get gaussian beam 
  IF NOT DEFINED(skyres) THEN skyres = SKY_RESOLUTION()
  ;;bl = GAUSSBEAM(skyres, MAX(cl.L), 2*polar+1)
  bl = GAUSSBEAM(skyres, MAX(cl.L), 3)

  ;; Get pixel window function 
  IF model.global.sky_pixwindow NE 0 THEN BEGIN
    tmpwl = HEALPIXWINDOW(skynside)
    tmpwl = tmpwl(*,0:2*polar)
    nwl=N_ELEMENTS(tmpwl[*,0])-1 < MAX(cl.L)
    wl = REPLICATE(0.d0, MAX(cl.L)+1, 2*polar+1)
    wl(0:nwl,*) = tmpwl(0:nwl,*)    
    DELVARX, tmpwl
    bl=bl[cl.L,*]*wl[cl.L,*]
  ENDIF

  ;; Apply beam to fields containing T
  wht = WHERE(cl.field EQ 'T', nwht, complement=whnt, ncomplement=nwhnt)
  IF nwht EQ 1 THEN BEGIN
     cmd = 'cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(wht[0]+1)+' = cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(wht[0]+1)+'*bl[cl.L,0]*bl[cl.L,0]' & PRINT, cmd
     ok = EXECUTE(cmd) & IF NOT ok THEN STOP
     IF nwhnt GT 0 THEN BEGIN
        FOR ii=0, nwhnt-1 DO BEGIN
           IF wht[0] LT whnt[ii] THEN BEGIN
              cmd = 'cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(whnt[ii]+1)+' = cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(whnt[ii]+1)+'*bl[cl.L,0]' & PRINT, cmd
              ok = EXECUTE(cmd) & IF NOT ok THEN STOP
           ENDIF ELSE BEGIN
              cmd = 'cl.c_'+NUMBER2STRING(whnt[ii]+1)+'_'+NUMBER2STRING(wht[0]+1)+' = cl.c_'+NUMBER2STRING(whnt[ii]+1)+'_'+NUMBER2STRING(wht[0]+1)+'*bl[cl.L,0]' & PRINT, cmd
              ok = EXECUTE(cmd) & IF NOT ok THEN STOP
           ENDELSE
        ENDFOR
     ENDIF
  ENDIF

  ;; Apply beam to fields containing E
  wht = WHERE(cl.field EQ 'E', nwht, complement=whnt, ncomplement=nwhnt)
  IF nwht EQ 1 THEN BEGIN
     cmd = 'cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(wht[0]+1)+' = cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(wht[0]+1)+'*bl[cl.L,1]*bl[cl.L,1]' & PRINT, cmd
     ok = EXECUTE(cmd) & IF NOT ok THEN STOP
     IF nwhnt GT 0 THEN BEGIN
        FOR ii=0, nwhnt-1 DO BEGIN
           IF wht[0] LT whnt[ii] THEN BEGIN
              cmd = 'cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(whnt[ii]+1)+' = cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(whnt[ii]+1)+'*bl[cl.L,1]' & PRINT, cmd
              ok = EXECUTE(cmd) & IF NOT ok THEN STOP
           ENDIF ELSE BEGIN
              cmd = 'cl.c_'+NUMBER2STRING(whnt[ii]+1)+'_'+NUMBER2STRING(wht[0]+1)+' = cl.c_'+NUMBER2STRING(whnt[ii]+1)+'_'+NUMBER2STRING(wht[0]+1)+'*bl[cl.L,0]' & PRINT, cmd
              ok = EXECUTE(cmd) & IF NOT ok THEN STOP
           ENDELSE
        ENDFOR
     ENDIF
  ENDIF

  ;; Apply beam to fields containing B
  wht = WHERE(cl.field EQ 'B', nwht, complement=whnt, ncomplement=nwhnt)
  IF nwht EQ 1 THEN BEGIN
     cmd = 'cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(wht[0]+1)+' = cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(wht[0]+1)+'*bl[cl.L,2]*bl[cl.L,2]' & PRINT, cmd
     ok = EXECUTE(cmd) & IF NOT ok THEN STOP
     IF nwhnt GT 0 THEN BEGIN
        FOR ii=0, nwhnt-1 DO BEGIN
           IF wht[0] LT whnt[ii] THEN BEGIN
              cmd = 'cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(whnt[ii]+1)+' = cl.c_'+NUMBER2STRING(wht[0]+1)+'_'+NUMBER2STRING(whnt[ii]+1)+'*bl[cl.L,2]' & PRINT, cmd
              ok = EXECUTE(cmd) & IF NOT ok THEN STOP
           ENDIF ELSE BEGIN
              cmd = 'cl.c_'+NUMBER2STRING(whnt[ii]+1)+'_'+NUMBER2STRING(wht[0]+1)+' = cl.c_'+NUMBER2STRING(whnt[ii]+1)+'_'+NUMBER2STRING(wht[0]+1)+'*bl[cl.L,0]' & PRINT, cmd
              ok = EXECUTE(cmd) & IF NOT ok THEN STOP
           ENDELSE
        ENDFOR
     ENDIF
  ENDIF  

  ;;--------------------------------------------------------------------------------------------------------------------;;
  ;; Get constraining alm and Nl if necessary
  ;;--------------------------------------------------------------------------------------------------------------------;;

  IF model.cmb.cmb_constrained EQ 'yes' THEN BEGIN

     ;; Get NILC5 alm and noise power spectrum, at the resolution of
     ;; the PSM sky, with the proposer window function
     ;;--------------------------------------------------------------
     lmaxnilc = (SKY_LMAX() < 1000)
     almobs = NILC5_ALM(beam_out=SKY_RESOLUTION(),lmax=lmax,pixwin_out=model.global.sky_pixwindow)
     nlnilc = NILC5_NL(LINDGEN(lmaxnilc+1L),beam_out=SKY_RESOLUTION(),pixwin_out=model.global.sky_pixwindow)
     nlobs = CL_STRUCTURE(nfields=1,fieldnames='T',units='uK_CMB',lmin=0, lmax=lmaxnilc, double=double)
     nlobs.c_1_1[*] = nlnilc
  ENDIF
  

  ;;--------------------------------------------------------------------------------------------------------------------;;
  ;; Generate CMB alm 
  ;;--------------------------------------------------------------------------------------------------------------------;;

  ;; Set sky resolution
  skyres = SKY_RESOLUTION()
  ;; Set sky nside
  skynside = SKY_NSIDE()
  ;; Set sky lmax
  skylmax = SKY_LMAX()
  ;; Set sky pixelisation
  skypix = SKY_PIX()

  CASE skypix.pixtype OF

    'HEALPIX': BEGIN 

       IF model.cmb.cmb_constrained EQ 'no' THEN BEGIN
          PSM_CL2ALM, cl, alm, lmin=2, lmax=MAX(cl.L) < skylmax, double=double
       ENDIF ELSE BEGIN
          PSM_CL2ALM, cl, alm, lmin=2, lmax=MAX(cl.L) < skylmax, almobs=almobs, nlobs=nlobs, double=double
       ENDELSE
       WRITE_CMB_ALM, alm(*,*,0:2*polar), alm_writefile, cmb
       *cmb.e1.ampl.value = alm_writefile

       ;;---------------------------------------;;
       ;; Write potential files if necessary
       ;;---------------------------------------;;
       IF model.cmb.cmb_lensing EQ 'ilens' THEN BEGIN
          ;; Write unlensed cmb map
          cmb_unlensedalm_file = cmbdir + param.names.cmb_unlensed_alm_file
          ;; cmb_unlensedmap_file = cmbdir + param.names.cmb_unlensed_map_file          
          ;; PSM_ALM2MAP, cmb_unlensedalm_file, cmb_unlensedmap_file, hpx_nside=skynside, polar=polar, double=double 
          ;; Write lensing potential alm and map files
          WRITE_ALM, alm(*,*,2*polar+1), potential_alm_file
          DELVARX, alm
          PSM_ALM2MAP, potential_alm_file, potential_map_file, hpx_nside=skynside, double=double
       ENDIF
    END ;; End of case of healpix maps

  ENDCASE
  ;;--------------------------------------------------------------------------------------------------------------------;;
  ;; 
  ;;--------------------------------------------------------------------------------------------------------------------;;  
END