pro writefits, filename, data, header, heap, Append = Append,  $
       compress = compress, CheckSum = checksum, NaNValue = NaNvalue

  On_error, 2
  compile_opt idl2

  if N_params() LT 2 then begin
       print,'Syntax - WRITEFITS, filename, data,[ header, /APPEND, /CHECKSUM]'
       return
  endif

; Get information about data

  siz = size( data )
  naxis = siz[0]                    ;Number of dimensions
  if naxis GT 0 then nax = siz[ 1:naxis ]              ;Vector of dimensions
  lim = siz[ naxis+2 ]              ;Total number of data points
  type = siz[naxis + 1]             ;Data type

;Create a primary or image extension header if not supplied by the user

        if N_elements(header) LT 2 then begin
                if keyword_set(append) then mkhdr, header, data, /IMAGE  $
                                       else mkhdr, header, data, /EXTEND
        endif else if naxis GT 0 then $
              check_FITS, data, header, /UPDATE, /FITS

; Remove any STSDAS/random group keywords from the primary header

  hdr = header
  if not keyword_set( APPEND) then begin
         simple = 'SIMPLE  =                    T / Written by IDL:  ' $
                        + systime()
         hdr[0] =  simple + string( replicate(32b,80-strlen(simple) ) )
         sxdelpar, hdr, [ 'GCOUNT', 'GROUPS', 'PCOUNT', 'PSIZE' ]
  endif

; If necessary,convert unsigned to signed.    Do not destroy the original data

  if naxis NE 0 then begin

        unsigned = (type EQ 12) or (type EQ 13)
        if  unsigned then begin
             if type EQ 12 then begin
                     sxaddpar,hdr,'BZERO',32768,'Data is Unsigned Integer', $
                              before = 'DATE'
                     newdata = fix(data - 32768)
             endif else if type EQ 13 then begin
                    sxaddpar,hdr,'BZERO',2147483648,'Data is Unsigned Long', $
                              before = 'DATE'
                    newdata = long(data - 2147483648)
             endif
         endif

; For floating or double precision test for NaN values to write

  NaNtest = keyword_set(NaNvalue) and ( (type EQ 4) or (type EQ 5) )
  if NaNtest then begin
     NaNpts = where( data EQ NaNvalue, N_NaN)
     if (N_NaN GT 0) then begin
         if type EQ 4 then data[NaNpts]  = !Values.F_NaN    $
     else if type EQ 8 then data[NaNpts] = !Values.D_NaN
     endif
  endif
  endif


; Open file and write header information

        if keyword_set( APPEND) then begin
            if (strmid( hdr[0],0,8 ) NE 'XTENSION') then begin
                   message, $
            'ERROR - "XTENSION" must be first keyword in header extension',/CON
                  return
            endif
            test = file_search(filename, COUNT = n)
             if n EQ 0 then  begin       ;Create default primary header
                 mkhdr,h0,0,/exten
                 writefits,filename,0,h0, checksum = checksum
                 openu, unit, filename, /GET_LUN, /swap_if_little_endian
             endif else begin
            openu, unit, filename, /GET_LUN, /swap_if_little_endian
            mrd_hread, unit, hprimary
            extend = where( strmid(hprimary,0,8) EQ 'EXTEND  ', Nextend)
            if Nextend EQ 0 then $
               message,'WARNING - EXTEND keyword not found in primary FITS header',/CON
            endelse

            file = fstat(unit)
            nbytes  = file.size
            point_lun, unit, nbytes
            npad = nbytes mod 2880
            if npad NE 0 then writeu, unit, replicate(32b, 2880 - npad)

    endif else begin

        ext = ''
        if keyword_set(COMPRESS) then begin
            if strlowcase(strmid(filename,2,3,/reverse)) NE '.gz' $
               then ext = '.gz'
        endif else compress = 0


       openw, unit, filename + ext, /GET_LUN, /swap_if_little_endian, $
                             compress = compress

    endelse

; Determine if an END line occurs, and add one if necessary

       endline = where( strmid(hdr,0,8) EQ 'END     ', Nend)
     if Nend EQ 0 then begin

 message,'WARNING - An END statement has been appended to the FITS header',/INF
     hdr = [ hdr, 'END' + string( replicate(32b,77) ) ]
     endline = N_elements(hdr) - 1

   endif

; Add any CHECKSUM keywords if desired or already present

    do_Checksum = keyword_set(checksum)
    if not do_checksum then test = sxpar(hdr,'CHECKSUM',count=do_checksum)

     if do_checksum then begin
               if N_elements(heap) GT 0 then $
             FITS_ADD_CHECKSUM, hdr, [data,heap] else $
                 FITS_ADD_CHECKSUM, hdr, data
               endline = where( strmid(hdr,0,8) EQ 'END     ', Nend)
       endif
       nmax = endline[0] + 1

; Convert to byte and force into 80 character lines

       bhdr = replicate(32b, 80l*nmax)
       for n = 0l, endline[0] do bhdr[80*n] = byte( hdr[n] )
       npad = 80l*nmax mod 2880
       writeu, unit, bhdr
       if npad GT 0 then writeu, unit,  replicate(32b, 2880 - npad)

; Write data
       if naxis EQ 0 then goto, DONE
        bitpix = sxpar( hdr, 'BITPIX' )
        nbytes = long64( N_elements( data)) * (abs(bitpix) / 8 )
        npad = nbytes mod 2880

        if unsigned then writeu, unit, newdata $
                    else writeu, unit, data

; Write optional heap area
        if N_elements(heap) GT 0 then begin
              theap = sxpar(hdr,'THEAP', Count=N_Theap)
              if N_Theap GT 0 then begin
                  offset = theap - nbytes
                  if offset GT 0 then begin
                      writeu, unit, bytarr(offset)
                      npad = (npad + offset) mod 2880
                  endif
                  writeu, unit, heap
                  npad = (npad + N_elements(heap)) mod 2880
              endif
         endif

; ASCII Tables padded with blanks (32b) otherwise pad with zeros
        if keyword_set( APPEND) then begin
             exten = sxpar( header, 'XTENSION')
             padnum =  exten EQ 'TABLE   ' ? 32b : 0b
        endif else padnum = 0b

        if npad GT 0 then writeu, unit, replicate( padnum, 2880 - npad)
DONE:
        free_lun, unit

  return
  end