pro mkhdr, header, im, naxisx, IMAGE = image, EXTEND = extend
 On_error,2

 npar = N_params()
 if npar LT 1 then begin
   print,'Syntax:  MKHDR, header, [ im, /IMAGE, /EXTEND ]'
   print,'    or   MKHDR, header, [ type, naxisx, /IMAGE, /EXTEND ]'
   print,'   header - output FITS header to be created'
   return
 endif

 if (npar eq 1) then begin               ;Prompt for keyword values
    read,'Enter number of dimensions (NAXIS): ',naxis
    s = lonarr(naxis+2)
    s[0] = naxis
    if ( naxis GT 0 ) then begin       ;Make sure not a dummy header
    for i = 1,naxis do begin       ;Get dimension of each axis
      keyword = 'NAXIS' + strtrim(i,2)
      read,'Enter size of dimension '+ strtrim(i,2) + ' ('+keyword+'): ',nx
      s[i] = nx
    endfor
  endif

  print,'Allowed datatypes are (1) Byte, (2) 16 bit integer, (3) 32 bit integer'
  print,'                      (4) 32bit floating, (5) 64 bit double precision'
  print,'                  or (14) 64bit integer'
  read,'Enter datatype: ',stype
  s[s[0] + 1] = stype

 endif else $
     if ( npar EQ 2 ) then s = size(im) $  ;Image array supplied
          else  s = [ N_elements(naxisx),naxisx, im ] ;Keyword values supplied

 stype = s[s[0]+1]              ;Type of data
        case stype of
        0:      message,'ERROR: Input data array is undefined'
        1:      bitpix = 8
        2:      bitpix = 16
        3:      bitpix = 32
        4:      bitpix = -32
        5:      bitpix = -64
        6:      message,'Complex types not allowed as FITS primary arrays'
        7:      bitpix = 8
       12:      bitpix = 16
       13:      bitpix = 32
       14:      bitpix = 64
        else:   message,'ERROR: Illegal Image Datatype'
        endcase

 header = strarr(s[0] + 7) + string(' ',format='(a80)')      ;Create empty array
 header[0] = 'END' + string(replicate(32b,77))

 if keyword_set( IMAGE) then $
    sxaddpar, header, 'XTENSION', 'IMAGE   ',' IMAGE extension' $
 else $
    sxaddpar, header, 'SIMPLE', 'T',' Written by IDL:  '+ systime()

 sxaddpar, header, 'BITPIX', bitpix, ' Number of bits per data pixel'
 sxaddpar, header, 'NAXIS', S[0],' Number of data axes'       ;# of dimensions

 if ( s[0] GT 0 ) then begin
   for i = 1, s[0] do sxaddpar,header,'NAXIS' + strtrim(i,2),s[i]
 endif

 if keyword_set( IMAGE) then begin
     sxaddpar, header, 'PCOUNT', 0, ' No Group Parameters'
     sxaddpar, header, 'GCOUNT', 1, ' One Data Group'
 endif else begin
     if keyword_set( EXTEND) or (s[0] EQ 0) then $
          sxaddpar, header, 'EXTEND', 'T', ' FITS data may contain extensions'
     Get_date, dte                       ;Get current date as CCYY-MM-DD
     sxaddpar, header, 'DATE', dte, $
           ' Creation UTC (CCCC-MM-DD) date of FITS header'
  endelse

 if stype EQ 12 then sxaddpar, header,'O_BZERO',32768, $
            ' Original Data is Unsigned Integer'
 if stype EQ 13 then sxaddpar, header,'O_BZERO',2147483648, $
            ' Original Data is Unsigned Long'
 header = header[0:s[0]+7]

 if not keyword_set(IMAGE) then begin   ;Add FITS definition for primary header
     sxaddpar,header,'COMMENT ', $
      "FITS (Flexible Image Transport System) format is defined in 'Astronomy"
     sxaddpar,header,'COMMENT ', $
      "and Astrophysics', volume 376, page 359; bibcode 2001A&A...376..359H"
 endif
 end