@codes/SUB_interpolate_vertical @codes/matchup/SUB_plot_profiles_single_height pro generate_diff, data_all, HEIGHT_prf, BENDINGANGLE_diff BENDINGANGLE_diff = fltarr(n_elements(data_all), n_elements(HEIGHT_prf))+!values.f_nan for profile_loop = 0, n_elements(data_all)-1 do begin ; if profile_loop ne 5 then continue BENDINGANGLE = *(data_all[profile_loop].BENDINGANGLE) HEIGHT = *(data_all[profile_loop].HEIGHT_ba) ; --- remove all the obvious bad files, with large negative BA --- range_id = where(height.data le 45, range_num) if range_num eq 0 then continue BENDINGANGLE_data=BENDINGANGLE.data BENDINGANGLE_data=BENDINGANGLE_data[range_id] if total(~finite(BENDINGANGLE_data)) gt range_num * 0.2 or $ total(BENDINGANGLE_data lt 0) gt 0 then begin print, 'bad:', total(~finite(BENDINGANGLE_data) or BENDINGANGLE_data lt 0), range_num,'|', profile_loop, n_elements(data_all) ; stop continue endif ; --- remove end --- ; --- remove all the obvious bad files, with large positive BA --- range_id = where(height.data ge 35 and height.data lt 45, range_num) if range_num eq 0 then continue BENDINGANGLE_data=BENDINGANGLE.data BENDINGANGLE_data=BENDINGANGLE_data[range_id] if total(BENDINGANGLE_data gt 0.001) gt range_num * 0.2 then begin print, 'bad2:', total(BENDINGANGLE_data gt 0.001), range_num,'|', profile_loop, n_elements(data_all) ; stop continue endif ; --- remove end --- BENDINGANGLE_fm = *(data_all[profile_loop].BENDINGANGLE_fm) HEIGHT_fm = *(data_all[profile_loop].HEIGHT_ba_fm) ; interpolate high resolution to low resolution status = SUB_interpolate_vertical( alog10(BENDINGANGLE.data * 1D6), HEIGHT.data, HEIGHT_fm.data, variable_1) if n_elements(variable_1) lt 10 then continue BENDINGANGLE_interpol = 10 ^ (variable_1 - 6) fractional_diff = 100*(BENDINGANGLE_interpol - BENDINGANGLE_fm.data)/BENDINGANGLE_fm.data ; print, fractional_diff ; stop status = SUB_interpolate_vertical( fractional_diff, HEIGHT_fm.data, HEIGHT_prf, variable_2) ;print, variable_2 ;stop ; SUB_plot_profiles_single_height, 1e3*BENDINGANGLE.data, HEIGHT.data, $ ; 1e3*BENDINGANGLE_fm.data, HEIGHT_fm.data, $ ; variable_2, HEIGHT_prf, $ ; 'BENDINGANGLE', $ ; 'test', $ ; 'codes/test_height_'+string(profile_loop,format='(i3.3)')+'.png' ; stop if n_elements(variable_2) lt 10 then continue BENDINGANGLE_diff[profile_loop, *] = variable_2 endfor end ;------------------------------------ ; .r codes/SUB_generate_fm_std_monthly ; SUB_generate_fm_std_monthly, 'cosmic2', 'ROPP' ; SUB_generate_fm_std_monthly, 'cosmic2', 'UMD' ; SUB_generate_fm_std_monthly, 'cicero', 'ROPP' ; SUB_generate_fm_std_monthly, 'cicero', 'UMD' ; SUB_generate_fm_std_monthly, 'spire', 'ROPP' <--- 2021-10 not done yet ; SUB_generate_fm_std_monthly, 'paz', 'ROPP' ; SUB_generate_fm_std_monthly, 'kompsat5', 'ROPP' ; SUB_generate_fm_std_monthly, 'cosmic', 'ROPP' ; SUB_generate_fm_std_monthly, 'metopa', 'ROPP' ; SUB_generate_fm_std_monthly, 'metopb', 'ROPP' <-- std too large, use metopa ; SUB_generate_fm_std_monthly, 'metopc', 'ROPP' <-- std too large, use metopa ; SUB_generate_fm_std_monthly, 'tsx', 'ROPP' ; SUB_generate_fm_std_monthly, 'tdx', 'ROPP' <-- std too large, use tsx pro SUB_generate_fm_std_monthly, mission, publisher ; ./ropp/plot_profiles_monthly_std.pro if mission eq 'cosmic2' then year_month='2019-10' $ ; ROPP/UMD else if mission eq 'cicero' then year_month='2021-04' $ ; ROPP/UMD else if mission eq 'spire' then year_month='2021-10' $ ; ROPP else if mission eq 'paz' then year_month='2019-10' $ ; ROPP else if mission eq 'kompsat5' then year_month='2015-10' $ ; ROPP else if mission eq 'metopa' then year_month='2021-02' $ ; ROPP else if mission eq 'metopb' then year_month='2021-02' $ ; ROPP else if mission eq 'metopc' then year_month='2021-02' $ ; ROPP else if mission eq 'tsx' then year_month='2021-02' $ ; ROPP else if mission eq 'tdx' then year_month='2021-02' $ ; ROPP else if mission eq 'cosmic' then year_month='2014-05' ; ROPP if mission eq 'paz' then mission2='PNOTS' $ else if mission eq 'kompsat5' then mission2='KOMPSAT' $ else if mission eq 'metopa' then mission2='GRAS' $ else if mission eq 'metopb' then mission2='GRAS' $ else if mission eq 'metopc' then mission2='GRAS' $ else if mission eq 'tsx' then mission2='IGOR' $ else if mission eq 'tdx' then mission2='IGOR' $ else mission2=strupcase(mission) ;file_patten = 'input/'+mission2+'/'+publisher+'dry_noref/from_jun/'+mission+'_'+year_month+'-??.sav' file_patten = 'input/'+mission2+'/'+publisher+'dry_noref/'+mission+'_'+year_month+'-??.sav' ;file_patten = 'input/'+mission2+'/'+publisher+'dry_noref/'+mission+'_'+year_month+'-01.sav' print, file_patten file_list = file_search(file_patten, count=count) data_all = [] for file_loop = 0, count-1 do begin print,file_list[file_loop] restore,file_list[file_loop] ;data_all = [data_all, data[where(data.is_bad eq 0, /null)]] data_all = [data_all, data] print, size(data_all) endfor height_min=8. height_max=45. height_interval=0.2 height_num = fix((height_max-height_min)/height_interval)+1 height_prf = findgen(height_num)*height_interval+height_min ; --- remove obvious outlier first --- generate_diff, data_all, HEIGHT_prf, $ BENDINGANGLE_diff ; output std_standard_0 = stddev(BENDINGANGLE_diff, dimension=1, /nan) num_all = n_elements(data_all) ;stop print, 'size(BENDINGANGLE_diff): ',size(BENDINGANGLE_diff) print, 'size(std_standard_0): ',size(std_standard_0) if publisher eq 'ROPP' or publisher eq 'UMD' then begin ; ropp include too many outliers, run generate_diff one more time std_limit=2 ; is_bad = bytarr(num_all)+1 for profile_loop = 0, num_all-1 do begin if total(finite(BENDINGANGLE_diff[profile_loop, *])) gt 0 and $ total(abs(BENDINGANGLE_diff[profile_loop, *]) gt abs(std_limit*std_standard_0)) eq 0 then $ is_bad[profile_loop] = 0 endfor data_all = data_all[where(is_bad eq 0)] ; --- then generate monthly std using non-outlier profiles --- generate_diff, data_all, HEIGHT_prf, BENDINGANGLE_diff std_standard_1 = stddev(BENDINGANGLE_diff, dimension=1, /nan) this_id = where(round(height_prf) eq 25) this_id = median(this_id) this_range = [this_id-1, this_id, this_id+1] STD_STANDARD_2=STD_STANDARD_1 STD_STANDARD_2[this_range]=STD_STANDARD_1[this_id+2] ; for ROPP only ; ---------------------------- ;ID height STD_STANDARD_1 ;97 24.4000 5.68725 ;98 24.6000 5.82798 ;99 24.8000 16.6040 ;100 25.0000 132.846 <----- ;101 25.2000 74.3548 ;102 25.4000 15.5577 ;103 25.6000 11.7343 ;104 25.8000 11.2513 endif else begin std_standard_1 = std_standard_0 std_standard_2 = std_standard_0 endelse num_good = n_elements(data_all) ; STD_STANDARD_2 is the one filename='codes/std_'+mission+'_'+publisher+'_'+year_month+'.sav' save,height_prf,num_all, num_good, std_standard_0, std_standard_1, STD_STANDARD_2, filename=filename ; -------- ; generate plot ; codes/ropp/plot_profiles_monthly_std.pro stop end