; Filename:  read_sws2a.pro
;
; Usage:
;
;   To run this program, use the following command:
;
;   IDL> read_sws2a
;
; Description:
;  
;   This file contains one (1) IDL procedure and four (4) functions
;   used to read the SeaWinds Level 2A data in Hierarchical Data 
;   Format (HDF).  The routines are as follows.
;   
;   1. correct_uint8: a function to correct variables 
;                     stored as 8-bit unsigned integers.
;
;   2. correct_uint16: a function to correct variables 
;                      stored as 16-bit unsigned integers.
;
;   3. get_vdata: a function to read the timetags which 
;                 are stored as VDATA.
;
;   4. get_sds: a function to read the data which are 
;               stored as Scientific Data Sets.
;
;   5. read_sws2a: the main procedure. Calls get_vdata and 
;                    get_sds.
;
; Notes:
;
; 1. The directory on your local system which contains the SeaWinds
;    L2A data must be input in program read_sws2a.pro.  (Search on
;    "data_dir" to find this line quickly.)
;
; 2. Due to the size of the SDSs contained in these files, the L2A 
;    data are read in slabs.  If you would like to read larger slabs,
;    change variable "slab_size".  Setting slab_size to 1702 will allow
;    *most* L2A SDSs to be read in their entirety.
;
; 3. The correction of Unsigned Integers is currently hard-coded for
;    versions of IDL prior to 5.2.
;
; 4. Please send all comments and questions concerning these routines
;    to sws@podaac.jpl.nasa.gov.
;
;06/07/2002, Ted Lungu
;
; based on QuikSCAT version by K.L. Perry, C.S. Hsu, R.S. Dunbar
;
; Modifications:
;
y
;
;==================================================================
;   Copyright (c) 2002, California Institute of Technology
;==================================================================

;==================================================================
; correct_uint8: a function to correct data which are stored
;                as 8-bit unsigned integers for versions of
;                IDL prior to 5.2
;==================================================================
function correct_uint8,sds_data
    w=where(sds_data lt 0,cnt)
    if (cnt gt 0) then begin
      sds_data(w)=sds_data(w)+256
    endif
return,sds_data
end

;==================================================================
; correct_uint16: a function to correct data which are stored
;                 as 16-bit unsigned integers for versions of
;                 IDL prior to 5.2
;==================================================================
function correct_uint16,sds_data
    w=where(sds_data lt 0,cnt)
    if (cnt gt 0) then begin
      sds_data(w)=sds_data(w)+65536.
      w2=where(sds_data lt 0,cnt2)
    endif
return,sds_data
end

;==================================================================
; get_vdata: a function to read the timetags which are stored as
;            VDATA.
;==================================================================
function get_vdata,fid,vdata_name
  result=HDF_VD_FIND(fid,vdata_name)
  vid=HDF_VD_ATTACH(fid,result)
  result=HDF_VD_READ(vid,vdata)
  HDF_VD_DETACH,vid
return,vdata
end

;==================================================================
; get_sds: a function to read the data which are stored
;          as Scientific Data Sets.  This function also
;          multiplies the data by the calibration and 
;          subtracts the offset.
;==================================================================
function get_sds,sd_id,sds_name,slab_start,slab_size

    index=HDF_SD_NAMETOINDEX(sd_id,sds_name)
    sds_id=HDF_SD_SELECT(sd_id,index)
    HDF_SD_GETINFO,sds_id,ndims=ndims,dims=dims,caldata=cal,type=data_type

; make sure edge, start and stride have correct array sizes
    edge=intarr(ndims)
    start=intarr(ndims)
    stride=intarr(ndims)

    for i=0,ndims-1 do begin
      edge(i)=dims(i)
      start(i)=0
      stride(i)=1
    endfor
    edge(ndims-1)=slab_size
    start(ndims-1)=slab_start

    HDF_SD_GETDATA,sds_id,data,stride=stride,start=start,count=edge

;Correct Unsigned Integers
;;  note: Versions of IDL prior to 5.2 do not handle unsigned 
;;        integers properly.  These versions also do not identify 
;;        DFNT_UINT's using HDF_SD_GETINFO.  Therefore, the 
;;        following hard code will be included until IDL 5.2 is 
;;        "standard". --KLP.
;;

;; UINT8
    if ((sds_name eq "num_sigma0_per_cell") or $
        (sds_name eq "num_wvc_tb_in") or $
        (sds_name eq "num_wvc_tb_out") or $
        (sds_name eq "cell_index")) then begin
      data=correct_uint8(float(data))
    endif

;; UINT16
    if ((sds_name eq "mean_wvc_tb_in") or $
	(sds_name eq "mean_wvc_tb_out") or $
	(sds_name eq "std_dev_wvc_tb_in") or $
	(sds_name eq "std_dev_wvc_tb_out") or $
        (sds_name eq "cell_lon") or $
        (sds_name eq "cell_azimuth") or $
        (sds_name eq "kp_beta") or $
        (sds_name eq "sigma0_qual_flag") or $
        (sds_name eq "sigma0_mode_flag") or $
        (sds_name eq "surface_flag")) then begin
      data=correct_uint16(float(data))
    endif

; Apply the scale and offset
    rdata = data*cal.cal - cal.offset
    HDF_SD_ENDACCESS,sds_id

return,rdata
end

;==================================================================
; read_sws2a: the main procedure in this file.  Calls are made 
;               to get_vdata to read the timetags and get_sds to 
;               read the data.  The results are then printed to 
;               the screen.
;==================================================================

pro read_sws2a

;Select a Level 2A file
;***** Change data_dir to suit your local system
  data_dir="./"
  filename=pickfile(/READ,path=data_dir,filter='SW_S2A*', $
                    title='SeaWinds L2A')

  print,' '
  print,'FILENAME: ',filename
  print,' '

;Read the Time Tags contained in the VDATA
  fid=HDF_OPEN(filename,/READ)
  wvc_row_time= get_vdata(fid,'wvc_row_time')
  HDF_CLOSE,fid

;Select wind vector cell rows to be read
  read,'Enter the first and last wvc row rec numbers to be read: ',irec1,irec2
  if (irec1 gt irec2) then begin
    itmp=irec1
    irec1=irec2
    irec2=itmp
  endif
  if ((irec1 lt 1) or (irec2 gt 1702)) then begin
    print,'ERROR: wvc rows must be between 1 and 1702'
    stop
  endif

;Read the Scientific Data Sets
  sd_id=HDF_SD_START(filename,/READ)

;; The L2A SDSs are read in slabs of size 1.  An example of
;; reading the SDSs in their entirety is shown below.
;;
;;  index=HDF_SD_NAMETOINDEX(sd_id,'row_number')
;;  sds_id=HDF_SD_SELECT(sd_id,index)
;;  HDF_SD_GETINFO,sds_id,ndims=rank,dims=dims,label=name, $
;;                  type=data_type,caldata=cal
;;  slab_size=strtrim(dims(0),2)   
;;  row_number=get_sds(sd_id,'row_number',ir,slab_size)
;;  for ir=irec1-1,irec2-1 do begin
;;    print,'TIME: ',string(wvc_row_time(*,ir))
;;    print,'WVC Row Number: ',row_number(ir)
;;  endfor

  slab_size=1

; Subract 1 from irec1 and irec2 to adjust for IDL running 
; from 0 instead of 1 

  for ir=irec1-1,irec2-1 do begin
    row_number= get_sds(sd_id,'row_number',ir,slab_size)
    num_sigma0= get_sds(sd_id,'num_sigma0',ir,slab_size)
    num_sigma0_per_cell=get_sds(sd_id,'num_sigma0_per_cell',ir,slab_size) 
    num_wvc_tb_in=get_sds(sd_id,'num_wvc_tb_in',ir,slab_size)
    num_wvc_tb_out=get_sds(sd_id,'num_wvc_tb_out',ir,slab_size)
    mean_wvc_tb_in=get_sds(sd_id,'mean_wvc_tb_in',ir,slab_size)
    mean_wvc_tb_out=get_sds(sd_id,'mean_wvc_tb_out',ir,slab_size)
    std_dev_wvc_tb_in=get_sds(sd_id,'std_dev_wvc_tb_in',ir,slab_size)
    std_dev_wvc_tb_out=get_sds(sd_id,'std_dev_wvc_tb_out',ir,slab_size)
    cell_lat= get_sds(sd_id,'cell_lat',ir,slab_size)
    cell_lon= get_sds(sd_id,'cell_lon',ir,slab_size)
    cell_azimuth= get_sds(sd_id,'cell_azimuth',ir,slab_size)
    cell_incidence= get_sds(sd_id,'cell_incidence',ir,slab_size)
    sigma0= get_sds(sd_id,'sigma0',ir,slab_size)
    sigma0_attn_amsr= get_sds(sd_id,'sigma0_attn_amsr',ir,slab_size)
    sigma0_attn_map= get_sds(sd_id,'sigma0_attn_map',ir,slab_size)
    kp_alpha= get_sds(sd_id,'kp_alpha',ir,slab_size)
    kp_beta= get_sds(sd_id,'kp_beta',ir,slab_size)
    kp_gamma= get_sds(sd_id,'kp_gamma',ir,slab_size)
    sigma0_qual_flag= get_sds(sd_id,'sigma0_qual_flag',ir,slab_size)
    sigma0_mode_flag= get_sds(sd_id,'sigma0_mode_flag',ir,slab_size)
    surface_flag= get_sds(sd_id,'surface_flag',ir,slab_size)
    cell_index= get_sds(sd_id,'cell_index',ir,slab_size)

; Print results to screen
    print,' '
    print,'TIME: ',string(wvc_row_time(*,ir))
    print,'WVC Row Number: ',row_number
    print,'Number of Sigma0: ',num_sigma0
    print,'Number of Sigma0 Per Cell:'
    print,num_sigma0_per_cell,format='(4(16i3/),1(12i3/))'
    print,'Number of Horizontal Brightness Temperatures'
    print,num_wvc_tb_in,format='(4(16i3/),1(12i3/))'
    print,'Mean of the Horizontal Brightness Temperatures'
    print,mean_wvc_tb_in,format='(4(16f7.2/),1(12f7.2/))'
    print,'Standard Deviation of the Horizontal Brightness Temperatures'
    print,std_dev_wvc_tb_in,format='(4(16f7.2/),1(12f7.2/))'
    print,'Number of Vertical Brightness Temperatures'
    print,num_wvc_tb_out,format='(4(16i3/),1(12i3/))'
    print,'Mean of the Vertical Brightness Temperatures'
    print,mean_wvc_tb_out,format='(4(16f7.2/),1(12f7.2/))'
    print,'Standard Deviation of the Vertical Brightness Temperatures'
    print,std_dev_wvc_tb_out,format='(4(16f7.2/),1(12f7.2/))'

    if (num_sigma0(0) gt 0) then begin
      print,' WVC   Lat   Lon    Azi  IncAng Sigma0  Atten  Kp_alpha   Kp_beta   Kp_gamma  Qual  Mode Surf'

      for j=0,num_sigma0(0)-1 do begin
        print,cell_index(j),cell_lat(j),cell_lon(j),cell_azimuth(j), $
        cell_incidence(j),sigma0(j),sigma0_attn_map(j),kp_alpha(j), $
        kp_beta(j),kp_gamma(j),sigma0_qual_flag(j),sigma0_mode_flag(j), $
        surface_flag(j), $
        format='(i3,5(1x,f6.2),1x,f6.3,1x,f9.6,2(1x,e10.4),x,i4,1x,i6,1x,i5)'
      endfor
    endif
  endfor

  HDF_SD_END,sd_id
end