; File: /home/stk/usr/code/src.lib/fit_reduce_new.pro
; Last Modification: 05-JUL-2010
; Author: Dieter Andre'
;
; modified to deal with a variable number of range gates
; up to 300
;
; modified to read both old and new fit files [uncompressed]
;
; modified to use the stk routines
;
; modified to work both for hp-ux and linux (swap_endian)
;
; Read a fit file and produce a fitred file
; that only contains important variables
;
; Since data for the years 1993 to 2001 were produced using hp-ux,
; it was decided to continue producing them in the same big-endian
; format under x86-linux, which is little-endian.
;
; on a 32 bit system arch = x86, an a 64 bit system arch = x86_64
; 20100521 DAndre
;
; added new radars [fhe, fhw, fir, mcm]; 220100705 DAndre


;@/home/stk/idl/lib/time.pro
;@/home/stk/idl/lib/datamap.pro
;@/home/stk/idl/lib/fit.pro
;@/home/stk/idl/lib/rprm.pro
;@/home/stk/idl/lib/oldfit.pro
;@/home/stk/idl/lib/genlib.pro


; Definition of the fit_red data structure

PRO fitdef_red, num_ranges, fred

  fred = { rec_time: long(0), $
	p: { parms_red, $
		st_id: 0, $
		year: 0, $
		month: 0, $
		day: 0, $
		hour: 0, $
		minut: 0, $
		sec: 0, $
		lagfr: 0, $
		smsep: 0, $
		nrang: 0, $
		frang: 0, $
		rsep: 0, $
		bmnum: 0, $
		xcf: 0, $
		tfreq: 0, $
		scan: 0, $
		cp: 0 $
		}, $
	qflag: INTARR( num_ranges), $
	pwr_l: FLTARR( num_ranges), $
	vel: FLTARR( num_ranges), $
	width_l: FLTARR( num_ranges), $
	gscat: INTARR( num_ranges), $
	x_qflag: INTARR( num_ranges), $
	phi0: FLTARR( num_ranges), $
	elev: FLTARR( num_ranges) $
	}
  RETURN
END ; fitdef_red


; Read a fitfile and produce a fitred file
; write fitred file to outdir, if given, otherwise .
; flag values as bad, where verr >= v

PRO fit_reduce, fitfile_name, OUTDIR= outdir, BAD=bad

  ; check if files are uncompressed
  IF ( STRPOS( fitfile_name, 'gz') NE (-1) ) THEN BEGIN
    print,'gzip compressed files cannot be read: ', fitfile_name
    GOTO, nochmal
  ENDIF
  IF ( STRPOS( fitfile_name, 'bz2') NE (-1) ) THEN BEGIN
    print,'bzip2 compressed files cannot be read: ', fitfile_name
    GOTO, nochmal
  ENDIF
  ; check if this is a fit file
  IF (( STRPOS( fitfile_name, 'fitacf') EQ (-1) ) AND ( STRPOS( fitfile_name, 'fit') EQ (-1) ) )THEN BEGIN
    print,'This is neither an old style nor a new style fit file: ', fitfile_name
    GOTO, nochmal
  ENDIF
  
  new_file= ( STRPOS( fitfile_name, 'fitacf') NE (-1) )
  
  IF new_file THEN BEGIN
    ifileptr= FitOpen( fitfile_name, /READ)
    IF (ifileptr LE 0) THEN BEGIN
      print,'Could not open the file', fitfile_name
      GOTO, nochmal
    ENDIF
  ENDIF ELSE BEGIN
    ifileptr= OldFitOpen( fitfile_name)
    IF (ifileptr.fitunit LE 0) THEN BEGIN
      print,'Could not open the file', fitfile_name
      GOTO, nochmal
    ENDIF
  ENDELSE
 
  ; make sure you start from the beginning of the file
  yr= 93
  mo= 1
  dy= 1
  hr= 0
  mn= 0
  sc= 0
  IF new_file THEN BEGIN
    status= FitSeek( ifileptr, yr, mo, dy, hr, mn, sc, atme=atme)
    status= FitRead( ifileptr, prm, fit)
  ENDIF ELSE BEGIN
    status= OldFitSeek( ifileptr, yr, mo, dy, hr, mn, sc, atme=atme)
    status= OldFitRead( ifileptr, prm, fit)
  ENDELSE

  num_rang= prm.nrang > 75
  nrm= num_rang - 1
  fitdef_red, num_rang, fred_str

  print, status, prm.bmnum, prm.time.hr, prm.time.mt, prm.time.sc

  IF KEYWORD_SET( outdir) THEN BEGIN
    out_name= outdir + strmid( fitfile_name, 0, rstrpos( fitfile_name,'.')+1) + 'fitred'
  ENDIF ELSE BEGIN
    out_name= strmid( fitfile_name, 0, rstrpos( fitfile_name,'.')+1) + 'fitred'
  ENDELSE
  OPENW, fr_un, out_name, /GET_LUN

  swap_it= (!version.arch EQ 'x86') OR (!version.arch EQ 'x86_64')	; swap endian
  WHILE (status EQ 0) DO BEGIN
    ; transfer values, omit bad records
    fred_str.rec_time= cnv_mdhms_sec( prm.time.yr, prm.time.mo, prm.time.dy, prm.time.hr, prm.time.mt, prm.time.sc)
    IF (fred_str.rec_time GT 0L) THEN BEGIN
      fred_str.p.st_id=  prm.stid
      fred_str.p.year=   prm.time.yr
      fred_str.p.month=  prm.time.mo
      fred_str.p.day=    prm.time.dy
      fred_str.p.hour=   prm.time.hr
      fred_str.p.minut=  prm.time.mt
      fred_str.p.sec=    prm.time.sc
      fred_str.p.lagfr=  prm.lagfr
      fred_str.p.smsep=  prm.smsep
      fred_str.p.nrang=  prm.nrang
      fred_str.p.frang=  prm.frang
      fred_str.p.rsep=   prm.rsep
      fred_str.p.bmnum=  prm.bmnum
      fred_str.p.tfreq=  prm.tfreq
      fred_str.p.scan=   prm.scan
      fred_str.p.cp=     prm.cp
      fred_str.qflag=    fit.qflg[ 0: nrm]
      ; to account for a bug in fitacf
      IF KEYWORD_SET( bad) THEN BEGIN
        good_fit= WHERE( (fit.qflg[ 0: nrm] EQ 1) AND (fit.gflg[ 0: nrm] EQ 0), igood)
	IF ( igood GT 0) THEN BEGIN
;          bad_fit= WHERE( (fit.v_e[ good_fit] GE ABS( fit.v[ good_fit])) OR (fit.p_l[ good_fit] LT 0.0), ibad)
;          bad_fit= WHERE( (fit.v_e[ good_fit] GE ABS( fit.v[ good_fit])) AND ( ABS( fit.v[ good_fit]) GT 3200.0), ibad)
          bad_fit= WHERE( (fit.v_e[ good_fit] GE ABS( fit.v[ good_fit])), ibad)
          IF (ibad NE 0) THEN BEGIN
	    fred_str.qflag[ good_fit[ bad_fit]]= 0
	  ENDIF
	ENDIF
      ENDIF
      fred_str.pwr_l=    fit.p_l[ 0: nrm]
      fred_str.vel=      fit.v[ 0: nrm]
      fred_str.width_l=  fit.w_l[ 0: nrm]
      fred_str.gscat=    fit.gflg[ 0: nrm]
      fred_str.x_qflag=  fit.x_qflg[ 0: nrm]
      fred_str.phi0=     fit.phi0[ 0: nrm]
      fred_str.elev=     fit.elv[ 0: nrm]
      IF swap_it THEN fred_str= swap_endian( fred_str)
      WRITEU, fr_un, fred_str
    ENDIF    
    IF new_file THEN BEGIN
      status= FitRead( ifileptr, prm, fit)
    ENDIF ELSE BEGIN
      status= OldFitRead( ifileptr, prm, fit)
    ENDELSE
    ;print, status, fit_str.p.bmnum, fit_str.p.hour, fit_str.p.minut, fit_str.p.sec
  ENDWHILE
  CLOSE, fr_un
  FREE_LUN, fr_un
  nochmal:
  IF new_file THEN BEGIN
    CLOSE, ifileptr
    FREE_LUN, ifileptr
  ENDIF ELSE BEGIN
    status= OldFitClose( ifileptr)
  ENDELSE
END ; fit_reduce


; Read data from a fitred file
; The data structure is an array 16x870 of fit_red beams
; Each fit_red is put into the array at its beam position
; and when the last beam is reached, the record number is increased
; Therefore read_fred only works for normal_scan and other
; well behaved radar modes
; The number of records was limited to 870 because of hp-ux memory limitations.
; For fast_scan the program will stop reading shortly before 15h
; This program also assumes, that the fitred files are in 1 day chunks.

PRO read_fred,	fredfile_name,		$ ; read from this fitred file
		fred_rec_ptr,		$ ; pointer to the data read
		START_HOUR= start_hour,	$ ; if set, read starting from this hour of day, otherwise 0
		STOP_HOUR= stop_hour,	$ ; if set, stop reading at this hour of day, otherwise 24
		NS=ns,			$ ; if set, read only normal_scan and fast_scan data
		GZ=gz			  ; if set, read from a gzipped fitred file [works only from IDL 5.3 and up]

  swap_it= (!version.arch EQ 'x86') OR (!version.arch EQ 'x86_64')	; swap endian for x86
  IF swap_it THEN BEGIN
    num_recs= 1441 ; one day of fast_scan for x86
  ENDIF ELSE BEGIN
    num_recs= 870 ; enough for a day with 100 s integration time; for hp_ux due to memory limit
  ENDELSE

  ; read the first record to find the actual number of range gates
  fitdef_red, 75, fred_str
  IF KEYWORD_SET( gz) THEN BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN, /COMPRESS
  ENDIF ELSE BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN
  ENDELSE
  READU, fr_un, fred_str
  IF swap_it THEN fred_str= swap_endian( fred_str)
  ; all previous records [before rkn have been created with 75 range gates]
  num_rang= fred_str.p.nrang > 75
  CLOSE, fr_un
  FREE_LUN, fr_un
  fred_str= 0
  
  fitdef_red, num_rang, fred_str
  fred_arr= REPLICATE( fred_str, 16, num_recs)
  IF KEYWORD_SET( gz) THEN BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN, /COMPRESS
  ENDIF ELSE BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN
  ENDELSE
  READU, fr_un, fred_str
  IF swap_it THEN fred_str= swap_endian( fred_str)
  
  ; this does not work, if the mode changes during the day
  IF KEYWORD_SET( ns) AND ( MAX( ABS( fred_str.p.cp)) GT 151) THEN BEGIN ; use normal_scan and fast_scan only
    irec= 0
    GOTO, finish    
  ENDIF
  IF NOT KEYWORD_SET( stop_hour) THEN stop_hour= 24
  IF KEYWORD_SET( start_hour) THEN BEGIN
    WHILE ( NOT eof( fr_un) ) AND ( fred_str.p.hour LT start_hour) DO BEGIN
      READU, fr_un, fred_str
      IF swap_it THEN fred_str= swap_endian( fred_str)
    ENDWHILE
  ENDIF
  ; last radar taken into account: Tiger r 14
  IF	(fred_str.p.st_id EQ 7)		$ ; Kodiak
	OR (fred_str.p.st_id EQ 5)	$ ; Saskatoon
	OR (fred_str.p.st_id EQ 1)	$ ; Goose Bay
	OR (fred_str.p.st_id EQ 9)	$ ; Iceland East
	OR (fred_str.p.st_id EQ 4)	$ ; Halley Bay
	OR (fred_str.p.st_id EQ 13)	$ ; Syowa East
	OR (fred_str.p.st_id EQ 18)	$ ; Unwin
	OR (fred_str.p.st_id EQ 20)	$ ; McMurdo
	OR (fred_str.p.st_id EQ 21)	$ ; Falkland
	OR (fred_str.p.st_id EQ 64)	$ ; Inuvik
	OR (fred_str.p.st_id EQ 205)	$ ; Fort Hays East
	THEN $
    last_beam= 15 $ ; forward scan
  ELSE $
    last_beam= 0 ; backward scan
  irec= 0
  WHILE ( NOT eof( fr_un) ) AND (irec LT num_recs) AND (fred_str.p.hour LT stop_hour) DO BEGIN
    ; this DOES work, if the mode changes during the day
    IF KEYWORD_SET( ns) AND ( MAX( ABS( fred_str.p.cp)) GT 151) THEN BEGIN ; use normal_scan and fast_scan only
      irec= 0
      GOTO, finish    
    ENDIF
    ; !!! THIS should not be necessary, but some Kodiak files have beam number -1 !!!
    fred_arr[ fred_str.p.bmnum > 0, irec]= fred_str
    IF (fred_str.p.bmnum EQ last_beam) THEN $
      irec= irec + 1
    READU, fr_un, fred_str
    IF swap_it THEN fred_str= swap_endian( fred_str)
  ENDWHILE
  finish:
  CLOSE, fr_un
  FREE_LUN, fr_un
  IF (irec GT 0) THEN BEGIN
    fred_arr= TEMPORARY( fred_arr[ *, 0: irec - 1])
    fred_rec_ptr= PTR_NEW( fred_arr, /NO_COPY)
  ENDIF ELSE BEGIN
    fred_rec_ptr= PTR_NEW()
  ENDELSE
  fred_str= 0 
  fred_arr= 0 
END; read_fred


; select data for one beam from a fitred file

PRO select_fred, fredfile_name,		$ ; fitred filename
		fred_rec_ptr,		$ ; pointer to the data read
		wanted_beam,		$ ; get data from this beam
		START_HOUR= start_hour,	$ ; if set, read starting from this hour of day, otherwise 0
		STOP_HOUR= stop_hour,	$ ; if set, stop reading at this hour of day, otherwise 24
		GZ= gz			  ; if set, read from a gzipped fitred file [works only from IDL 5.3 and up]

  swap_it= (!version.arch EQ 'x86') OR (!version.arch EQ 'x86_64')	; swap endian
  ; num_recs= 2000 ; enough for a beam a day
  num_recs= 25000 ; increased to handle a 24hour Themis run
  
  ; read the first record to find the actual number of range gates
  fitdef_red, 75, fred_str
  IF KEYWORD_SET( gz) THEN BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN, /COMPRESS
  ENDIF ELSE BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN
  ENDELSE
  READU, fr_un, fred_str
  IF swap_it THEN fred_str= swap_endian( fred_str)
  ; all previous records [before rkn have been created with 75 range gates]
  num_rang= fred_str.p.nrang > 75
  CLOSE, fr_un
  FREE_LUN, fr_un
  fred_str= 0
  
  fitdef_red, num_rang, fred_str
  fred_arr= REPLICATE( fred_str, num_recs)
  IF KEYWORD_SET( gz) THEN BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN, /COMPRESS
  ENDIF ELSE BEGIN
    OPENR, fr_un, fredfile_name, /GET_LUN
  ENDELSE
  READU, fr_un, fred_str
  IF swap_it THEN fred_str= swap_endian( fred_str)
  IF NOT KEYWORD_SET( stop_hour) THEN stop_hour= 24
  IF KEYWORD_SET( start_hour) THEN BEGIN
    WHILE ( NOT eof( fr_un) ) AND ( fred_str.p.hour LT start_hour) DO BEGIN
      READU, fr_un, fred_str
      IF swap_it THEN fred_str= swap_endian( fred_str)
    ENDWHILE
  ENDIF
  irec= 0
  WHILE ( NOT eof( fr_un) ) AND (irec LT num_recs)  AND (fred_str.p.hour LT stop_hour) DO BEGIN
    IF (fred_str.p.bmnum EQ wanted_beam) THEN BEGIN
      fred_arr[ irec]= fred_str
      irec= irec + 1
    ENDIF
    READU, fr_un, fred_str
    IF swap_it THEN fred_str= swap_endian( fred_str)
  ENDWHILE
  IF (irec GT 0) THEN BEGIN
    fred_arr= TEMPORARY( fred_arr[ 0: irec - 1])
  ENDIF ELSE BEGIN
    fred_arr= TEMPORARY( fred_arr[ 0])
  ENDELSE
  CLOSE, fr_un
  FREE_LUN, fr_un
  fred_rec_ptr= PTR_NEW( fred_arr, /NO_COPY)
  fred_str= 0
  fred_arr= 0 
END; select_fred


PRO fred_ascii, shr, ehr, FITRED_FILE= fitred_file, BEAM_LIST=beam_list

  ; ask for a fitred file to open, if not given as parameter
  IF NOT KEYWORD_SET( fitred_file) THEN BEGIN
    fitred_directory= ''
    fitred_filename= pickfile( /READ, PATH= '/jukebox/fitred', FILTER= '*.fitred*', GET_PATH=fitred_directory)
    CD, fitred_directory
    last_slash= RSTRPOS( fitred_filename, '/')
    fitred_filename= STRMID( fitred_filename, last_slash+1, STRLEN( fitred_filename))
  ENDIF ELSE BEGIN
    fitred_filename= fitred_file
  ENDELSE
    
  fred_rec_ptr= PTR_NEW()
  IF ( RSTRPOS( fitred_filename, '.gz') NE -1) THEN BEGIN
    read_fred, fitred_filename, fred_rec_ptr, start_hour= shr, stop_hour= ehr, /GZ
  ENDIF ELSE BEGIN
    read_fred, fitred_filename, fred_rec_ptr, start_hour= shr, stop_hour= ehr
  ENDELSE
  
  IF NOT PTR_VALID( fred_rec_ptr) THEN BEGIN
    PRINT, 'No data in ' + fitred_filename
    STOP
  ENDIF

  IF NOT KEYWORD_SET( beam_list) THEN beam_list= INDGEN( 16)
  nbm= N_ELEMENTS( beam_list)

  p0= (*fred_rec_ptr)[ 0, 0].p
  fas_file= STRCOMPRESS( STRING( '~/fas_', p0.year, '_', p0.month, '_', p0.day, '_', shr, '_', ehr, '_', p0.st_id, '.txt'), /REMOVE_ALL)
  OPENW, fas_un, fas_file, /GET_LUN
  PRINTF, fas_un, 'File structure info:'
  PRINTF, fas_un, 'Beam header:'
  PRINTF, fas_un, 'Year Month Day Hour Minute Second Tx-Freq[MHz] First_Range[km] Range_Separation[km]'
  PRINTF, fas_un, 'Beam_Number Number_of_good_Range_Gates Station_Id'
  PRINTF, fas_un, 'For each range gate:'
  PRINTF, fas_un, 'Range_Gate_Number Power_l[dB] Velocity[m/s] Spectral_Width[m/s] Elevation-Angle[deg]'
  PRINTF, fas_un, 'Ground_Scatter_Flag[T means ground scatter]'
  PRINTF, fas_un, 'Invalid elevation angles are set to 99'
  ntm= SIZE( (*fred_rec_ptr))
  ntm= ntm[ 2]
  gsc_let= [ ' F', ' T']
  FOR itm= 0, ntm-1 DO BEGIN
    FOR jbi= 0, nbm-1 DO BEGIN
      ibi= beam_list[ jbi]
      one_beam= (*fred_rec_ptr)[ ibi, itm]
      good_pts= WHERE(	( one_beam.qflag EQ 1) AND ( one_beam.pwr_l GT 3.0), ngp)
      IF (ngp GT 0) THEN BEGIN
	PRINTF, fas_un, STRCOMPRESS( STRING( one_beam.p.year, one_beam.p.month, one_beam.p.day, $
		one_beam.p.hour, one_beam.p.minut, one_beam.p.sec, one_beam.p.tfreq/ 1000.0, $
		one_beam.p.frang, one_beam.p.rsep, ibi, ngp, one_beam.p.st_id ))
	FOR igp= 0, ngp - 1 DO BEGIN
	  rgt= good_pts[ igp]
	  t_elev= one_beam.elev[ rgt]
	  IF (t_elev EQ 0.0) THEN t_elev= 99
	  PRINTF, fas_un, STRCOMPRESS( STRING( rgt, FIX( one_beam.pwr_l[ rgt]), FIX( one_beam.vel[ rgt]), $
		FIX( one_beam.width_l[ rgt]), FIX( t_elev), gsc_let[ one_beam.gscat[ rgt]] ))
	ENDFOR
      ENDIF
    ENDFOR
  ENDFOR

  CLOSE, fas_un
  FREE_LUN, fas_un
END ; fred_ascii