mwrfits.pro

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

    chk_and_upd, header, 'GCOUNT', 1, 'Required value'
    chk_and_upd, header, 'TFIELDS', nfld, 'Number of columns in table'

    ;
    ; Handle the special cases.
    ;
    if keyword_set(logical_cols) then begin
       nl = n_elements(logical_cols)
       for i = 0, nl-1 do begin
           icol = logical_cols[i]
	   if types[icol-1] ne 'A' and types[icol-1] ne 'B' then begin
              print,'WARNING: Invalid attempt to create Logical column:',icol
                    goto, next_logical
           endif
           types[icol-1] = 'L'
  next_logical:
       endfor
    endif
       
    if keyword_set(bit_cols) then begin
       nb = n_elements(bit_cols)
       if nb ne n_elements(nbit_cols) then begin
           print,'WARNING: Bit_cols and Nbit_cols not same size'
           print,'         No bit columns generated.'
          goto, after_bits
       endif
       for i = 0, nb-1 do begin
           nbyte = (nbit_cols[i]+7)/8
           icol = bit_cols[i]
           if types[icol-1] ne 'B'  or (dims[icol-1] ne nbyte) then begin
              print,'WARNING: Invalid attempt to create bit column:',icol
                    goto, next_bit
           endif
           types[icol-1] = 'X'
           tdims[icol-1] = ''
           dims[icol-1] = nbit_cols[i]
  next_bit:
       endfor
  after_bits:
    endif



    ; Write scaling info as needed.
    if n_elements(offsets) gt 0 then begin
        w = where(offsets gt 0)

        for i=0, n_elements(w) - 1 do begin
            key = 'TSCAL'+strcompress(string(w[i])+1,/remo)
           chk_and_upd, header, key, 1
        endfor
    
        for i=0, n_elements(w) - 1 do begin
           key = 'TZERO'+strcompress(string(w[i]+1),/remo)
           chk_and_upd, header, key, offsets[w[i]]
        endfor
    
        if not keyword_set(no_comment) then begin
            key = 'TSCAL'+strcompress(string(w[0])+1,/remo)
           fxaddpar, header, 'COMMENT', ' ', before=key
           fxaddpar, header, 'COMMENT', ' *** Unsigned integer column scalings ***', before=key
           fxaddpar, header, 'COMMENT', ' ', before=key
        endif
    endif

    ; Now add in the TFORM keywords
    for i=0, nfld-1 do begin
       if dims[i] eq 1 then begin
           form = types[i]
       endif else begin
           form=strcompress(string(dims[i]),/remove) + types[i]
        endelse
       
       tfld = 'TFORM'+strcompress(string(i+1),/remove)
       
       ; Check to see if there is an existing value for this keyword.
       ; If it has the proper value we will not modify it.
       ; This can matter if there is optional information coded
       ; beyond required TFORM information.
              
       oval = fxpar(header, tfld)
       oval = strcompress(string(oval),/remove_all)
       if (oval eq '0')  or  (strmid(oval, 0, strlen(form)) ne form) then begin
           chk_and_upd, header, tfld, form
       endif
    endfor

    if (not keyword_set(no_comment)) then begin
        fxaddpar, header, 'COMMENT', ' ', before='TFORM1'
        fxaddpar, header, 'COMMENT', ' *** Column formats ***', before='TFORM1'
        fxaddpar, header, 'COMMENT', ' ', before='TFORM1'
    endif

    ; Now write TDIM info as needed.
    for i=nfld-1, 0,-1 do begin
        if tdims[i] ne '' then begin
            fxaddpar, header, 'TDIM'+strcompress(string(i+1),/remo), tdims[i],after=tfld
        endif
    endfor

    w=where(tdims ne '')
    if w[0] ne -1 and not keyword_set(no_comment) then begin
        fxaddpar, header, 'COMMENT', ' ', after=tfld
        fxaddpar, header, 'COMMENT', ' *** Column dimensions (2 D or greater) ***', after=tfld
        fxaddpar, header, 'COMMENT', ' ', after=tfld
    endif

    for i=0, nfld-1 do begin
        if tdims[i] ne '' then begin
            chk_and_upd, header, 'TDIM'+strcompress(string(i+1),/remo), tdims[i]
        endif
    endfor

    if n_elements(vtypes) gt 0 then begin
        fxaddpar, header, 'THEAP', nbyte*n_elements(input), 'Offset of start of heap'
        offset = 0L
        for i=0,n_elements(vtypes)-1 do begin
           if vtypes[i].status then offset = offset + vtypes[i].offset*vtypes[i].ilen
        endfor
        fxaddpar, header, 'PCOUNT', offset, 'Size of heap'
    endif

    ;
    ; Last add in the TTYPE keywords if desired.
    ;
    if not no_types then begin
       for i=0, nfld - 1 do begin
           key = 'TTYPE'+strcompress(string(i+1),/remove)
           if not keyword_set(use_colnums) then begin
               value= mwr_checktype(tags[i],alias=alias)+' '
           endif else begin
               value = 'C'+strmid(key,5,2)
           endelse
           chk_and_upd, header, key, value
       endfor
       
        if (not keyword_set(no_comment)) then begin
           fxaddpar, header, 'COMMENT', ' ', before='TTYPE1'
           fxaddpar, header, 'COMMENT', ' *** Column names *** ',before='TTYPE1'
           fxaddpar, header, 'COMMENT', ' ',before='TTYPE1'
       endif
    endif

    if (not keyword_set(no_comment)) then begin
        fxaddpar, header, 'COMMENT', ' ', after='TFIELDS'
        fxaddpar, header, 'COMMENT', ' *** End of mandatory fields ***', after='TFIELDS'
        fxaddpar, header, 'COMMENT', ' ', after='TFIELDS'
    endif

    ; Write to the output device.
    mwr_header, lun, header

end

; Modify the structure to put the pointer column in.
function mwr_retable, input, vtypes


    offset = 0L
    tags = tag_names(input);
;Create an output structure identical to the input structure but with pointers replaced
; by a 2 word lonarr to point to the heap area

      if vtypes[0].status then begin
        output = CREATE_STRUCT(tags[0],lonarr(2))
      endif else begin
         output = CREATE_STRUCT(tags[0],input[0].(0))
      endelse
      for i=1, n_elements(tags) -1 do begin
         if vtypes[i].status then begin
           output = CREATE_STRUCT(temporary(output), tags[i], lonarr(2))
         endif else begin
           output = CREATE_STRUCT(temporary(output), tags[i], input[0].(i))
         endelse
      endfor
      output = replicate(temporary(output), N_elements(input) )
      struct_assign, input, output      ;Available since V5.1

    for i=0, n_elements(tags)-1 do begin
       if vtypes[i].status then begin
           for j=0, n_elements(input)-1 do begin
              ptr = input[j].(i)
              if ptr_valid(ptr) then begin
                  sz = size(*ptr)
                  if sz[0] eq 0 then xsz = 1 else xsz= sz[1]

                  output[j].(i)[0] = xsz
                  output[j].(i)[1] = offset
                  
                  offset = offset + vtypes[i].ilen*xsz
              endif
           endfor
       endif
    endfor
    return,output
end

; Write the heap data.
function mwr_writeheap, lun, vtypes

    offset = 0L
    
    for i=0, n_elements(vtypes)-1 do begin
       if vtypes[i].status then begin
           
           itype = vtypes[i].itype
           unsigned = mwr_unsigned_offset(itype)
           
           ptrs = vtypes[i].data
           
           for j=0,n_elements(ptrs)-1 do begin
              if ptr_valid(ptrs[j]) then begin
                  if (unsigned gt 0) then begin
                     *ptrs[j] = *ptrs[j] + unsigned
                  endif

                      writeu, lun, *ptrs[j]
                  
                  sz = size(*ptrs[j])
                  xsz = 1 > sz[1]
                  offset = offset + xsz * vtypes[i].ilen
              endif
           endfor
       endif
    endfor

    return, offset
    
end

; Write the brinary table.
pro mwr_tabledat, lun, input, header, vtypes
    ;
    ; file              -- unit to which data is to be written.
    ; Input              -- IDL structure
    ; Header       -- Filled header

    nfld = n_tags(input)

    ; Any special processing?

    typ = intarr(nfld)
    for i=0, nfld-1 do begin
        
        typ[i] = size(input.(i),/type)
	    if (typ[i] eq 7) then begin

             dim = size(input.(i),/dimen) >1
             siz = max(strlen(input.(i))) > 1
	     input.(i) = $
	        strmid( input.(i) + string(replicate(32b, siz)), 0, siz)

       endif
 
       unsigned = mwr_unsigned_offset(typ[i])
       if (unsigned gt 0) then begin
           input.(i) = input.(i) + unsigned
       endif
       
    endfor

    if n_elements(vtypes) gt 0 then begin
          
      
        input = mwr_retable(input, vtypes)
    endif

    ; Write the data segment.
    ;
    writeu, lun, input

    nbyte = long(fxpar(header, 'NAXIS1'))
    nrow  = n_elements(input)

    heap = 0
    if n_elements(vtypes) gt 0 then begin
        heap = mwr_writeheap(lun, vtypes)
    endif

    siz   = nbyte*nrow + heap
    padding = 2880 - (siz mod 2880)
    if padding eq 2880 then padding = 0

    ;
    ; If necessary write the padding.
    ;
    if padding gt 0 then begin
        pad = bytarr(padding)  ; Should be null-filled by default.
        writeu, lun, pad
    endif

end


; Scale parameters for GROUPed data.
pro mwr_pscale, grp, header, pscale=pscale, pzero=pzero

; This function assumes group is a 2-d array.

    if not keyword_set(pscale) and not keyword_set(pzero) then return

    if not keyword_set(pscale) then begin
        pscale = dblarr(sizg[1])
        pscale[*] = 1.
    endif

    w = where(pzero eq 0.d0)

    if w[0] ne 0 then begin
        print, 'MWRFITS  Warning: PSCALE value of 0 found, set to 1.'
        pscale[w] = 1.d0
    endif

    if keyword_set(pscale) then begin
        for i=0L, sizg[1]-1 do begin
            key= 'PSCAL' + strcompress(string(i+1),/remo)
            chk_and_upd, header, key, pscale[i]
        endfor
    endif

    if not keyword_set(pzero) then begin
        pzero = dblarr(sizg[1])
        pzero[*] = 0.
    endif else begin
        for i=0L, sizg[1]-1 do begin
            key= 'PZERO' + strcompress(string(i+1),/remo)
            chk_and_upd, header, key, pscale[i]
        endfor
    endelse

    for i=0L, sizg[1]-1 do begin
        grp[i,*] = grp[i,*]/pscale[i] - pzero[i]
    endfor

end


; Find the appropriate scaling parameters.
pro mwr_findscale, flag, array, nbits, scale, offset, error


    error = 0
    if n_elements(flag) eq 2 then begin
         scale  = double(flag[0])
        offset = double(flag[1])
    endif else if n_elements(flag) eq 1 and flag[0] ne 1 then begin
         minmum = min(array, max=maxmum)
        offset = 0.d0
        scale  = double(flag[0])
    endif else if n_elements(flag) ne 1 then begin
         print, 'MWRFITS Error: Invalid scaling parameters.'
        error  = 1
        return
    endif else begin
        
         minmum = min(array, max=maxmum)
        scale  = (maxmum-minmum)/(2.d0^nbits)
        amin   = -(2.d0^(nbits-1))
        if (amin gt -130) then amin = 0  ; looking for -128
        offset = minmum - scale*amin
        
    endelse
    return
end

; Scale and possibly convert array according to information
; in flags.
pro mwr_scale, array, scale, offset, lscale=lscale, iscale=iscale,  $
   bscale=bscale, null=null


    ; First deallocate scale and offset
    if n_elements(scale)  gt 0 then xx = temporary(scale)
    if n_elements(offset) gt 0 then xx = temporary(offset)

    if not keyword_set(lscale) and not keyword_set(iscale) and  $
       not keyword_set(bscale) then return

    siz = size(array)
    if keyword_set(lscale) then begin

        ; Doesn't make sense to scale data that can be stored exactly.
        if siz[siz[0]+1] lt 4 then return
        amin = -2.d0^31
        amax = -(amin + 1)
    
        mwr_findscale, lscale, array, 32, scale, offset, error

    endif else if keyword_set(iscale) then begin
        if siz[siz[0]+1] lt 3 then return
        amin = -2.d0^15
        amax = -(amin + 1)
    
        mwr_findscale, iscale, array, 16, scale, offset, error

    endif else begin
        if siz[siz[0]+1] lt 2 then return
    
        amin = 0
        amax = 255
    
        mwr_findscale, bscale, array, 8, scale, offset, error
    endelse

    ; Check that there was no error in mwr_findscale
    if error gt 0 then return

    if scale le 0.d0 then begin
        print, 'MWRFITS Error: BSCALE/TSCAL=0'
        return
    endif

    array = round((array-offset)/scale)

    w=where(array lt 0)
    w = where(array gt amax)
    if w[0] ne -1 then begin
        if keyword_set(null) then array[w] = null else array[w]=amax
    endif

    w = where(array lt amin)
    if w[0] ne -1 then begin
        if keyword_set(null) then array[w] = null else array[w] = amin
    endif