PRO PSM_ALM2MAP, alm, map, output_beam=output_beam, ordering=ordering, hpx_nside=hpx_nside, $
                 hpx_pixelwindow=hpx_pixelwindow, polar=polar, double=double, overwrite=overwrite, nohdr=nohdr, help=help

  IF KEYWORD_SET(help) NE 0 THEN BEGIN
     PRINT, ''
     PRINT, PSM_STRING_TO_LENGTH('-',LINE_LENGTH(),character='-')
     PRINT, 'PSM_ALM2MAP: procedure to generate a map from alm coefficients'
     PSM_INFO_MESSAGE, 'SYNTAX: PSM_ALM2MAP, alm, map, output_beam=, ordering=, hpx_nside='
     PSM_INFO_MESSAGE, '                     /hpx_pixelwindow, /polar, /double, /nohdr, /overwrite'
     PSM_INFO_MESSAGE, "output_beam is the structure describing the output beam, if any."
     PSM_INFO_MESSAGE, "ordering is assumed to be 'RING' by default if not set"
     PSM_INFO_MESSAGE, "hpx_nside is the HEALPix nside parameter for the output map (required)"
     PSM_INFO_MESSAGE, 'set /hpx_pixelwindow to multiply the alm by the healpix pixel window function'
     PSM_INFO_MESSAGE, 'set /polar for polarised inverse transform'
     PSM_INFO_MESSAGE, 'set /double to force double precision for output alm'
     PSM_INFO_MESSAGE, 'set /overwrite to overwrite any already existing output'
     PSM_INFO_MESSAGE, 'set /nohdr to skip the writing of the correct header in the alm file'
     PRINT, PSM_STRING_TO_LENGTH('-',LINE_LENGTH(),character='-')
     RETURN
  ENDIF

  IF NOT KEYWORD_SET(ordering) THEN ordering='RING'
  
  IF PSM_DEBUG_VALUE() LE 1 THEN ON_ERROR, 2

  IF KEYWORD_SET(help) NE 0 THEN BEGIN
     PRINT, PSM_STRING_TO_LENGTH('-',LINE_LENGTH(),character='-')
     PRINT, 'PSM_ALM2MAP: procedure to calculate a map from alm coefficients'
     PSM_INFO_MESSAGE, 'SYNTAX: PSM_ALM2MAP, alm, map, /double'
     PSM_INFO_MESSAGE, 'set /double to force double precision for output map'
     PRINT, PSM_STRING_TO_LENGTH('-',LINE_LENGTH(),character='-')
     RETURN
  ENDIF

  IF NOT KEYWORD_SET(hpx_nside) THEN BEGIN
     MESSAGE, 'Need to define hpx_nside keyword', /info
     GOTO, closing
  ENDIF

  almtype = SIZE(alm, /type)


  ;; Test format of input alms, write alm file if needed
  ;;----------------------------------------------------
  IF almtype EQ 7 THEN BEGIN
     almfile = alm
     erasealmfile = 0
  ENDIF

  IF almtype LE 5 THEN BEGIN

     almfile = PSM_TMPFILENAME(extension='.fits')
     WRITE_ALM, alm, almfile
     erasealmfile = 1
  ENDIF

  ;; Test that the file exists
  exists = FILE_TEST(almfile)
  IF NOT exists THEN MESSAGE, 'Input alm file not found'
 
  ;; Test the type of the output map
  defmap = DEFINED(map)
  IF NOT defmap THEN BEGIN
     outname = PSM_TMPFILENAME(extension='.fits')
     readanderasemap = 1
  ENDIF ELSE BEGIN
     maptype = SIZE(map, /type)
     IF maptype NE 7 THEN BEGIN
        IF maptype LE 5 THEN BEGIN
           outname = PSM_TMPFILENAME(extension='.fits')
           readanderasemap = 1
           IF DEFINED(map) THEN BEGIN
              IF NOT KEYWORD_SET(overwrite) THEN BEGIN
                 MESSAGE, 'Variable for map already exists - set /overwrite to replace it by new map'
                 ;RETURN
              ENDIF ELSE BEGIN
                 DELVARX, map   ;saves memory
              ENDELSE
           ENDIF
        ENDIF ELSE MESSAGE, 'Case not implemented'
     ENDIF ELSE BEGIN
        outname = map
        readanderasemap = 0
        ;; Erase old file if it exists and if requested
        outputmapfileexists = FILE_TEST(map)
        IF outputmapfileexists THEN BEGIN
           IF NOT KEYWORD_SET(overwrite) THEN MESSAGE, map + ': File exists - set /overwrite keyword if necessary' $
           ELSE FILE_DELETE, outname, /quiet
        ENDIF
     ENDELSE
  ENDELSE

  ;; Write configuration file for alm2map_cxx
  ;;------------------------------------------
  conf = PSM_TMPFILENAME(extension='.conf')
  GET_LUN, unit
  OPENW, unit, conf
  IF NOT KEYWORD_SET(hpx_pixelwindow) THEN BEGIN
     PRINTF, unit,'pixel_window = false'       
  ENDIF ELSE BEGIN
     healpix_data=!HEALPIX_DIR+'/data/'
     IF PSM_VERBOSITY() GE 2 THEN PSM_INFO_MESSAGE, 'HEALPIX data directory: '+healpix_data
     PRINTF, unit,'pixel_window = true'
     PRINTF, unit,'healpix_data = '+ healpix_data
  ENDELSE
 
  ;; Find input and output beams, rebeam into a temporary file if needed
  filebeam = FITSFILE_BEAM(almfile)
  IF NOT KEYWORD_SET(output_beam) THEN output_beam = filebeam
  
  norebeam = STRUCT_EQUAL(filebeam, output_beam)
  IF norebeam THEN BEGIN
     almname = almfile
  ENDIF ELSE BEGIN
     temporaryalmname = PSM_TMPFILENAME(extension='.fits')
     almname = temporaryalmname
     CHANGE_ALM_RESOLUTION, almfile, almname, beam_in=filebeam, beam_out=output_beam, double=double
  ENDELSE
  
  ;; Output beam set to 0 (beam conversion done beforehand)
  PRINTF, unit, 'fwhm_arcmin = 0'  
  
  ;; lmax
  lmax = FITSFILE_LMAX(almfile)
  IF lmax EQ -1 THEN BEGIN
     ;; The lmax info is missing in the alm file. Get it from the size
     ;; of the input array
     PSM_DEBUG_MESSAGE, 'The lmax info is missing in the alm file', routine='PSM_ALM2MAP'
  ENDIF
  PRINTF, unit, 'nlmax = '+strtrim(lmax,2)
 
  PRINTF, unit, 'infile = '+ almname
  PRINTF, unit, 'nside = '+ strtrim(string(hpx_nside),2) 
  PRINTF, unit, 'outfile = ' + outname
  IF NOT DEFINED(polar) THEN polar =ISPOLARISED(almname,/alm)
  IF polar THEN BEGIN
     FITS_INFO, almname, n_ext=n_ext, /silent
     IF n_ext EQ 1 THEN BEGIN
        IF PSM_VERBOSITY() GE 1 THEN MESSAGE, 'Input file not polarised', /info
        IF PSM_CAREFULNESS() GE 2 THEN PSM_DEBUG_MESSAGE, 'Do not set /polar keyword when the input is not polarised', routine='PSM_ALM2MAP'
        PRINTF, unit, 'polarisation = false'
     ENDIF ELSE PRINTF, unit, 'polarisation = true'    
  ENDIF ELSE PRINTF, unit, 'polarisation = false'
  IF KEYWORD_SET(double) THEN $
     PRINTF, unit, 'double_precision = true'
  CLOSE, unit
  FREE_LUN, unit
  
  ;; Launch the C++ executable
  ;;---------------------------

  SPAWN, 'alm2map_cxx ' +conf
  PRINT, ' '

  ;; Check that it worked
  ;;----------------------
  ok = FILE_TEST(outname)
  IF NOT ok THEN STOP
  
  ;; Erase  the configuration file and the temporary alm file (if exists)
  ;;----------------------------------------------------------------------
  ;SPAWN, 'rm ' + conf
  FILE_DELETE, conf, /quiet
  IF NOT norebeam THEN FILE_DELETE, temporaryalmname, /quiet
  

  ;; If the output is not a file
  ;;-----------------------------
  IF readanderasemap EQ 1 THEN BEGIN

     hdmap = HEADFITS(outname,exten=1)
        
     SPLIT_FITSHDR, hdmap, hout=hout, hend=hend, hpsm=hpsm
     EXTRACT_PSMHDRBLOCK, hpsm, blockname='base', whereblock=whereblock
     IF whereblock[0] EQ -1 THEN hpsm = PSM_BASE_HDR(dataform='MAP', /setval)
     UPDATEHDRPAR, hpsm, 'PSM_DTFM', 'MAP'
     PSM_ID_KEYGEN, dtkey
     UPDATEHDRPAR, hpsm, 'PSM_DTID', dtkey
     PSM_ID_KEYGEN, flkey
     UPDATEHDRPAR, hpsm, 'PSM_FLID', flkey

     ;; Replace the alm header block by an appropriate map header block
     EXTRACT_PSMHDRBLOCK, hpsm, blockname='alm', wherenotblock=wherenotblock
     IF wherenotblock[0] NE -1 THEN hpsm = hpsm[wherenotblock]
     maphdrblock = PSM_MAP_HDR(pixtype='HEALPIX',lmax=lmax)
     UPDATE_PSMHDRBLOCK, hpsm, maphdrblock
     beamhdrblock = PSM_BEAM_HDR(output_beam)
     UPDATE_PSMHDRBLOCK, hpsm, beamhdrblock

     ;; Remove file ID from the header of the parent map, if any
     oldfile_id = SXPAR(hout,'PSM_FLID',count=c) 
     IF c EQ 1 THEN SXDELPAR, hout, 'PSM_FLID'
        
     ;; Create output header by concatenation
     hdmap = [hout,hpsm,hend]     
     REORDER_FITSHDR, hdmap

     MODFITS, outname, 0, hdmap, EXTEN_NO = 1

     READ_MAP, outname, map
     FILE_DELETE, outname, /quiet

     ;; If the output is expected to be nested, then reorder
     IF ordering eq 'NESTED' THEN map = REORDER(TEMPORARY(map), /r2n)

  ;; If the output is a file
  ;;-------------------------   
  ENDIF ELSE BEGIN
     
     fits_info, outname, n_ext=n, /silent
     IF n GT 1 THEN PSM_DEBUG_MESSAGE, 'More extensions than expected in output file', routine='PSM_ALM2MAP'

     hdalm = HEADFITS(almfile,exten=1)
     hdmap = HEADFITS(outname,exten=1)
     
     ;; If the output is expected to be nested, then reorder the output file
     ;;----------------------------------------------------------------------
     IF ordering eq 'NESTED' THEN BEGIN
        READ_MAP, outname, map
        map = REORDER(TEMPORARY(map), /r2n)
        UPDATEHDRPAR, hdmap, 'ORDERING', 'NESTED'
        FILE_DELETE, outname, /quiet
        WRITE_FITS_MAP, outname, map, ordering='NESTED'
     ENDIF

     ;; Fix the header
     ;;----------------
     SXDELPAR, hdmap, 'NLMAX'
     SXDELPAR, hdmap, 'NMMAX'
     SXDELPAR, hdmap, 'PIXWIN'
     SXDELPAR, hdmap, 'MAX_LPOL'
     SXDELPAR, hdmap, 'MAX_MPOL'
     wh = WHERE(STRTRIM(hdmap, 2) NE 'COMMENT End of imported HDU', nwh)
     hdmap = hdmap(wh)
     wh = WHERE(STRTRIM(hdmap, 2) NE 'COMMENT -----------------------', nwh)
     hdmap = hdmap(wh)
     wh = WHERE(STRTRIM(hdmap, 2) NE 'COMMENT  ** alm2map_cxx 1.0 **', nwh)
     hdmap = hdmap(wh)
     wh = WHERE(STRTRIM(hdmap, 2) NE 'COMMENT imported from HDU 2 of', nwh)
     hdmap = hdmap(wh)
     wh = WHERE(STRTRIM(hdmap, 2) NE 'COMMENT '+almname, nwh)
     hdmap = hdmap(wh)
     SPLIT_FITSHDR, hdmap, hout=hout, hend=hend, hpsm=hpsm    

     ;; If there is no PSM header in the map, then do not create it if /nohdr is set
     ;;------------------------------------------------------------------------------
     IF N_ELEMENTS(hpsm) LT 2 THEN BEGIN
        IF KEYWORD_SET(nohdr) THEN GOTO, finish
     ENDIF
     
     ;; If there is already a PSM header in hdmap, fix it, otherwise create it
     ;;------------------------------------------------------------------------
     IF N_ELEMENTS(hpsm) LT 2 THEN hpsm = [PSMHDR_HEADLINE(),PSMHDR_ENDLINE()]
     
     ;; Update the base header block
     EXTRACT_PSMHDRBLOCK, hpsm, blockname='base', whereblock=whereblock
     IF whereblock[0] NE -1 THEN BEGIN
        UPDATEHDRPAR, hpsm, 'PSM_DTFM', 'MAP'
        PSM_ID_KEYGEN, dtkey
        UPDATEHDRPAR, hpsm, 'PSM_DTID', dtkey
        PSM_ID_KEYGEN, flkey
        UPDATEHDRPAR, hpsm, 'PSM_FLID', flkey
     ENDIF ELSE BEGIN
        PSM_DEBUG_MESSAGE, 'Insert a base header block', routine='PSM_ALM2MAP'
        hpsm = PSM_BASE_HDR(dataform='MAP', /setval)
     ENDELSE
     
     ;; Replace the alm header block by an appropriate map header block
     EXTRACT_PSMHDRBLOCK, hpsm, blockname='alm', wherenotblock=wherenotblock
     IF wherenotblock[0] NE -1 THEN hpsm = hpsm[wherenotblock]
     maphdrblock = PSM_MAP_HDR(pixtype='HEALPIX',lmax=lmax)
     UPDATE_PSMHDRBLOCK, hpsm, maphdrblock
     beamhdrblock = PSM_BEAM_HDR(output_beam)
     UPDATE_PSMHDRBLOCK, hpsm, beamhdrblock
     
     ;; Create output header by concatenation
     hdmap = [hout,hpsm,hend]
     
     ;; Update units in standard fits header
     ;;--------------------------------------
     unit = SXPAR(hdalm,'TUNIT2', count=count)
     nfields = SXPAR(hdmap,'TFIELDS')
     IF count EQ 1 THEN FOR i=1,nfields DO SXADDPAR, hdmap, 'TUNIT'+NUMBER2STRING(i), unit
     
     REORDER_FITSHDR, hdmap
     MODFITS, outname, 0, hdmap, EXTEN_NO = 1
     
  ENDELSE
  
finish:
  IF erasealmfile EQ 1 THEN FILE_DELETE, almfile
  
closing:
END