;+
; PROCEDURE:
;       kgy_lrs_load
; PURPOSE:
;       loads Kaguya LRS natural wave data
; CALLING SEQUENCE:
;       timespan,'2008-01-01',2
;       kgy_lrs_load
; KEYWORDS:
;       types: 'NPW', 'WFC', or ['NPW','WFC'] (Def: ['NPW','WFC'])
;       version: data version (Def: '010')
;       append: if set, append to pre-loaded tplot variables (Def: clear old data)
;       wfcdatadir: local directory in which WFC files are stored
;       files: local files to read
;              if set, does not download files from the data archive site
; CREATED BY:
;       Yuki Harada on 2016-09-02
;
; $LastChangedBy: haraday $
; $LastChangedDate: 2021-07-12 00:51:24 -0700 (Mon, 12 Jul 2021) $
; $LastChangedRevision: 30121 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_5_0/projects/kaguya/lrs/kgy_lrs_load.pro $
;-

pro kgy_lrs_load, files=files, version=version, trange=trange, append=append, types=types, wfcdatadir=wfcdatadir, _extra=_extra


if ~keyword_set(version) then version = '010' ;- incapable of automatic version search
version2 = strmid(version,1,1)+'.'+strmid(version,2,1)
if ~keyword_set(append) then store_data,'kgy_lrs_*',/delete
if ~keyword_set(types) then types = ['NPW','WFC'] else types = strupcase(types)

;;; wfcdatadir = root_data_dir()+'kaguya/lrs/WFC-H/CDF2/'
tr = timerange(trange)

;;; retrieve files
if ~keyword_set(files) then begin
   files = ''

   ;;; NPW
   if total(strmatch(types,'NPW')) then begin
      pf = 'sln-l-lrs-5-npw-spectrum-v'+version2+'/YYYYMMDD/data/LRS_NPW_V'+version+'_YYYYMMDD.cdf'
;      pf = 'LRS_NPW_V'+version+'_YYYYMMDD' ;- obsolete
      f = kgy_file_retrieve(pf,trange=trange,/public,/valid, _extra=_extra)
      if total(strlen(f)) gt 0 then files = [files,f]
   endif

   ;;; WFC, skip odd hours
   if total(strmatch(types,'WFC')) then begin
      if keyword_set(wfcdatadir) then begin ;- find local files
         pf = wfcdatadir+time_intervals(tf='LRS_WFC_V'+version+'_YYYYMMDDhhmmss.cdf',/hourly,trange=tr)
         f = file_search(pf)
      endif else begin
         pf = 'sln-l-lrs-4-wfc-spectrum-v'+version2+'/YYYYMMDD/data/LRS_WFC_V'+version+'_YYYYMMDDhhmmss.cdf'
;      pf = 'LRS_WFC_V'+version+'_YYYYMMDDhhmmss' ;- obsolete
         f = kgy_file_retrieve(pf,trange=trange,/public,/hourly,/valid,/skipodd, _extra=_extra)
      endelse
      if total(strlen(f)) gt 0 then files = [files,f]
   endif
endif


;;; select files to read
w = where(file_test(files),nw)
if nw eq 0 then begin
   dprint,'No valid files'
   return
endif
files = files[w]


;;; read files
for ifile=0,n_elements(files)-1 do begin
   fname = files[ifile]
   cdfi = cdf_load_vars(fname,varformat='*',/convert_int1_to_int2)
   id = cdf_open(fname)
   dprint,'reading '+fname

   ;;; rename old tplot variable
   tn = tnames('kgy_lrs_*',ntn)
   for itn=0,ntn-1 do tplot_rename,tn[itn],tn[itn]+'_oldload'

   ;;; NPW
   if total(strmatch(fname,'*NPW*')) then begin
      ;;; varnames: Epoch Mode Frequency RX1 RX2

      ;;; get times
      idx = where( cdfi.vars.name eq 'Epoch' , nw )
      if nw ne 1 then continue
      v = cdfi.vars[idx]
      attr = *v.attrptr
      if v.numrec eq 0 then continue ;- no data
      data = *v.dataptr
      cdf_epoch,reform(data),yr,mo,dy,hr,mn,sc,ml,/BREAK
      times = time_double( string(yr,format='(i4.4)')+'-' $
                           +string(mo,format='(i2.2)')+'-' $
                           +string(dy,format='(i2.2)')+'/' $
                           +string(hr,format='(i2.2)')+':' $
                           +string(mn,format='(i2.2)')+':' $
                           +string(sc,format='(i2.2)')+'.' $
                           +string(ml,format='(i3.3)')+'.'  )

      ;;; get freq
      idx = where( cdfi.vars.name eq 'Frequency' , nw )
      if nw ne 1 then continue
      v = cdfi.vars[idx]
      data = *v.dataptr
      attr = *v.attrptr
      freq = data
      frequnits = attr.units

      ;;; get RX1
      idx = where( cdfi.vars.name eq 'RX1' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            data = float( *v.dataptr )
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then data[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_npw_rx1', $
                       data={x:times,y:data,v:freq}, $
                       dlim={ytitle:'LRS NPW!cE-field RX1!cFrequency!c['+frequnits+']',$
                             ylog:1,yticklen:-.01, $
                             yrange:minmax(freq),ystyle:1, $
                             spec:1,zlog:0,ztitle:attr.units, $
                             zrange:[-200,-120]}
         endif
      endif

      ;;; get RX2
      idx = where( cdfi.vars.name eq 'RX2' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            data = float( *v.dataptr )
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then data[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_npw_rx2', $
                       data={x:times,y:data,v:freq}, $
                       dlim={ytitle:'LRS NPW!cE-field RX2!cFrequency!c['+frequnits+']',$
                             ylog:1,yticklen:-.01, $
                             yrange:minmax(freq),ystyle:1, $
                             spec:1,zlog:0,ztitle:attr.units, $
                             zrange:[-200,-120]}
         endif
      endif

      ;;; get modes
      idx = where( cdfi.vars.name eq 'Mode' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            data = *v.dataptr
            store_data,'kgy_lrs_npw_mode',data={x:times,y:data}, $
                       dlim={ytitle:'LRS NPW!cMode',psym:1}
         endif
      endif

   endif                        ;- NPW


   ;;; WFC
   if total(strmatch(fname,'*WFC*')) then begin
      ;;; varnames: Ex Ey Frequency Epoch PDC-TI Mode Gain PostGap
      ;;; cf. muro-shuuron.pdf

      ;;; get times
      idx = where( cdfi.vars.name eq 'Epoch' , nw )
      if nw ne 1 then continue
      v = cdfi.vars[idx]
      attr = *v.attrptr
      if v.numrec eq 0 then continue ;- no data
      data = *v.dataptr
      cdf_epoch,reform(data),yr,mo,dy,hr,mn,sc,ml,/BREAK
      times = time_double( string(yr,format='(i4.4)')+'-' $
                           +string(mo,format='(i2.2)')+'-' $
                           +string(dy,format='(i2.2)')+'/' $
                           +string(hr,format='(i2.2)')+':' $
                           +string(mn,format='(i2.2)')+':' $
                           +string(sc,format='(i2.2)')+'.' $
                           +string(ml,format='(i3.3)')+'.'  )

      ;;; get freq
      idx = where( cdfi.vars.name eq 'Frequency' , nw )
      if nw ne 1 then continue
      v = cdfi.vars[idx]
      data = *v.dataptr
      attr = *v.attrptr
      freq = data
      frequnits = attr.units

      ;;; get Gain
      idx = where( cdfi.vars.name eq 'Gain' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            if v.numrec eq n_elements(times) then data = *v.dataptr $
            else begin
               cdf_varget,id,'Gain',data,rec_count=n_elements(times)
            endelse
            ;;; 3-4bit( 0x00:40dB 0x01:20dB 0x02:20dB 0x03:0dB )
            ;;; a = [0,0,0,0,1,1,0,0] & print,a.frombits()
            gain = ishft( data and 12b , -2 ) ;- 00001100 = 12b
            w = where( gain eq 0 , nw )
            if nw gt 0 then gain[w] = 40
            w = where( gain eq 1 or gain eq 2 , nw )
            if nw gt 0 then gain[w] = 20
            w = where( gain eq 3 , nw )
            if nw gt 0 then gain[w] = 0
            gain = reform(float(gain))
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then gain[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_wfc_gain', $
                       data={x:times,y:gain}, $
                       dlim={ytitle:'LRS WFC!cgain',psym:1,yrange:[0,40]}
         endif
      endif

      ;;; get Ex
      idx = where( cdfi.vars.name eq 'Ex' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            if v.numrec eq n_elements(times) then data = float( *v.dataptr ) $
            else begin
               cdf_varget,id,'Ex',data,rec_count=n_elements(times)
               data = float(transpose(data))
            endelse
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then data[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_wfc_Ex_dB', $
                       data={x:times,y:data,v:freq}, $
                       dlim={ytitle:'Ex!cFrequency!c['+frequnits+']', $
                             ylog:1,yticklen:-.01, $
                             yrange:minmax(freq)>1,ystyle:1, $
                             spec:1,zlog:0,ztitle:attr.units, $
                             zrange:[-20,60]}
            if n_elements(gain) eq n_elements(times) then begin
               ;;; gain correction
               ex = data - rebin(gain,n_elements(times),n_elements(freq))
               ex2 = ex - 15.5  ;- correct for sensor to ADC gain
               ex2[*,0:127] += 20. ;- WFCH-LOW
               ex2[*,128:*] += 10. ;- WFCH-HIGH
               store_data,'kgy_lrs_wfc_Ex_phys', $
                          data={x:times,y:ex2,v:freq}, $ ;- see WFCH1_calData() and WFCH1_getdbgain in wfch1_v2.c
                          dlim={ytitle:'Ex!cFrequency!c['+frequnits+']', $
                                ylog:1,yticklen:-.01, $
                                yrange:minmax(freq)>1,ystyle:1, $
                                spec:1,zlog:0,ztitle:'dB'+'!4l!X'+'V/m', $
                                zrange:[-60,40]}

               dfreq = freq * 0.
               ;;; see muro-shuuron.pdf
               ;;; 0-127: 0-9.77 kHz, WFC-H LOW
               dfreq[0:127] = (freq[1] - freq[0])
               ;;; 128-350: 9.77 kHz-1 MHz, WFC-H HIGH
               w = where( freq gt 9.77 and freq le 78.13 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               w = where( freq gt 78.13 and freq le 117.19 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               w = where( freq gt 117.19 and freq le 156.25 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               w = where( freq gt 156.27 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               dfreq *= 1e3
               pow = 10.^(ex2/10.-12.) $ ;- (V/m)^2/Hz
                     / transpose(rebin(dfreq,n_elements(freq),n_elements(times)))
               store_data,'kgy_lrs_wfc_Ex_phys2',data={x:times,y:pow,v:freq}, $
                          dlim={ytitle:'Ex!cFrequency!c['+frequnits+']', $
                                ylog:1,yticklen:-.01, $
                                yrange:minmax(freq)>1,ystyle:1, $
                                spec:1,zlog:1,ztitle:'(V/m)!u2!n/Hz', $
                                zrange:[1e-20,1e-10],dfreq:dfreq}
            endif
         endif
      endif

      ;;; get Ey
      idx = where( cdfi.vars.name eq 'Ey' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            if v.numrec eq n_elements(times) then data = float( *v.dataptr ) $
            else begin
               cdf_varget,id,'Ey',data,rec_count=n_elements(times)
               data = float(transpose(data))
            endelse
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then data[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_wfc_Ey_dB', $
                       data={x:times,y:data,v:freq}, $ ;- raw dB w/o gain correction
                       dlim={ytitle:'Ey!cFrequency!c['+frequnits+']', $
                             ylog:1,yticklen:-.01, $
                             yrange:minmax(freq)>1,ystyle:1, $
                             spec:1,zlog:0,ztitle:attr.units, $
                             zrange:[-20,60]}
            if n_elements(gain) eq n_elements(times) then begin
               ;;; gain correction
               ey = data - rebin(gain,n_elements(times),n_elements(freq))
               ey2 = ey - 15.5  ;- correct for sensor to ADC gain
               ey2[*,0:127] += 20. ;- WFCH-LOW
               ey2[*,128:*] += 10. ;- WFCH-HIGH
               store_data,'kgy_lrs_wfc_Ey_phys', $
                          data={x:times,y:ey2,v:freq}, $ ;- see WFCH1_calData() and WFCH1_getdbgain in wfch1_v2.c
                          dlim={ytitle:'Ey!cFrequency!c['+frequnits+']', $
                                ylog:1,yticklen:-.01, $
                                yrange:minmax(freq)>1,ystyle:1, $
                                spec:1,zlog:0,ztitle:'dB'+'!4l!X'+'V/m', $
                                zrange:[-60,40]}

               dfreq = freq * 0.
               ;;; see muro-shuuron.pdf
               ;;; 0-127: 0-9.77 kHz, WFC-H LOW
               dfreq[0:127] = (freq[1] - freq[0])
               ;;; 128-350: 9.77 kHz-1 MHz, WFC-H HIGH
               w = where( freq gt 9.77 and freq le 78.13 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               w = where( freq gt 78.13 and freq le 117.19 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               w = where( freq gt 117.19 and freq le 156.25 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               w = where( freq gt 156.27 , nw )
               dfreq[w] = median( freq[w] - shift(freq[w],1) )
               dfreq *= 1e3
               pow = 10.^(ey2/10.-12.) $ ;- (V/m)^2/Hz
                     / transpose(rebin(dfreq,n_elements(freq),n_elements(times)))
               store_data,'kgy_lrs_wfc_Ey_phys2',data={x:times,y:pow,v:freq}, $
                          dlim={ytitle:'Ey!cFrequency!c['+frequnits+']', $
                                ylog:1,yticklen:-.01, $
                                yrange:minmax(freq)>1,ystyle:1, $
                                spec:1,zlog:1,ztitle:'(V/m)!u2!n/Hz', $
                                zrange:[1e-20,1e-10],dfreq:dfreq}
               ;;; obsolete:
               ;; ey2 = ey-15.5
               ;; store_data,'kgy_lrs_wfc_Ey_phys', $
               ;;            data={x:times,y:ey2,v:freq}, $ ;- see WFCH1_calData_v1() and WFCH1_getdbgain in wfch1_v2.c
               ;;            dlim={ytitle:'Ey!cFrequency!c['+frequnits+']', $
               ;;                  ylog:1,yticklen:-.01, $
               ;;                  yrange:minmax(freq)>1,ystyle:1, $
               ;;                  spec:1,zlog:0,ztitle:'dB'+'!4l!X'+'V/m', $
               ;;                  zrange:[-80,30]}
            endif
         endif
      endif

      ;;; get Mode
      idx = where( cdfi.vars.name eq 'Mode' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            if v.numrec eq n_elements(times) then data = *v.dataptr $
            else begin
               cdf_varget,id,'Mode',data,rec_count=n_elements(times)
            endelse
; Amount of information is 1 byte.: 1-2bit: Channel( 0x01:X 0x02:Y
; 0x03:X-Y )  3-4bit: Frequency Band( 0x01:~10kHz 0x02:~1MHz
; 0x03:~10kHz ~1MHz )  5-6bit: Observation Mode( 0x00:WAVE 0x01:FFT
; 0x02:PHASE )
            xymode = reform(float( data and 3b )) ;- 00000011
            fband = reform(float( ishft( data and 12b , -2 ) )) ;- 00001100
            omode = reform(float( ishft( data and 48b , -4 ) )) ;- 00110000
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then begin
                  xymode[wnan] = !values.f_nan
                  fband[wnan] = !values.f_nan
                  omode[wnan] = !values.f_nan
               endif
            endif
            store_data,'kgy_lrs_wfc_xymode', $
                       data={x:times,y:xymode}, $
                       dlim={ytitle:'LRS WFC!cXY Ch',psym:1}
            store_data,'kgy_lrs_wfc_fband', $
                       data={x:times,y:fband}, $
                       dlim={ytitle:'LRS WFC!cFreq. Band',psym:1}
            store_data,'kgy_lrs_wfc_omode', $
                       data={x:times,y:omode}, $
                       dlim={ytitle:'LRS WFC!cMode',psym:1}
         endif
      endif

      ;;; get PDC-TI
      idx = where( cdfi.vars.name eq 'PDC-TI' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            if v.numrec eq n_elements(times) then data = float( *v.dataptr ) $
            else begin
               cdf_varget,id,'PDC-TI',data,rec_count=n_elements(times)
               data = float(transpose(data))
            endelse
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then data[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_wfc_pdc-ti', $
                       data={x:times,y:reform(data)}, $
                       dlim={ytitle:'LRS WFC!cPDC-TI',psym:1}
         endif
      endif

      ;;; get PostGap
      idx = where( cdfi.vars.name eq 'PostGap' , nw )
      if nw eq 1 then begin
         v = cdfi.vars[idx]
         attr = *v.attrptr
         if v.numrec gt 0 then begin
            if v.numrec eq n_elements(times) then data = float( *v.dataptr ) $
            else begin
               cdf_varget,id,'PostGap',data,rec_count=n_elements(times)
               data = float(reform(data))
            endelse
            if tag_exist(attr,'fillval') then begin
               wnan = where( data eq attr.fillval , nwnan )
               if nwnan gt 0 then data[wnan] = !values.f_nan
            endif
            store_data,'kgy_lrs_wfc_postgap', $
                       data={x:times,y:reform(data)}, $
                       dlim={ytitle:'LRS WFC!cPost Gap',psym:1}
         endif
      endif

   endif                        ;- WFC


   ;;; concat tplot variable
   tn = tnames('kgy_lrs_*_oldload',ntn)
   for itn=0,ntn-1 do begin
      get_data,tn[itn],data=dold,dlim=dlold
      newtn = strmid(tn[itn],0,strlen(tn[itn])-8)
      get_data,newtn,data=dnew,dlim=dlnew,dtype=dtype
      if dtype eq 0 then tplot_rename,tn[itn],newtn else begin
         newx = [ dold.x, dnew.x ]
         newy = [ dold.y, dnew.y ]
         if tag_exist(dold,'v') and tag_exist(dnew,'v') then begin
            if size(dold.v,/n_dim) ne size(dnew.v,/n_dim) then begin
               dprint,'v tag dims do not match: ',tn[itn],' ',newtn
               continue
            endif
            if size(dold.v,/n_dim) eq 1 then begin
               if total(abs(dold.v-dnew.v)) eq 0 then newv = dold.v else begin
                  newv = [ transpose(rebin(dold.v,n_elements(dold.v),n_elements(dold.x))), transpose(rebin(dnew.v,n_elements(dnew.v),n_elements(dnew.x))) ]
               endelse
            endif else newv = [ dold.v, dnew.v ]
            store_data,newtn,data={x:newx,y:newy,v:newv},dlim=dlnew
         endif else begin ;- no v tag
            store_data,newtn,data={x:newx,y:newy},dlim=dlnew
         endelse
      endelse
   endfor
   tn = tnames('kgy_lrs_*_oldload',ntn)
   if ntn gt 0 then store_data,tn,/delete


endfor                          ;- ifile


end