;+
; Read density from uh line for phasef data production.
; Adopted from rbsp_efw_get_flag_values.
;-

pro rbsp_efw_phasef_read_density_gen_file, time, probe=probe, filename=file


;---Check inputs.
    if n_elements(file) eq 0 then begin
        errmsg = handle_error('No output file ...')
        return
    endif

    if n_elements(probe) eq 0 then begin
        errmsg = handle_error('No input probe ...')
        return
    endif

    if n_elements(time) eq 0 then begin
        errmsg = handle_error('No input time ...')
        return
    endif

    if n_elements(boom_pair) eq 0 then boom_pair = '12'

    ; No density limit. Apply limit when reading density.
    dmin = 0
    dmax = 1e31


;---Constants and settings.
    secofday = 86400d
    errmsg = ''

    prefix = 'rbsp'+probe+'_'
    if size(time[0],/type) eq 7 then time = time_double(time)
    date = time[0]-(time[0] mod secofday)
    time_range = date+[0,secofday]
    timespan, time_range[0], total(time_range*[-1,1]), /second
    boom_pairs = ['12','34','13','14','23','24']

    time_step = 1.
    common_times = make_bins(time_range+[0,-1]*time_step, time_step)
    ntime = n_elements(common_times)

    ; Load waveform data
    vsvy_var = prefix+'efw_vsvy_interp'
    rbsp_efw_phasef_read_vsvy, time_range, probe=probe
    copy_data, prefix+'efw_vsvy', vsvy_var
    interp_time, vsvy_var, common_times
    vsvy = get_var_data(vsvy_var)

    fillval = !values.f_nan
;    flag_var = prefix+'efw_flags'
;    if check_if_update(flag_var, time_range) then begin
;        rbsp_efw_read_flags, time_range, probe=probe
;    endif
;    flags = get_var_data(flag_var, at=common_times, limits=lim)
;    flag_names = lim.labels

    foreach boom_pair, boom_pairs do begin
        ii = fix(strmid(boom_pair,0,1))-1
        jj = fix(strmid(boom_pair,1,1))-1
        sumpair = total(vsvy[*,[ii,jj]],2)*0.5
        vsc_var = prefix+'sc_potential_'+boom_pair
        store_data, vsc_var, common_times, sumpair

        ; Calculate density and remove bad values
        ; Determine density from sc potential.
        density_var = prefix+'density_'+boom_pair
        rbsp_efw_density_fit_from_uh_line, vsc_var, probe, newname=density_var, $
            dmin=dmin, dmax=dmax, setval=fillval

        ; When rbsp_efw_density_fit_from_uh_line does not work, the times are all 0.
        get_data, density_var, times
        if min(times) eq 0 and max(times) eq 0 then begin
            store_data, density_var, common_times, fltarr(ntime)+fillval
        endif

;        key = 'charging_'+boom_pair
;        index = (where(flag_names eq key))[0]
;        charging_flag = flags[*,index]
;        goo = where(charging_flag eq 1, count)
;        if count ne 0 then begin
;
;            get_data, density_var, c_genommon_times, dens, filename=file
;            dens[goo] = fillval
;            store_data, density_var, common_times, dens
;        endif
    endforeach


;---Save data.
    path = fgetpath(file)
    if file_test(path,/directory) eq 0 then file_mkdir, path
    data_file = file
    if file_test(data_file) eq 1 then file_delete, data_file  ; overwrite old files.

    ginfo = dictionary($
        'TITLE', 'RBSP EFW density calibrated according to the upper hybrid line', $
        'TEXT', 'Generated by Sheng Tian at the University of Minnesota, adopted from rbsp_efw_get_flag_values' )
    cdf_save_setting, ginfo, filename=file
    save_vars = prefix+'density_'+boom_pairs
    foreach save_var, save_vars do begin
        get_data, save_var, times, data
        store_data, save_var, times, float(data), limits={units:'cm!U-3'}
    endforeach
    stplot2cdf, save_vars, istp=1, filename=file, time_var='epoch'

end

;time_range = time_double(['2013-01-01','2013-01-02'])
;probe = 'a'
;file = join_path([homedir(),'test.cdf'])
;if file_test(file) eq 1 then file_delete, file
;tic
;rbsp_efw_phasef_read_density_gen_file, time_range, probe=probe, filename=file
;toc

stop
probes = ['a','b']
root_dir = join_path([diskdir('data'),'rbsp'])
secofday = constant('secofday')
foreach probe, probes do begin
    prefix = 'rbsp'+probe+'_'
    rbspx = 'rbsp'+probe
    time_range = rbsp_efw_phasef_get_valid_range('density_uh', probe=probe)
    days = make_bins(time_range+[0,-1]*secofday, secofday)
    foreach day, days do begin
        str_year = time_string(day,tformat='YYYY')
        path = join_path([root_dir,'efw_phasef','density_uh_v02',rbspx,str_year])
        base = prefix+'efw_density_uh_'+time_string(day,tformat='YYYY_MMDD')+'_v02.cdf'
        file = join_path([path,base])
;if file_test(file) eq 1 then continue
        rbsp_efw_phasef_read_density_gen_file, day, probe=probe, filename=file
    endforeach
endforeach


end
