st_diskread.pro

basic median filter simulation

pro st_diskread, infiles, DUMP = dump
;+
; NAME: 
;       ST_DISKREAD
;
; PURPOSE:  
;       Read HST FITS formatted disk files and reconstruct GEIS (STSDAS) files.
;
; CALLING SEQUENCE:  
;       ST_DISKREAD, infiles
;
; INPUT PARAMETER:
;       infiles - (scalar string) input disk files to be converted into GEIS
;                       files. Wildcards are allowed.
; FILES CREATED:
;
;   GEIS files:
;         The GEIS file is reconstructed from each input Fits file. The 
;       output filename is composed from the rootname of the observation
;       and the appropriate GEIS file extension (i.e. d0h/d, c0h/d, etc.).
;   Tables:
;         If input file is a fits table, the output is an SDAS table.
;
; EXAMPLES:
;       a) Reconstruct the GEIS file for disk FITS file z29i020ct*.fits.
;               st_diskread,'z29i020ct*.fits'
;
; PROCEDURES CALLED:
;       ST_DISK_DATA, ST_DISK_TABLE, ST_DISK_GEIS
;       FTSIZE,SXPAR(),TAB_CREATE, TAB_WRITE
; HISTORY: 
;       10/17/94        JKF/ACC - taken from ST_TAPEREAD.
;       11/02/94        JKF/ACC - added /block on open statement to
;                                 handle files with 512 bytes/record.
;       12/6/95         JKF/ACC - include new jitter files...replaces
;                                               st_read_jitter.pro.
;       03/5/96         W. Landsman, change FORRD to READU, remove Version 1
;                               type codes, add message facility
;       05/20/00        W. Landsman, remove obsolete !ERR calls, new calling
;                               sequence to FTINFO
;       09/2006        W. Landsman, remove obsolete keywords to OPEN
;
;****************************************************************************
;       Converted to IDL V5.0   W. Landsman   September 1997
;-

 On_error,2

 if n_params() lt 1 then begin
        print,'Syntax - ST_DISKREAD, infiles'
        return
 endif
 !ERROR = 0
 if not keyword_set(DUMP) then dump = 0
;
; Search for names of input disk FITS files.
;
   file_list = file_search(infiles,count=count)
   if count le 0 then $                                              
        message,' No files found: '+ infiles $
   else message,/INF, $
        'Number of files to process: ' + strtrim(count,2)
;
; Loop on files
;
   for file = 0,count-1 do begin
        openr,unit,file_list[file],/get_lun
;
; read data header and data
;
        st_disk_data,unit,h,data,fname,gcount,dimen,opsize,nbytes,itype
        if !ERROR NE 0 then return
;
; read optional table extension
;
        st_disk_table,unit,htab,tab,table_available
        if !ERROR NE 0 then return
;
; Finished reading the input dataset at this point. Now process the information
; and create the output datasets.
;
;       GEIS file or trailer text file
;

        if sxpar(h,'naxis') gt 0 then begin
                st_disk_geis,h,data,htab,tab,table_available, $
                        fname,gcount,dimen,opsize,nbytes,itype  ;GEIS file
                if !ERROR NE 0 then return
                if dump gt 0 then $
                        print,format='(t5,i4,t15,a)',file+1,strlowcase(fname)
        end else begin                  ;either a text trailer or jitter table

           outname = strtrim(sxpar(htab,'extname'),2)
           if outname eq strtrim(0,2) then $
                outname= strtrim(sxpar(h,'filename')) 

           if  table_available then begin               

                outname = strtrim(sxpar(htab,'extname'))
                s=size(tab) & nl=s[2]                           
                name=strtrim(sxpar(htab,'extname'))             ;file name
                ;
                ;  What type of table?
                ;     - trailer file - ascii table
                ;     - jitter data  - sdas table
                ;
                if strpos(strlowcase(name),'jit') eq -1 then begin; text trailer
                  ;
                  ;     Special case NAME: PODPS/IRAF uses j7 as special 
                  ;     character, so that a file with z0j7<...> will be 
                  ;     created as z0.<...> ( . is substituted for j7 ).
                  ;     To avoid: Check file name for ., if found replace
                  ;     with j7.
                  ;
                  invalid_char = strpos(name,'.')
                  if invalid_char lt 5 then begin
                        message,' Warning: Invalid filename found: '+name ,/cont
                        name = strmid(name,0,invalid_char) + 'j7' + $
                                 strmid(name,invalid_char+1,strlen(name)) 
                        message,'   Filename will be changed to: '+ name,/cont
                  end       
                 
                  openw,ounit,name,/get_lun
                  for i = 0,nl-1 do printf,ounit,strtrim(string(tab[*,i]))
                  free_lun,ounit
                  if dump gt 0 then $
                        print,format='(t5,i4,t15,a)',file+1,strlowcase(name)
                end else begin                                  ; jitter table
                  ;
                  ; Convert from FITS to SDAS table
                  ;
                  ftsize,htab,tab,ncols,nrows,tfields
                  tab_create,tcb,otab,tfields,nrows,ncols/2
                  ftinfo,htab,ft_str
                  fname = ft_str.ttype
                  for j= 0, tfields-1 do begin
                        val=ftget(ft_str,tab,j+1)     ; extract column
                        tab_put,strtrim(fname[i]),val,tcb,otab
                  end
                  tab_write,outname,tcb,otab,htab
                  if dump gt 0 then $
                        print,format='(t5,i4,t15,a,a)',file+1, $
                                strlowcase(outname)," jitter table "
                end
           end else $
                if dump gt 0 then $
                        print,format='(t5,i4,t15,a,a)',file+1, $
                                strlowcase(outname)," (No data found)
        end     
        free_lun,unit
   endfor
return
end
;
pro st_disk_data,unit,h,data,name,gcount,dimen,opsize,nbytes,itype
;**************************************************************************
;+
; NAME:
;       ST_DISK_DATA 
;
; PURPOSE:
;       Routine to read next header and data array from an HST FITS disk file.
;       This is a subroutine of ST_DISKREAD and not intended for stand alone 
;       use.
;
;CALLING SEQUENCE:
;       st_disk_data,unit,h,data,name,gcount,dimen,opsize,nbytes,itype
;
;INPUTS:
;       unit - logical unit number.
;
;OUTPUTS:
;       h - FITS header
;       data - data array
;       name - file name
;       gcount - number of groups
;       dimen - data dimensions
;       opsize - parameter blocks size
;       nbytes - bytes per data group
;       itype - idl data type
;
; Notes:
;       This is not a standalone program. Use ST_DISKREAD.
;
; PROCEDURES CALLED:
;       GETTOK(), SXPAR()
; HISTORY:
;       10/17/94        JKF/ACC         - taken from ST_TAPE_DATA.
;
;***************************************************************************
;-
        On_error,2
;
; read fits header
;
        h = strarr(500)
        nhead = 0
        while 1 do begin
            buf=bytarr(2880)
            readu,unit,buf
        
            for i=0,35 do begin
                st = string(buf[i*80:i*80+79])
                h[nhead]=st
                if strtrim(strmid(st,0,8)) eq 'END' then goto,fini
                nhead=nhead+1
            endfor
        endwhile
fini:
;
; get keywords from header needed to read data
;
        bitpix = sxpar(h,'bitpix', Count = N_bitpix)

        if N_bitpix EQ 0 then begin
            message,/CON,'ERROR - BITPIX missing from FITS header'
            return
        endif

        naxis = sxpar(h,'naxis', Count = N_naxis)
        if N_naxis EQ 0 then begin
            message,/CON,'ERROR- NAXIS missing from FITS header'
            return
        endif
        if naxis eq 0 then return               ;NO data to read
;
; get scale factors
;
        bscale = sxpar(h,'bscale', Count = N_bscale)
        if N_bscale EQ 0 then bscale=1.
        bzero = sxpar(h,'bzero', Count = N_bzero)
        if N_bzero EQ 0 then bzero=0.
        iraf_bp = sxpar(h,'IRAF-B/P')           ;Geis file bitpix
        if iraf_bp ne 64 then begin
                bscale = float(bscale)
                bzero = float(bzero)
            end else begin
                bscale = double(bscale)
                bzero = double(bzero)
        end
;
; determine output bitpix
;
        obitpix = abs(bitpix)
        if (bscale ne 1.0) or (bzero ne 0.0) then obitpix = 32
        if iraf_bp eq 64 then obitpix = 64 
;
; get dimensions
;
        dimen = lonarr(naxis)
        npoints = 1L
        for i=0,naxis-1 do begin
            dimen[i]=sxpar(h,'naxis'+strtrim(i+1,2))
            if dimen[i] le 0 then begin
                message,/CON,'ERROR- Invalid data dimension'
                return
            endif
            npoints = npoints*dimen[i]
        endfor
;
; determine group count
;
        gcount = sxpar(h,'sdasmgnu')>1
        if gcount gt 1 then begin
                naxis = naxis-1
                dimen = dimen[0:naxis-1]     
                if n_elements(dimen) eq 1 then dimen = lonarr(1)+dimen
                npoints = npoints/gcount
        endif
;
; determine orignal psize in bytes
;
        opsize = sxpar(h,'opsize', Count = N_opsize)
        if N_opsize EQ 0 then opsize = 0
        opsize = opsize/8
;
; set up data array
;
        case bitpix of
           8: data = make_array(dimen=dimen,/byte)
          16: data = make_array(dimen=dimen,/int)
          32: data = make_array(dimen=dimen,/long)
          64: data = make_array(dimen=dimen,/double)
         -32: data = make_array(dimen=dimen,/float)
         -64: data = make_array(dimen=dimen,/double)

          else: begin
                message,/CON,'ERROR - Invalid BITPIX value'
                return
                end
        endcase
;
; determine file name
;
        ;
        ; Keyword IRAFNAME has been changed to FILENAME in new style 
        ;       PODPS keywords (JHB 11-2-91)
        ;
        name = sxpar(h,'FILENAME', Count = N_filename)
        if N_filename EQ 0 then begin
                name = sxpar(h,'IRAFNAME', Count = N_irafname)
                if N_irafname EQ 0 then $
                        message,' Keyword(IRAFNAME) missing from data header'+ $
                        '...ABORTING '
        endif

        ;
        ; Special case NAME: PODPS/IRAF uses j7 as special
        ; character, so that a file with z0j7<...> will be
        ; created as z0.<...> ( . is substituted for j7 ).
        ; To avoid: Check file name for ., if found replace
        ; with j7.
        ; Special case code added by JKF/ACC 12/30/91
        ;
        invalid_char = strpos(name,'.')
        if invalid_char lt 5 then begin
            message,' Warning: Invalid filename found: '+name ,/cont
            name = strmid(name,0,invalid_char) + 'j7' + $
                    strmid(name,invalid_char+1,strlen(name))
            message,'   Filename will be changed to: '+ name,/cont
            end

        name = strtrim(gettok(name,'.') +'.'+ gettok(name,'.'),2)
        pos = strpos(name,'_cvt')               ;take out _cvt
        if pos gt 4 then name = strmid(name,0,pos) + $
                                strmid(name,pos+4,strlen(name)-pos-4)
        dname = name
        strput,dname,'d',strlen(name)-1 ;change last character to a d
;
; determine number of blocks in the file
;
        bytes_per_point = obitpix/8
        in_bytes_per_point = abs(bitpix)/8
        nbytes = bytes_per_point * npoints
        nblocks = ((nbytes + opsize)*gcount + 511)/512
;
; open output data file
;
        close,1
        openw,1,dname
;
; create output assoc variable
;
        if (bzero eq 0) and (bscale eq 1) and (bitpix gt 0) then begin
                s = size(data) & itype = s[s[0]+1] ; idl data type
                tmp_data = make_array( dimen=dimen, type= itype )

           end else begin   

                if obitpix eq 32 then begin
                        tmp_data =  make_array(dimen=dimen,/float)
                        itype = 4
                   end else begin
                        tmp_data =  make_array(dimen=dimen,/double)
                        itype = 5
                end
        end 
;
; read data
;

        pointer = 2880          ;byte pointer in current 2880 byte disk record
              
        for group=0,gcount-1 do begin           ;loop on groups
            pos = 0                             ;current pointer in data array
            while pos lt npoints do begin
                if pointer ge 2880 then begin
                   readu,unit,buf
                   case bitpix of
                        16: byteorder,buf,/NtoHS
                        32: byteorder,buf,/NtoHL
                        -32: byteorder,buf,/XDRTOF
                        -64: byteorder,buf,/XDRTOD
                        ELSE:
                   endcase
                   pointer = 0
                endif
                words_needed = (npoints-pos)
                bytes_needed = words_needed*in_bytes_per_point
                bytes_to_take = (2880-pointer) < bytes_needed
                words_to_take = bytes_to_take/in_bytes_per_point

                case bitpix of
                        8: data[pos]=buf[pointer:bytes_to_take-1]
                        16: data[pos]=fix(buf,pointer,words_to_take)
                        32: data[pos]=long(buf,pointer,words_to_take)
                        64: data[pos]=double(buf,pointer,words_to_take)
                       -32: data[pos]=float(buf,pointer,words_to_take)   ;IEEE
                       -64: data[pos]=double(buf,pointer,words_to_take)  ;IEEE
                endcase
                pos = pos + words_to_take