;+ ; Fit E_measure-E_model per maneuver. ; ; time_range. Optional, used because not all data are available. ;- function rbsp_phasef_read_fit_coef, time_range, probe=probe prefix = 'rbsp'+probe+'_' fillval = !values.f_nan xyz = constant('xyz') rbspx = 'rbsp'+probe root_dir = join_path([default_local_root(),'rbsp',rbspx,'e_fit']) if file_test(root_dir) eq 0 then file_mkdir, root_dir ; Load the time of the sections. txt_file = join_path([root_dir,prefix+'maneuver_list_for_efit.txt']) section_times = rbsp_phasef_read_fit_times(probe=probe, filename=txt_file) nsection = n_elements(section_times)*0.5 if n_elements(time_range) eq 0 then time_range = minmax(section_times) ; The file to save data. sav_file = join_path([root_dir, prefix+'efit_coef.sav']) if file_test(sav_file) eq 0 then begin fit_list = list() for section_id=0, nsection-1 do begin section_time_range = reform(section_times[section_id,*]) fit_info = dictionary() fit_info.time_range = section_time_range fit_list.add, fit_info endfor save, fit_list, filename=sav_file endif restore, filename=sav_file ;---Loop through the sections to see if update is needed. is_updated = 0 foreach fit_info, fit_list do begin section_time_range = fit_info.time_range if section_time_range[0] lt time_range[0] then continue if section_time_range[1] gt time_range[1] then continue if fit_info.haskey('y') and fit_info.haskey('z') then continue ;---Load data. rbsp_efw_phasef_prepare_residue_removal, section_time_range, probe=probe dis = snorm(get_var_data(prefix+'r_mgse', times=times)) store_data, prefix+'dis', times, dis orbit_time_step = total(times[0:1]*[-1,1]) index = where(dis le 2) perigee_times = times[time_to_range(index, time_step=1)] store_data, prefix+'perigee_times', 0, perigee_times perigee_times = get_var_data(prefix+'perigee_times') nfit_section = 10 index = where(perigee_times[*,0] ge section_time_range[0] and $ perigee_times[*,1] le section_time_range[1], count) if count eq 0 or count lt nfit_section+1 then continue ;---Prepare fitting. rbsp_read_e_model, section_time_range, probe=probe, datatype='e_model_related' b_mgse = get_var_data(prefix+'b_mgse', times=times) v_mgse = get_var_data(prefix+'v_mgse') vcoro_mgse = get_var_data(prefix+'vcoro_mgse') u_mgse = (v_mgse-vcoro_mgse)*1e-3 ndim = 3 nrec = n_elements(times) xxs = fltarr(ndim,nrec,ndim) ; For Ey. xxs[0,*,1] = u_mgse[*,0]*b_mgse[*,1] xxs[1,*,1] = -u_mgse[*,0]*b_mgse[*,0]-u_mgse[*,2]*b_mgse[*,2] xxs[2,*,1] = u_mgse[*,2]*b_mgse[*,1] ; For Ez. xxs[0,*,2] = u_mgse[*,0]*b_mgse[*,2] xxs[1,*,2] = u_mgse[*,1]*b_mgse[*,0] xxs[2,*,2] = -u_mgse[*,0]*b_mgse[*,0]-u_mgse[*,1]*b_mgse[*,1] rbsp_efw_phasef_read_e_uvw, section_time_range, probe=probe interp_time, prefix+'e_uvw', times e_uvw = get_var_data(prefix+'e_uvw') e_mgse = cotran(e_uvw, times, 'uvw2mgse', probe=probe) ; ; Remove eclipse. ; pad_time = 0. ; rbsp_read_eclipse_flag, section_time_range, probe=probe ; flag = get_var_data(prefix+'eclipse_flag', times=uts) ; index = where(flag eq 1, count) ; fillval = !values.f_nan ; if count ne 0 then begin ; nan_times = uts[time_to_range(index,time_step=1)] ; nnan_time = n_elements(nan_times)*0.5 ; for ii=0,nnan_time-1 do begin ; index = lazy_where(times,'[]', nan_times[ii,*]+[-1,1]*pad_time, count=count) ; if count eq 0 then continue ; e_mgse[index,*] = fillval ; endfor ; endif store_data, prefix+'de_mgse', times, e_mgse yys = get_var_data(prefix+'de_mgse') store_data, prefix+'fit_data', times, xxs, yys ;---Do fit. ; Find the last several perigees. fit_time_ranges = perigee_times[index[count-1-nfit_section:count-1],*] nfit_time_range = n_elements(fit_time_ranges)*0.5 get_data, prefix+'fit_data', times, xxs, yys time_index = [] for jj=0,nfit_time_range-1 do begin index = lazy_where(times, '[]', fit_time_ranges[jj,*], count=count) if count eq 0 then continue time_index = [time_index,index] endfor xxs = xxs[*,time_index,*] yys = yys[time_index,*] uts = times[time_index] ; Calc fit coefs. fit_coef = fltarr(ndim+1,ndim)+fillval fit_index = [1,2] foreach jj, fit_index do begin yy = yys[*,jj] xx = xxs[*,*,jj] index = where(finite(yy) and finite(snorm(transpose(xx))), count) if count lt 10 then continue res = regress(xx[*,index],yy[index], sigma=sigma, const=const) fit_coef[*,jj] = [res[*],const] yfit = reform(xxs[*,*,jj] ## res)+const store_data, prefix+'de'+xyz[jj]+'_fit', uts, [[yy],[yfit],[yy-yfit]], $ limits={colors:sgcolor(['red','green','blue']),labels:['old','fit','new']} endforeach ; Save fit coefs. foreach jj, fit_index do fit_info[xyz[jj]] = reform(fit_coef[*,jj]) is_updated = 1 endforeach if is_updated then begin save, fit_list, filename=sav_file log_file = join_path([data_dir, prefix+'maneuver_fit_coef.txt']) tab = ' ' if file_test(log_file) eq 1 then file_delete, log_file ftouch, log_file foreach fit_info, fit_list do begin msg = '' msg += strjoin(time_string(fit_info.time_range,tformat='YYYY-MM-DD/hh:mm:ss'),' ')+tab the_y = fit_info.haskey('y')? fit_info['y']: 0 the_z = fit_info.haskey('z')? fit_info['z']: 0 msg += strjoin(string(the_y,format='(F13.10)'),' ')+tab msg += strjoin(string(the_z,format='(F13.10)'),' ')+tab lprmsg, msg, log_file endforeach endif good_fit_list = list() foreach fit_info, fit_list do begin if ~fit_info.haskey('y') then continue if ~fit_info.haskey('z') then continue good_fit_list.add, fit_info endforeach return, good_fit_list end foreach probe, ['a','b'] do fit_list = rbsp_phasef_read_fit_coef(probe=probe) stop time_range = time_double(['2018-01-01','2019-01-01']) time_range = time_double(['2013-01-01','2019-01-01']) foreach probe, ['a','b'] do fit_list = rbsp_phasef_read_fit_coef(time_range, probe=probe) end