;+
; Generate spice product file for a day.
;
; time. A time in UT sec. Only time[0] is used to determine the day.
; probe=. A string of 'a' or 'b'.
; filename=. A string to set the output file name.
;
;-
pro rbsp_efw_phasef_read_spice_gen_file, time, probe=probe, $
    filename=file, errmsg=errmsg

;---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


;---Constants and settings.
    secofday = 86400d
    deg = 180d/!dpi
    rad = !dpi/180
    re = 6378d & re1 = 1d/re
    errmsg = ''
    common_time_step = 1d      ; sec.


    ; Derived settings.
    date = time[0]-(time[0] mod secofday)
    time_range = date+[0,secofday]
    prefix = 'rbsp'+probe+'_'


;---Init file.
    out_dir = fgetpath(file)
    if file_test(out_dir,/directory) eq 0 then file_mkdir, out_dir
    data_file = file
    if file_test(data_file) eq 1 then file_delete, data_file  ; overwrite old files.
    cdf_id = cdf_create(data_file)
    ginfo = dictionary($
        'TITLE', 'RBSP position and boom direction from SPICE kernel', $
        'TEXT', 'Generated by Sheng Tian at the University of Minnesota' )
    cdf_save_setting, ginfo, filename=cdf_id


;---Save common_times.
    ; Need to overlapping with next day for interpolation purpose.
    common_times = make_bins(time_range, common_time_step)
    utname = 'Epoch'
    data = stoepoch(common_times,'unix')
    settings = dictionary($
        'FIELDNAM', 'Epoch', $
        'UNITS', 'ms', $
        'VAR_TYPE', 'support_data' )
    cdf_save_var, utname, value=data, filename=cdf_id, cdf_type='CDF_EPOCH'
    cdf_save_setting, settings, var=utname, filename=cdf_id



;---Load spice kernels for all times.
    defsysv,'!rbsp_spice', exists=flag
    if flag eq 0 then rbsp_load_spice_kernel

    ; Prepare epochs used for q_uvw2gse.
    scid = strupcase(prefix+'science')
    cspice_str2et, time_string(common_times[0], tformat='YYYY-MM-DDThh:mm:ss.ffffff'), epoch0
    epochs = epoch0+common_times-common_times[0]


    ; Get r_coord, v_coord, and q_uvw2coord.
    ;foreach coord, ['gsm','gse'] do begin
    foreach coord, ['gse'] do begin
        cap_coord = strupcase(coord)
        suffix = '_'+coord

        rbsp_load_spice_state, probe=probe, coord=cap_coord, times=common_times, /no_spice_load


        ; r_coord.
        vname = prefix+'r'+suffix
        tplot_rename, prefix+'state_pos_'+coord, vname
        data = get_var_data(vname)*re1
        store_data, vname, common_times, data
        settings = dictionary($
            'FIELDNAM', 'R '+cap_coord, $
            'UNITS', 'Re', $
            'VAR_TYPE', 'data', $
            'DEPEND_0', utname )
        cdf_save_var, vname, value=data, filename=cdf_id
        cdf_save_setting, settings, var=vname, filename=cdf_id

        ; v_coord.
        vname = prefix+'v'+suffix
        tplot_rename, prefix+'state_vel_'+coord, vname
        r_coord = get_var_data(prefix+'r'+suffix)*re
        data = r_coord
        for ii=0,2 do data[*,ii] = deriv(common_times, r_coord[*,ii])
        settings = dictionary($
            'FIELDNAM', 'V '+cap_coord, $
            'UNITS', 'km/s', $
            'VAR_TYPE', 'data', $
            'DEPEND_0', utname )
        cdf_save_var, vname, value=data, filename=cdf_id
        cdf_save_setting, settings, var=vname, filename=cdf_id

        ; q_uvw2coord.
        cspice_pxform, scid, cap_coord, epochs, muvw
        muvw = transpose(muvw)  ; [n,3,3].
        data = mtoq(muvw)
        vname = prefix+'q_uvw2'+coord
        settings = dictionary($
            'FIELDNAM', 'Q UVW2'+cap_coord, $
            'UNITS', '#', $
            'VAR_TYPE', 'data', $
            'DEPEND_0', utname )
        cdf_save_var, vname, value=data, filename=cdf_id
        cdf_save_setting, settings, var=vname, filename=cdf_id
        ; Save to use for DSC.
        store_data, prefix+'wsc_'+coord, common_times, muvw[*,*,2]
    endforeach


    ; Wsc_gse.
    vname = prefix+'wsc_gse'
    data = get_var_data(vname)
    settings = dictionary($
        'FIELDNAM', 'Wsc in GSE', $
        'UNITS', '#', $
        'VAR_TYPE', 'data', $
        'DEPEND_0', utname )
    cdf_save_var, vname, value=data, filename=cdf_id
    cdf_save_setting, settings, var=vname, filename=cdf_id


    ; MLT.
    vname = prefix+'mlt'
    r_gse = get_var_data(prefix+'r_gse')
    data = atan(r_gse[*,1],r_gse[*,0])*deg/15
    settings = dictionary($
        'FIELDNAM', 'MLT', $
        'UNITS', 'hr', $
        'VAR_TYPE', 'data', $
        'DEPEND_0', utname )
    cdf_save_var, vname, value=data, filename=cdf_id
    cdf_save_setting, settings, var=vname, filename=cdf_id

    ; MLat.
    vname = prefix+'mlat'
    cotrans, prefix+'r_gse', prefix+'r_gsm', /GSE2GSM
    cotrans, prefix+'r_gsm', prefix+'r_sm', /GSM2SM
    r_sm = get_var_data(prefix+'r_sm')
    mlat = atan(r_sm[*,2],sqrt(r_sm[*,0]^2+r_sm[*,1]^2))
    data = mlat*deg
    settings = dictionary($
        'FIELDNAM', 'MLat', $
        'UNITS', 'deg', $
        'VAR_TYPE', 'data', $
        'DEPEND_0', utname )
    cdf_save_var, vname, value=data, filename=cdf_id
    cdf_save_setting, settings, var=vname, filename=cdf_id

    ; lshell.
    vname = prefix+'lshell'
    dis = snorm(r_gse)
    data = dis/(cos(mlat)^2)
    settings = dictionary($
        'FIELDNAM', 'L', $
        'UNITS', '#', $
        'VAR_TYPE', 'data', $
        'DEPEND_0', utname )
    cdf_save_var, vname, value=data, filename=cdf_id
    cdf_save_setting, settings, var=vname, filename=cdf_id


    ; UVW to DSC, adopted from rbsp_uvw_to_dsc.
    ; Save sun_phase_[ssha,sshb]. This is the essential info for this rotation.
    ; From sun_phase (sphase) to rotation, check Line 198-204 in rbsp_uvw_to_dsc.pro.
    wsc_gse = get_var_data(prefix+'wsc_gse')
    ncommon_time = n_elements(common_times)
    ndim = 3
    ; The normal to the sun sensor triggering plane, nTsg, is defined as the
    ; cross product between X_GSE and the spin axis direction.
    x_gse = [1d,0,0] ## (dblarr(ncommon_time)+1)
    ntsg = sunitvec(vec_cross(x_gse, wsc_gse))
    ; The triggering sun sensor look direction perpendicular to the spin axis
    ; is defined by the cross product of the spin axis direction and nTsg
    tsg = sunitvec(vec_cross(wsc_gse,ntsg))
    ;foreach the_ssh, ['A','B'] do begin
    foreach the_ssh, ['A'] do begin
        ssh_id = 'RBSP'+strupcase(probe)+'_SSH_'+the_ssh
        cspice_pxform, 'GSE', ssh_id, epochs, gse2ssh
        tsg_ssh = dblarr(ncommon_time,3)
        for ii=0,ncommon_time-1 do tsg_ssh[ii,*] = gse2ssh[*,*,ii] ## tsg[ii,*]
        deg = constant('deg')
        sphase = atan(tsg_ssh[*,1],tsg_ssh[*,0])*deg
        index = where(sphase lt 0, count)
        if count ne 0 then sphase[index] += 360.
        vname = prefix+'sphase_ssh'+strlowcase(the_ssh)
        data = sphase
        settings = dictionary($
            'FIELDNAM', 'sphase of SSH'+the_ssh, $
            'UNITS', 'deg', $
            'VAR_TYPE', 'data', $
            'DEPEND_0', utname )
        cdf_save_var, vname, value=data, filename=cdf_id
        cdf_save_setting, settings, var=vname, filename=cdf_id
    endforeach


    ; Spin period and spin phase.
    timespan, time, 86400d, /seconds
    rbsp_load_state, probe=probe, dt=common_time_step, /no_spice_load

    vname = prefix+'spin_period'
    data = get_var_data(prefix+'spinper')
    settings = dictionary($
        'FIELDNAM', 'Spin period', $
        'UNITS', 'sec', $
        'VAR_TYPE', 'data', $
        'DEPEND_0', utname )
    cdf_save_var, vname, value=data, filename=cdf_id
    cdf_save_setting, settings, var=vname, filename=cdf_id

    vname = prefix+'spin_phase'
    data = get_var_data(prefix+'spinphase')
    settings = dictionary($
        'FIELDNAM', 'Spin phase', $
        'UNITS', 'deg', $
        'VAR_TYPE', 'data', $
        'DEPEND_0', utname )
    cdf_save_var, vname, value=data, filename=cdf_id
    cdf_save_setting, settings, var=vname, filename=cdf_id


    cdf_close, cdf_id

end


time = time_double(['2013-05-01'])
probe = 'b'


;timespan, time, 86400d, /seconds
;rbsp_read_emfisis, time+[0,86400], probe=probe, id='l2%magnetometer'
;rbsp_uvw_to_dsc, probe, prefix+'b_uvw', /debug
;stop

prefix = 'rbsp'+probe+'_'
file = join_path([homedir(),prefix+'spice_products_'+time_string(time,tformat='YYYY_MMDD')+'_v08.cdf'])
rbsp_efw_phasef_read_spice_gen_file, time, probe=probe, $
    filename=file, errmsg=errmsg
end
