fits_open.pro

basic median filter simulation

;
        n = 100
        xtension = strarr(n)
        extname = strarr(n)
        extver = lonarr(n)
        extlevel = lonarr(n)
        gcount = lonarr(n)
        pcount = lonarr(n)
        bitpix = lonarr(n)
        naxis  = lonarr(n)
        axis = lonarr(20,n)
        start_header = lon64arr(n)        ; starting byte in file for header
        start_data = lon64arr(n)          ; starting byte in file for data
        position = long64(0)             ; current byte position in file
        skip = long64(0)                 ; Amount to skip from current position
;
; read and process each header in the file if open for read or update
;
        extend_number = 0               ; current extension number being
                                        ; processed
 
        if open_for_read or open_for_update then begin
            main_header = 1             ; first header in file flag
            h = bytarr(80,36,/nozero)   ; read buffer
;
; loop on headers in the file
;
            repeat begin
              if skip GT 0 then point_lun,unit,position 
              start = position
;
; loop on header blocks
;
                first_block = 1         ; first block in header flag
                repeat begin

                    if position+2879 ge nbytes_in_file then begin
                        if extend_number eq 0 then begin
                                message = 'EOF encountered while reading header'
                                goto,error_exit
                        endif
                        print,'EOF encountered reading extension header'
                        print,'Only '+strtrim(extend_number-1,2) + $
                                ' extensions processed'
                        goto,done_headers
                    endif

                    readu,unit,h
                    position = position + 2880
                    hdr = string(h>32b)
                    endline = where(strmid(hdr,0,8) eq 'END     ',nend)
                    if nend gt 0 then hdr = hdr[0:endline[0]]
                    if first_block then begin
;
; check for valid header (SIMPLE keyword must be first for PDU and
; XTENSION keyword for the extensions.
;
                        header = hdr 
                        keyword = strmid(header[0],0,8)
                        if (extend_number eq 0) and $
                           (keyword ne 'SIMPLE  ') then begin
                                message = 'Invalid header, no SIMPLE keyword'
                                goto,error_exit
                        endif

                        if (extend_number gt 0) and $
                           (keyword ne 'XTENSION') then begin
                                print,'Invalid extension header encountered'
                                print,'XTENSION keyword missing'
                                print,'Only '+strtrim(extend_number-1,2) + $
                                        ' extensions processed'
                                goto,done_headers
                        endif

                    end else header = [header,hdr]

                    first_block = 0
                end until (nend gt 0)   

;
; print header if hprint set
;
                if keyword_set(hprint) then hprint,header
;
; end of loop on header blocks
;
; Increase size of vectors if needed
;
                if extend_number ge n then begin
                        xtension = [xtension,strarr(n)]
                        extname = [extname,strarr(n)]
                        extver = [extver,lonarr(n)]
                        extlevel = [extver,lonarr(n)]
                        gcount = [gcount,lonarr(n)]
                        pcount = [pcount,lonarr(n)]
                        bitpix = [bitpix,lonarr(n)]
                        naxis  = [naxis,lonarr(n)]
                        old_axis = axis
                        axis = lonarr(20,n*2)
                        axis[0,0] = old_axis
                        start_header = [start_header,lonarr(n)]
                        start_data = [start_data,lonarr(n)]
                        n = n*2
                end
;
; extract information from header
;
                xtension[extend_number] = strtrim(sxpar(header,'xtension'))
                st = sxpar(header,'extname', Count = N_extname)
                if N_extname EQ 0 then st = ''
                extname[extend_number] = strtrim(st,2)  
                extver[extend_number] = sxpar(header,'extver')          
                extlevel[extend_number] = sxpar(header,'extlevel')              
                gcount[extend_number] = sxpar(header,'gcount')
                pcount[extend_number] = sxpar(header,'pcount')
                bitpix[extend_number] = sxpar(header,'bitpix')
                nax = sxpar(header,'naxis')
                naxis[extend_number] = nax
                if nax gt 0 then for i=1,nax do $
                    axis[i-1,extend_number] = sxpar(header,'naxis'+strtrim(i,2))
                start_data[extend_number] = position    
                start_header[extend_number] = start
;
; if first header, save without FITS required keywords
;
                if extend_number eq 0 then begin
                    hmain = header
                    random_groups = sxpar(header,'groups')
                    sxdelpar,hmain,['SIMPLE','BITPIX','NAXIS','NAXIS1', $
                                    'NAXIS2','NAXIS3','NAXIS4','NAXIS5', $
                                    'NAXIS6','NAXIS7','NAXIS8','EXTEND', $
                                    'PCOUNT','GCOUNT','GROUPS','BSCALE', $
                                    'BZERO','NPIX1','NPIX2','PIXVALUE']
                        if (pcount[0] gt 0) then for i=1,pcount[0] do $
                        sxdelpar,hmain,['ptype','pscal','pzero']+strtrim(i,2)
                endif
;
; skip past data to go to next header
;
                ndata = long64(axis[0,extend_number])
                if naxis[extend_number] gt 1 then $
                                for i=2,naxis[extend_number] do $
                                    ndata = ndata*axis[i-1,extend_number]
                nbytes = (abs(bitpix[extend_number])/8) * $
                       (gcount[extend_number]>1)*(pcount[extend_number] + ndata)
                skip = (nbytes + 2879)/2880*2880
                position = position + skip

;
; end loop on headers
;           

                extend_number = extend_number + 1
            end until (position ge nbytes_in_file-2879)
        end
;
; point at end of file in /extend
;
done_headers:
        if open_for_update then point_lun,unit,nbytes_in_file
;
; number of extensions
;
        if open_for_write then nextend = -1 $
                          else nextend = extend_number - 1
;
; set up blank hmain if open for write
;
        if open_for_write then begin
                hmain = strarr(1)
                hmain[0] = 'END     '
        end
;
; create output structure for the file control block
;
        if open_for_write or open_for_update then begin
                fcb = {filename:fname,unit:unit,nextend:nextend, $
                        open_for_write:open_for_write + open_for_update*2}
           end else begin
                nx = nextend
                fcb = {filename:fname,unit:unit,nextend:nextend, $
                        xtension:xtension[0:nx],extname:extname[0:nx], $
                        extver:extver[0:nx],extlevel:extlevel[0:nx], $
                        gcount:gcount[0:nx],pcount:pcount[0:nx], $
                        bitpix:bitpix[0:nx],naxis:naxis[0:nx], $
                        axis:axis[*,0:nx], $
                        start_header:start_header[0:nx], $
                        start_data:start_data[0:nx],hmain:hmain, $
                        open_for_write:open_for_overwrite*3,$
                        last_extension:-1, $
                        random_groups:random_groups, $
                        nbytes: nbytes_in_file }
        end
        !err = 1            ;For obsolete users still using !err
        return
;
; error exit
;
ioerror: 
        message = !ERROR_STATE.msg
error_exit:
        free_lun,unit
        !err = -1
        if keyword_set(no_abort) then return
        message,' ERROR: '+message,/CON
        return
end