; Copyright (C) 2002-2003  NASA/TRMM Satellite Validation Office
;
; This program is free software; you can redistribute it and/or modify
; it under the terms of the GNU General Public License as published by
; the Free Software Foundation; either version 2 of the License, or
; (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
;**************************************************************************
;

pro utdate, year, month, day

; Return current date in UTC.

datestr = systime(/utc)
month = strmid(datestr,4,3)
moname=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
i = 0
while moname[i] ne month do i = i + 1
month = i + 1
day = fix(strmid(datestr,8,2))
year = fix(strmid(datestr,20,4))
end

;***************************;
;      rsl_radar_to_uf      ;
;***************************;

pro rsl_radar_to_uf, radar, uf_file, fields=fields

;***********************************************************************
; Write the data from a Radar structure to a file in Universal Format.
;
; Syntax:
;     rsl_radar_to_uf, radar, uf_file [, FIELDS=string_array]
;
; Arguments:
;     radar:    a Radar data structure.
;     uf_file:  a string expression containing the name of the output UF file
;               to be created.
; Keywords:
;     FIELDS    string array (or scalar) containing the fields to be written
;               to the output file.  Default is all fields.  Fields are in the
;               form of the 2-character field names used by RSL, such as
;               'DZ', 'VR', etc.
;
; Written by:  Bart Kelley, GMU, October 2002
;
; Based on the RSL program RSL_radar_to_uf by John Merritt of SM&A Corp.
;
; Method:
; In UF, each ray (record) contains data for all fields used, and rays are
; written in the order recorded (all rays for sweep 1, followed by all rays for
; sweep 2, and so on), whereas in the radar structure the data is organized by
; fields (called volumes), then by sweeps and rays.  The following pseudo code
; describes how the transformation from radar structure to UF is done.
;
; for each sweep in radar structure do
;     for each ray do
;        for each volume do
; 	   if ray exists for this volume, put its data into record
;          put data header info into record.
;        done for volumes
;        swap bytes if little endian machine.
;        write record with fortran record length descriptors.
;     done for ray
; done for sweep
;***********************************************************************

rsl_speed_of_light = 299792458.0
no_data = fix(-32768)
do_opthdr = 1

openw, lunit, uf_file, /get_lun ; Open UF file for writing.

; Prepare character strings for header: Convert them to integer so they can
; be stored in the array used for UF record, and swap bytes if this is a little
; endian machine.  Strings don't need to be byte swapped, but we do it now
; because the entire record will be byte swapped later, and they will be in
; correct order.

uf = fix(byte('UF'),0,1)
byteorder,uf,/sswap, /swap_if_little_endian
radarname = byte(radar.h.radar_name)
if n_elements(radarname) lt 8 then $
    radarname = [radarname,bytarr(8-n_elements(radarname))]
if strlen(radar.h.name) eq 0 then $
    radarsite = radarname $
else begin
    radarsite = byte(radar.h.name)
    if n_elements(radarsite) lt 8 then $
        radarsite = [radarsite,bytarr(8-n_elements(radarsite))]
endelse
radarname=fix(radarname,0,4)
radarsite=fix(radarsite,0,4)
byteorder,radarname,/sswap, /swap_if_little_endian
byteorder,radarsite,/sswap, /swap_if_little_endian
tz = fix(byte('UT'),0,1)
byteorder,tz,/sswap, /swap_if_little_endian

nvolumes = radar.h.nvolumes
nsweeps = max(radar.volume.h.nsweeps)

buf = intarr(nvolumes * n_elements(radar.volume.sweep[0].ray[0].range))
dat = buf
nrec = 0

for iswp = 0, nsweeps-1 do begin
    maxrays = max(radar.volume.sweep[iswp].h.nrays)
    for iray = 0, maxrays-1 do begin
	; Search fields for a ray to use for header data, one with nbins > 0.
        ; If none found, skip the record building/writing process for this ray.
	havevol = where(radar.volume.sweep[iswp].ray[iray].h.nbins) 
	if size(havevol,/n_dimensions) gt 0 then begin
	    ray = radar.volume[havevol[0]].sweep[iswp].ray[iray]
	endif else continue
	 
	; Put mandatory header into record.
	replicate_inplace,buf,0
	replicate_inplace,dat,0
	buf[0] = uf
	buf[6] = 1 
	buf[8] = 1
	buf[9] = iswp + 1
	buf[10:13] = radarname
	buf[14:17] = radarsite
	buf[18] = radar.h.latd
	buf[19] = radar.h.latm
	buf[20] = round(radar.h.lats * 64.)
	buf[21] = radar.h.lond
	buf[22] = radar.h.lonm
	buf[23] = round(radar.h.lons * 64.)


	; If ray contains lat and lon, get them from ray.
	lat = ray.h.lat
	lon = ray.h.lon
	if lat ne 0.0 then begin
	    sign = long(abs(lat)/lat)
	    lat = abs(lat)
	    latd=floor(lat)
	    latm=floor((lat-latd) * 60.)
	    lats=floor(((lat-latd) * 60. - latm) * 60.)
	    lats=round(lats * 64.)
	    buf[18] = sign * latd
	    buf[19] = sign * latm
	    buf[20] = sign * lats
	endif
	if lon ne 0.0 then begin
	    sign = long(abs(lon)/lon)
	    lon = abs(lon)
	    lond=floor(lon)
	    lonm=floor((lon-lond) * 60.)
	    lons=floor(((lon-lond) * 60. - lonm) * 60.)
	    lons=round(lons * 64.)
	    buf[21] = sign * lond
	    buf[22] = sign * lonm
	    buf[23] = sign * lons
	endif
	buf[24] = radar.h.height
	buf[25] = ray.h.year mod 100
	buf[26] = ray.h.month
	buf[27] = ray.h.day
	buf[28] = ray.h.hour
	buf[29] = ray.h.minute
	buf[30] = ray.h.sec
	buf[31] = tz
	buf[32] = round(ray.h.azimuth * 64.)
	buf[33] = round(ray.h.elev * 64.)
	if radar.h.scan_mode eq 'PPI' then scan_mode = 1 $
	else if radar.h.scan_mode eq 'RHI' then scan_mode = 3 $
	else if radar.h.scan_mode eq 'Manual' then scan_mode = 6
	buf[34] = scan_mode
	buf[35] = round(radar.volume[havevol[0]].sweep[iswp].h.elev * 64.)
	buf[36] = round(ray.h.azim_rate * 64.)
	if buf[36] eq 0 then buf[36] = round(ray.h.sweep_rate * (360./60.)*64.)
	utdate, year, month, day
	buf[37] = year mod 100
	buf[38] = month
	buf[39] = day
	tmpbuf = fix(byte('RSIDL0.0'),0,4)
	byteorder, tmpbuf, /sswap, /swap_if_little_endian
	buf[40:43] = tmpbuf
	buf[44] = no_data
	mandihdr_len = 45

	; As in the C version, write optional header only once.
	opthdr_len = 0
	if do_opthdr then begin
	    do_opthdr = 0
	    tmpbuf = fix(byte('TRMMGVUF'),0,4)
	    byteorder, tmpbuf, /sswap, /swap_if_little_endian
	    buf[45:48] = tmpbuf
	    buf[49:50] = no_data
	    buf[51] = ray.h.hour
	    buf[52] = ray.h.minute
	    buf[53] = ray.h.sec
	    tmpbuf = fix(byte('RADAR_UF'),0,4)
	    byteorder, tmpbuf, /sswap, /swap_if_little_endian
	    buf[54:57] = tmpbuf
	    buf[58] = 2
	    opthdr_len = 14
	endif
	locuhdr_len = 0
	do_locuhdr = 0

	; If we have DZ and VR, check to see if their azimuths are
	; different. If they are, we will store VR azimuth in Local
	; Use Header. This happens with NEXRAD WSR88D radars, which
	; run separate sweeps for DZ and VR/SW at lower elevations.

	dzvol = where(radar.volume.h.field_type eq 'DZ', dzcount)
	vrvol = where(radar.volume.h.field_type eq 'VR', vrcount)
	if dzcount eq 1 and vrcount eq 1 then begin
	    vr_az = radar.volume[vrvol].sweep[iswp].ray[iray].h.azimuth
	    if vr_az ne radar.volume[dzvol].sweep[iswp].ray[iray].h.azimuth $
	        then do_locuhdr = 1
	endif
	if do_locuhdr then begin
	    lu_start = mandihdr_len + opthdr_len
	    azlabel = fix(byte('AZ'),0)
	    byteorder, azlabel, /sswap, /swap_if_little_endian
	    buf[lu_start] = azlabel
	    buf[lu_start+1] = round(vr_az * 64.)
	    locuhdr_len = 2
	endif
	
        ; Start of data header.
	dh_start = mandihdr_len + opthdr_len + locuhdr_len
	nfields = 0
	indx = 0 ; data buffer index
	indxdh = dh_start + 3 ; data header index

	; Check for selected fields list.
	nselect = n_elements(fields)
	if nselect gt 0 then sel_fields = strupcase(fields)

        ; Loop through volumes (fields) for this ray and get data.
	for ivol = 0, nvolumes-1 do begin
	    ; If ray.h.nbins not 0, then the ray includes this field.
	    if radar.volume[ivol].sweep[iswp].ray[iray].h.nbins ne 0 then begin
		fieldname = radar.volume[ivol].h.field_type
		; If fields have been selected, check if this is one of them.
		if nselect gt 0 then begin
		    n = where(sel_fields eq fieldname, count)
		    if count eq 0 then continue
		endif
	        nfields = nfields + 1
		tmpbuf = fix(byte(fieldname),0,1)
	        byteorder, tmpbuf, /sswap, /swap_if_little_endian
		buf[indxdh] = tmpbuf 
		buf[indxdh + 1] = indx
	        indxdh = indxdh + 2
		ray = radar.volume[ivol].sweep[iswp].ray[iray]
		scale_factor = 100
		if fieldname eq 'PH' then scale_factor = 10
	        ; Make field header.
		dat[indx+1] = scale_factor
		dat[indx+2] = ray.h.range_bin1/1000
		dat[indx+3] = ray.h.range_bin1 - dat[indx+2]*1000
		dat[indx+4] = ray.h.gate_size
		dat[indx+5] = ray.h.nbins
		dat[indx+6] = round(ray.h.pulse_width*rsl_speed_of_light/1.0e6)
		dat[indx+7] = round(ray.h.beam_width * 64.)
		dat[indx+8] = dat[7]
		dat[indx+9] = round(ray.h.frequency * 64.)
		dat[indx+10] = 0 ; horizontal polarization
		dat[indx+11] = round(ray.h.wavelength * 100. * 64.)
		dat[indx+12] = ray.h.pulse_count
		dat[indx+13] = fix(byte('  '),0,1)
		dat[indx+14:indx+15] = no_data
		if fieldname eq 'DZ' or fieldname eq 'ZT' then $
		    dat[indx+16] = round(radar.volume[ivol].h.calibr_const * 100.) $
		else dat[indx+16] = fix(byte('  '),0,1)
		if ray.h.prf ne 0 then dat[indx+17] = round(1./ray.h.prf*1.0e6) $
		else dat[indx+17] = no_data
		dat[indx+18] = 16
		datastart = indx+19
		if fieldname eq 'VR' or fieldname eq 'VE' then begin
		    dat[indx+19] = ray.h.nyq_vel * 100.
		    dat[indx+20] = 1
		    datastart = datastart + 2
		endif
		dat[indx] = datastart
		; Put data for this field into buffer.
		nbins = ray.h.nbins
		raydata = ray.range[0:nbins-1] * float(dat[indx+1])
		s = where(ray.range[0:nbins-1] eq $
		    radar.volume[ivol].h.no_data_flag)
		if size(s,/n_dimensions) gt 0 then raydata[s] = no_data
		indx = datastart
		dat[indx:indx+nbins-1] = round(raydata)
		indx = indx + nbins
	    endif
	endfor ; volumes

	; Compute offsets and build UF record.
	datalen = indx
	buf[[dh_start,dh_start+2]] = nfields
	buf[dh_start+1] = 1 ; number of records for ray.
	datahdr_len = 3 + 2*nfields
	data_start = dh_start + datahdr_len
	buf[data_start:data_start + datalen-1] = dat[0:datalen-1]
	for i = 0, nfields-1 do begin
	    indxdh = dh_start + 4 + i*2
	    buf[indxdh] = buf[indxdh] + data_start + 1
	    fldhdr_start = buf[indxdh] - 1
	    buf[fldhdr_start] = buf[fldhdr_start] + data_start + 1
	endfor

	recordlen = data_start + datalen
	buf[1] = recordlen
	buf[2] = mandihdr_len + 1
	buf[3] = mandihdr_len + opthdr_len + 1
	buf[4] = dh_start + 1
	nrec = nrec + 1
	buf[5] = nrec
	buf[7] = radar.volume[0].sweep[iswp].ray[iray].h.ray_num
	    
	; Swap bytes if this is little endian machine.
	byteorder, buf,/sswap, /swap_if_little_endian
	recordlen_bytes = recordlen * 2L
        byteorder, recordlen_bytes , /lswap, /swap_if_little_endian
	; Write record to UF.
	writeu, lunit, recordlen_bytes ; Write fortran record length descriptor.
	writeu, lunit, buf[0:recordlen-1]
	writeu, lunit, recordlen_bytes
    endfor ; rays
endfor ; sweeps
free_lun, lunit
		
end