;+
; Read orbit number.
;
; time. A string or double (unix time) for the wanted date.
; probe=. A string 'a' or 'b'.
;-

pro rbsp_efw_phasef_read_orbit_num_gen_file, time, probe=probe, filename=file, $
    errmsg=errmsg, log_file=log_file


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

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

    if n_elements(time) eq 0 then begin
        errmsg = handle_error('No input time ...')
        lprmsg, errmsg, log_file
        return
    endif
    if size(time[0],/type) eq 7 then time = time_double(time)

;---Constants and settings.
    secofday = 86400d
    errmsg = ''
    prefix = 'rbsp'+probe+'_'
    rbspx = 'rbsp'+probe
    date = time[0]-(time[0] mod secofday)
    time_range = date+[0,secofday]
    time_step = 60
    common_times = make_bins(time_range, time_step)

    ; Read one more orbit to ensure we catch the full perigee.
    rbsp_read_orbit, time_range+[-9,0]*3600, probe=probe, coord='gse'
    r_var = prefix+'r_gse'
    uniform_time, r_var, time_step
    get_data, r_var, times, r_gse
    dis = snorm(r_gse)
    index = where(dis le 2, count)
    if count eq 0 then message, 'Invalid orbit data ...'
    perigee_time_ranges = times[time_to_range(index,time_step=1)]
    nperigee = n_elements(perigee_time_ranges)*0.5
    perigee_times = dblarr(nperigee)
    for perigee_id=0, nperigee-1 do begin
        index = lazy_where(times, '[]', perigee_time_ranges[perigee_id,*])
        min_dis = min(dis[index], min_index)
        perigee_times[perigee_id] = (times[index])[min_index]
    endfor
    index = where(perigee_times ge time_range[0], nperigee_time)
    if nperigee_time eq 0 then message, 'No perigee found ...'
    perigee_times = perigee_times[index]

    ; Orbit num increases when passing a perigee.
    ncommon_time = n_elements(common_times)
    orbit_nums = fltarr(ncommon_time)
    nperigee_ut = n_elements(perigee_times)
    for ii=0, nperigee_ut-1 do begin
        index = where(common_times gt perigee_times[ii], count)
        if count eq 0 then continue
        orbit_nums[index] += 1
    endfor
    
    ; Find the orbit num at the beginning of the day.
    orbnum_var = prefix+'orbit_num'
    if date eq time_double('2012-09-05') then begin
        orbit_num0 = 15
    endif else begin
        rbsp_efw_phasef_read_orbit_num, time_range-secofday, probe=probe
        orbit_num0 = get_var_data(orbnum_var, at=time_range[0])
    endelse
    
    orbit_nums += orbit_num0
    store_data, orbnum_var, common_times, orbit_nums


;---Save data.
    save_var = orbnum_var
    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 orbit number derived from perigee times', $
        'TEXT', 'Generated by Sheng Tian at the University of Minnesota' )
    cdf_save_setting, ginfo, filename=file
    get_data, save_var, times, data
    store_data, save_var, times, float(data), limits={units:'#'}
    stplot2cdf, save_var, istp=1, filename=file, time_var='epoch'

end


probe = 'a'
date = '2012-01-01'
date = '2012-09-25'
;date = '2012-09-05'
date = '2019-01-13'
date = '2016-01-01'
; First day.
date = '2012-09-05'
probe = 'b'
; A day when perigee time is at the beginning of the day.
date = '2012-10-07'
probe = 'b'

;; A day when perigee time is around the beginning of the day.
;date = '2012-09-08'
;probe = 'a'

file = join_path([homedir(),'test.cdf'])
rbsp_efw_phasef_read_orbit_num_gen_file, date, probe=probe, filename=file
end
