mwrfits.pro

来自「basic median filter simulation」· PRO 代码 · 共 1,696 行 · 第 1/4 页

    if keyword_set(lscale) then      array = long(array) $
    else if keyword_set(iscale) then array = fix(array)  $
    else                             array = byte(array)
    
end

; Write a header
pro mwr_header, lun, header

    ; Fill strings to at least 80 characters and then truncate.

    space = string(replicate(32b, 80))
    header = strmid(header+space, 0, 80)

    w = where(strmid(header,0,8) eq "END     ")

    if w[0] eq -1 then begin

       header = [header, strmid("END"+space,0,80)]
       
    endif else begin
        if (n_elements(w) gt 1) then begin 
           ; Get rid of extra end keywords;
           print,"MWRFITS Warning: multiple END keywords found."
           for irec=0L, n_elements(w)-2 do begin
              header[w[irec]] = strmid('COMMENT INVALID END REPLACED'+  $
                space, 0, 80)
           endfor
       endif

       ; Truncate header array at END keyword.
       header = header[0:w[n_elements(w)-1]]
    endelse

    nrec = n_elements(header)
    if nrec mod 36 ne 0 then header = [header, replicate(space,36 - nrec mod 36)]

    writeu, lun, byte(header)
end


; Move the group information within the data.
pro mwr_groupinfix, data, group, hdr

    siz = size(data)
    sizg = size(group)

    ; Check if group info is same type as data 

    if siz[siz[0]+1] ne sizg[3] then begin
        case siz[siz[0]+1] of
         1: begin
               mwr_groupscale, 127.d0, group, hdr
               group = byte(group)
           end
         2: begin
               mwr_groupscale, 32767.d0, group, hdr
               group = fix(group)
           end
         3: begin
               mwr_groupscale, 2147483647.d0, group, hdr
               group = long(group)
           end
         4: group = float(group)
         5: group = double(group)
      else: begin
                print,'MWRFITS Internal error: Conversion of group data'
               return
            end
        endcase
    endif

    nrow = 1
    for i=1, siz[0]-1 do begin
        nrow = nrow*siz[i]
    endfor

    data = reform(data, siz[siz[0]+2])
    for i=0L, siz[siz[0]] - 1 do begin
        if i eq 0 then begin
            gdata = group[*,0]
           gdata = reform(gdata)
            tdata = [ gdata , data[0:nrow-1]]
        endif else begin
            start = nrow*i
           fin = start+nrow-1
           gdata = group[*,i]
            tdata = [tdata, gdata ,data[start:fin]]
       endelse
    endfor

    data = temporary(tdata)
end

; If an array is being scaled to integer type, then
; check to see if the group parameters will exceed the maximum
; values allowed.  If so scale them and update the header.
pro mwr_groupscale, maxval, group, hdr

    sz = size(group)
    for i=0L, sz[1]-1 do begin
         pmax = max(abs(group[i,*]))
         if (pmax gt maxval) then begin
             ratio = pmax/maxval
            psc = 'PSCAL'+strcompress(string(i+1),/remo)
            currat = fxpar(hdr, psc)
            if (currat ne 0) then begin
                fxaddpar, hdr, psc, currat*ratio, 'Scaling overriden by MWRFITS'
            endif else begin
                fxaddpar, hdr, psc, ratio, ' Scaling added by MWRFITS'
            endelse
             group[i,*] = group[i,*]/ratio
         endif
    endfor
end
        
        
; Write out header and image for IMAGE extensions and primary arrays.
pro mwr_image, input, siz, lun, bof, hdr,       $
       null=null,                              $
       group=group,                            $
       pscale=pscale, pzero=pzero,             $
       lscale=lscale, iscale=iscale,              $
       bscale=bscale,                          $
        no_comment=no_comment,                  $
       silent=silent


    type = siz[siz[0] + 1]

    bitpixes=[8,8,16,32,-32,-64,-32,0,0,-64,0,0,16,32,64,64]

    ; Convert complexes to two element real array.

    if type eq 6 or type eq 9 then begin
 
        if not keyword_set(silent) then begin
            print, "MWRFITS Note: Complex numbers treated as arrays"
        endif
    
        array_dimen=(2)
        if siz[0] gt 0 then array_dimen=[array_dimen, siz[1:siz[0]]] 
        if siz[siz[0]+1] eq 6 then data = float(input,0,array_dimen)  $
        else data = double(input,0,array_dimen)

    ; Convert strings to bytes.
    endif else if type eq 7 then begin
        data = input
        len = max(strlen(input))
        if len eq 0 then begin
            print, 'MWRFITS Error: strings all have zero length'
           return
        endif

        for i=0L, n_elements(input)-1 do begin
            t = len - strlen(input[i])
           if t gt 0 then input[i] = input[i] + string(replicate(32B, len))
        endfor
    
        ; Note that byte operation works on strings in a special way
        ; so we don't go through the subterfuge we tried above.
    
        data = byte(data)
    
    endif else if n_elements(input) gt 0 then data = input

    ; Convert scalar to 1-d array.
    if siz[0] eq 0 and siz[1] ne 0 then data=(data)

    ; Do any scaling of the data.
    mwr_scale, data, scalval, offsetval, lscale=lscale, $
      iscale=iscale, bscale=bscale, null=null

    ; This may have changed the type.
    siz  = size(data)
    type = siz[siz[0]+1]


    ; If grouped data scale the group parameters.
    if keyword_set(group) then mwr_pscale, group, hdr, pscale=pscale, pzero=pzero

    if bof then begin
        chk_and_upd, hdr, 'SIMPLE', 'T','Primary Header created by MWRFITS v'+mwr_version()
        chk_and_upd, hdr, 'BITPIX', bitpixes[type]
        chk_and_upd, hdr, 'NAXIS', siz[0]
        chk_and_upd, hdr, 'EXTEND', 'T', 'Extensions may be present'
    endif else begin
        chk_and_upd, hdr, 'XTENSION', 'IMAGE','Image Extension created by MWRFITS v'+mwr_version()
        chk_and_upd, hdr, 'BITPIX', bitpixes[type]
        chk_and_upd, hdr, 'NAXIS', siz[0]
        chk_and_upd, hdr, 'PCOUNT', 0
        chk_and_upd, hdr, 'GCOUNT', 1
    endelse


    if keyword_set(group) then begin
        group_offset = 1
    endif else group_offset = 0

    if keyword_set(group) then begin
       chk_and_upd, hdr, 'NAXIS1', 0
    endif

    for i=1L, siz[0]-group_offset do begin
        chk_and_upd, hdr, 'NAXIS'+strcompress(string(i+group_offset),/remo), siz[i]
    endfor


    if keyword_set(group) then begin
        chk_and_upd, hdr, 'GROUPS', 'T'
        sizg = size(group)
        if sizg[0] ne 2 then begin
            print,'MWRFITS Error: Group data is not 2-d array'
           return
        endif
        if sizg[2] ne siz[siz[0]] then begin
            print,'MWRFITS Error: Group data has wrong number of rows'
           return
        endif
        chk_and_upd,hdr,  'PCOUNT', sizg[1]
        chk_and_upd, hdr, 'GCOUNT', siz[siz[0]]
    endif
    
    if n_elements(scalval) gt 0 then begin
    
        chk_and_upd, hdr, 'BSCALE', scalval
        chk_and_upd, hdr, 'BZERO', offsetval
    
    endif else begin
       
       ; Handle unsigned offsets
       bzero = mwr_unsigned_offset(type)
       if bzero gt 0 then begin
           chk_and_upd,hdr,'BSCALE', 1
           chk_and_upd, hdr, 'BZERO', bzero
           data = data + bzero
        endif
       
    endelse

    if keyword_set(group) then begin
        if keyword_set(pscale) then begin
            if n_elements(pscale) ne sizg[1] then begin
               print, 'MWRFITS Warning: wrong number of PSCALE values'
           endif else begin
                for i=1L, sizg[1] do begin
                    chk_and_upd, hdr, 'PSCALE'+strcompress(string(i),/remo)
               endfor
           endelse
        endif
        if keyword_set(pzero) then begin
            if n_elements(pscale) ne sizg[1] then begin
               print, 'MWRFITS Warning: Wrong number of PSCALE values'
           endif else begin
                for i=1L, sizg[1] do begin
                    chk_and_upd, hdr, 'PZERO'+strcompress(string(i),/remo)
               endfor
           endelse
        endif
    endif

    bytpix=abs(bitpixes[siz[siz[0]+1]])/8             ; Number of bytes per pixel.
    npixel = n_elements(data) + n_elements(group)     ; Number of pixels.

    if keyword_set(group) then mwr_groupinfix, data, group, hdr

    ; Write the FITS header
    mwr_header, lun, hdr

    ; This is all we need to do if input is undefined.
    if (n_elements(input) eq 0) or (siz[0] eq 0) then return

    ; Write the data.
    writeu, lun, data

    nbytes = bytpix*npixel
    filler = 2880 - nbytes mod 2880
    if filler eq 2880 then filler = 0
  
    ; Write any needed filler.
    if filler gt 0 then writeu, lun, replicate(0B,filler)
end


; Main routine -- see documentation at start
pro mwrfits, xinput, file, header,              $
        ascii=ascii,                            $
       separator=separator,                    $
       terminator=terminator,                  $
       create=create,                          $
       null=null,                              $
       group=group,                            $
       pscale=pscale, pzero=pzero,             $
       alias=alias,                            $
       use_colnum = use_colnum,                $
       lscale=lscale, iscale=iscale,              $
       bscale=bscale,                          $
       no_types=no_types,                      $
       silent=silent,                          $
       no_comment=no_comment,                  $
       logical_cols=logical_cols,              $
       bit_cols=bit_cols,                      $
       nbit_cols=nbit_cols,                    $
       status = status,                        $
       version=version


    ; Check required keywords.
    compile_opt idl2
    status = -1                     ;Status changes to 0 upon completion
    if (keyword_set(Version)) then begin
        print, "MWRFITS V"+mwr_version()+":  April 10, 2009"
    endif

    if n_elements(file) eq 0 then begin
        if (not keyword_set(Version)) then begin
            print, 'MWRFITS: Usage:'
            print, '    MWRFITS, struct_name, file, [header,] '
            print, '             /CREATE, /SILENT, /NO_TYPES, /NO_COMMENT, '
            print, '             GROUP=, PSCALE=, PZERO=,'
            print, '             LSCALE=, ISCALE=, BSCALE=,'
            print, '             LOGICAL_COLS=, BIT_COLS=, NBIT_COLS=,'
            print, '             ASCII=, SEPARATOR=, TERMINATOR=, NULL='
           print, '             /USE_COLNUM, ALIAS=, STATUS='
        endif
        return
    endif


    ; Save the data into an array/structure that we can modify.
 
    if n_elements(xinput) gt 0 then input = xinput

    on_ioerror, open_error

    ; Open the input file.
    ; If the create keyword is not specified we
    ; try to open the file readonly to see if it
    ; already exists and if so we append to it.
    ; An error implies the file does not exist.
    ;

    if  not keyword_set(create) then begin
        on_ioerror, not_found
        openr, lun, file, /get_lun,/swap_if_little
        free_lun, lun
        on_ioerror, open_error
        openu, lun, file, /get_lun, /append,/swap_if_little
        if not keyword_set(silent) then message,/inf,'Appending FITS extension to file ' + file
        bof = 0
        goto, finished_open
    endif

  not_found:
    on_ioerror, open_error
    openw, lun, file, /get_lun, /swap_if_little
    bof = 1
    on_ioerror, null

  finished_open:

    siz = size(input) 
     if siz[siz[0]+1] ne 8 then begin

        ; If input is not a structure then call image writing utilities.
        mwr_image, input, siz, lun, bof, header,    $
         null=null,                              $
         group=group,                            $
         pscale=pscale, pzero=pzero,             $
         lscale=lscale, iscale=iscale,              $
         bscale=bscale,                          $
         no_comment=no_comment,                  $
         silent=silent

    endif else if keyword_set(ascii) then begin

        if bof then mwr_dummy, lun
        ; Create an ASCII table.
        mwr_ascii, input, siz, lun, bof, header,     $
         ascii=ascii,                             $
         null=null,                               $
         use_colnum = use_colnum,                 $
         lscale=lscale, iscale=iscale,               $
         bscale=bscale,                           $
         no_types=no_types,                      $
         separator=separator,                     $
         terminator=terminator,                   $
         no_comment=no_comment,                   $
         alias=alias,                             $
         silent=silent

    endif else begin

        if bof then mwr_dummy, lun

        ; Create a binary table.
        mwr_tablehdr, lun, input, header, vtypes,    $
          no_types=no_types,                        $
          logical_cols = logical_cols,                    $
          bit_cols = bit_cols,                           $
          nbit_cols= nbit_cols,                     $
          alias=alias,                              $
          no_comment=no_comment,                    $
	  silent=silent
       
        mwr_tabledat, lun, input, header, vtypes

    endelse

    free_lun, lun
    status=0
    return
    
    ; Handle error in opening file.
  open_error:
    on_ioerror, null
    print, 'MWRFITS Error: Cannot open output: ', file
	 print,!ERROR_STATE.SYS_MSG
    if n_elements(lun) gt 0 then free_lun, lun
    
    return
end