diff --git a/dictionary.md b/dictionary.md index 5c6d394..0c13839 100644 --- a/dictionary.md +++ b/dictionary.md @@ -11,7 +11,7 @@ label it with !Q. Many of these keywords are re-used in other wrappers, but since this is the primary eppsilon wrapper they are defined here. -## input data definitions: +## Input data definitions: **folder_name**: **Required** This defines what folder the data live in. For FHD, it should be the top level folder for that run (which contains folders @@ -52,6 +52,19 @@ RTS it is ['xx', 'yy']. (e.g. 'dirty', 'model', 'res'). For FHD the default is all available types, for RTS it is ['res']. +**uvf_input**: This is a flag (valid values are 0/1, default=0) indicating +that the input cubes are uvf cubes, rather than image space cubes. This is only +supported for FHD inputs and only works for single obsid cubes. This is used +most often for simulation testing & exploration. + +**no_evenodd**: This is a flag (valid values are 0/1, default=0) indicating +that only one set of files (rather than both evens and odds) are expected to be +present. If not set, eppsilon will error fairly quickly rather than proceeding +with just one set of files. + + +## Frequencing selecting, flagging, and averaging keywords: + **freq_ch_range**: Specifies what range of frequency channels to calculate the power spectra for. Default is all available channels (this is modified if the coarse_harm_width keyword is used). @@ -60,18 +73,23 @@ power spectra for. Default is all available channels (this is modified if the spectrum. Errors will be generated if there are fewer than 3 unflagged frequencies. +**freq_flag_repeat**: Only used if freq_flags is used. Integer of number of repeats +of freq_flags needed to generate a full flag array. (Length of freq_flags) * freq_flag_repeat +must equal length of frequencies. + **freq_flag_name**: Only used if freq_flags is used. String to use in the output file names to identify the files with the frequency flagging applied. -**uvf_input**: This is a flag (valid values are 0/1, default=0) indicating -that the input cubes are uvf cubes, rather than image space cubes. This is only -supported for FHD inputs and only works for single obsid cubes. This is used -most often for simulation testing & exploration. +**freq_avg_factor**: Integer by which to average along the frequency axis. -**no_evenodd**: This is a flag (valid values are 0/1, default=0) indicating -that only one set of files (rather than both evens and odds) are expected to be -present. If not set, eppsilon will error fairly quickly rather than proceeding -with just one set of files. +**force_even_freqs**: Only used if freq_avg_factor is used. Forces flags to be ignored +while averaging the frequency array, thus resulting in evenly spaced frequencies. + +**freq_avg_bins**: Set of integers assigning each frequency to a bin. Must have the same +length as number of frequencies. Bins must start at zero and be contiguous. + +**freq_bin_name**: Only used if freq_avg_bins is set. String to use in the output file +names to identify files which have been averaged into bins. ## RTS specific keywords: diff --git a/ps_compare/compare_plot_prep.pro b/ps_compare/compare_plot_prep.pro index aac8510..8bb6294 100644 --- a/ps_compare/compare_plot_prep.pro +++ b/ps_compare/compare_plot_prep.pro @@ -1,12 +1,10 @@ pro compare_plot_prep, folder_names, obs_info, $ cube_types, pols, comp_type, compare_files, $ ps_foldernames = ps_foldernames, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, $ - ps_options = ps_options, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ binning_2d_options = binning_2d_options, binning_1d_options = binning_1d_options, $ plot_slices = plot_slices, slice_type = slice_type, fadd_2dbin = fadd_2dbin, $ - freq_ch_range = freq_ch_range, $ plot_filebase = plot_filebase, save_path = save_path, savefilebase = savefilebase, $ axis_type_1d = axis_type_1d, full_compare = full_compare @@ -21,9 +19,9 @@ pro compare_plot_prep, folder_names, obs_info, $ ;; density correction defaults & file naming for 2D & 1D files if n_elements(ps_options) eq 2 then begin if ps_options[0].wt_cutoffs ne ps_options[1].wt_cutoffs and $ - abs(ps_options[1].wt_cutoffs - ps_options[1].wt_cutoffs) le 1e-3 then begin + abs(ps_options[1].wt_cutoffs - ps_options[0].wt_cutoffs) le 1e-3 then begin - ps_options[1].wt_cutoffs = ps_options[1].wt_cutoffs + ps_options[1].wt_cutoffs = ps_options[0].wt_cutoffs endif if (ps_options[0].wt_cutoffs ne ps_options[1].wt_cutoffs $ @@ -36,19 +34,23 @@ pro compare_plot_prep, folder_names, obs_info, $ endif else begin n_wtcuts = 1 endelse + wt_cutoffs = [ps_options[0].wt_cutoffs, ps_options[1].wt_cutoffs] + wt_measures = [ps_options[0].wt_measures, ps_options[1].wt_measures] endif else begin n_wtcuts = 1 + wt_cutoffs = [ps_options.wt_cutoffs] + wt_measures = [ps_options.wt_measures] endelse kperp_density_names = strarr(n_wtcuts) - wh_cutoff0 = where(ps_options.wt_cutoffs eq 0, count_cutoff0, complement = wh_cutoff_n0, $ + wh_cutoff0 = where(wt_cutoffs eq 0, count_cutoff0, complement = wh_cutoff_n0, $ ncomplement = count_cutoff_n0) - wh_std = where(ps_options.wt_cutoffs eq 1 and ps_options.wt_measures eq 'min', count_std) + wh_std = where(wt_cutoffs eq 1 and wt_measures eq 'min', count_std) if count_cutoff0 gt 0 then kperp_density_names[wh_cutoff0] = '_nodensitycorr' if count_cutoff_n0 gt 0 then begin - kperp_density_names[wh_cutoff_n0] = '_kperp_density_' + ps_options[wh_cutoff_n0].wt_measures $ - + '_gt' + number_formatter(ps_options[wh_cutoff_n0].wt_cutoffs) + kperp_density_names[wh_cutoff_n0] = '_kperp_density_' + wt_measures[wh_cutoff_n0] $ + + '_gt' + number_formatter(wt_measures[wh_cutoff_n0]) endif if count_std gt 0 then kperp_density_names[wh_std] = '_dencorr' @@ -64,28 +66,22 @@ pro compare_plot_prep, folder_names, obs_info, $ same_density_tag = '' endelse - if n_elements(uvf_options1) gt 0 then begin - n_uvf_options = 2 - endif else begin - n_uvf_options = 1 - endelse - n_diffs = max([n_elements(obs_info.info_files), n_elements(cube_types), $ - n_elements(pols), n_elements(ps_options), n_uvf_options]) + n_elements(pols), n_elements(ps_options), n_elements(freq_options), n_elements(uvf_options)]) if n_diffs gt 2 then begin - message, 'only 1 or 2 each of [folder_names, ps_foldernames, ' + $ - 'obs_names, cube_types, pols, spec_window_types, wt_cutoffs, ave_removal] allowed' + message, 'only 1 or 2 each of [info_files, cube_types, pols, ps_options, freq_options, uvf_options] allowed' endif if (n_elements(obs_info.info_files) eq 2 or n_elements(ps_options) eq 2 $ - or n_uvf_options eq 2) and n_elements(cube_types) eq 0 $ - and n_elements(pols) eq 0 and n_elements(full_compare) eq 0 then full_compare = 1 - - if keyword_set(full_compare) and n_elements(info_files) eq 1 and n_elements(ps_options) eq 1 $ - and n_uvf_options eq 1 then begin + or n_elements(freq_options) eq 2 or n_elements(uvf_options) eq 2) $ + and n_elements(cube_types) eq 0 and n_elements(pols) eq 0 and n_elements(full_compare) eq 0 $ + then begin + full_compare = 1 + endif + if keyword_set(full_compare) and n_elements(obs_info.info_files) eq 1 and n_elements(ps_options) eq 1 $ + and n_elements(freq_options) eq 1 and n_elements(uvf_options) eq 1 then begin - message, 'full_compare can only be set if one of [folder names, obs names, ' + $ - 'spec windows, wt_cutoffs/measures, ave_removal values, image window name/size] is length 2.' + message, 'full_compare can only be set if one of [info_files, ps_options, freq_options, uvf_options] is length 2.' endif if comp_type eq 'ratio' and keyword_set(full_compare) then comp_type = 'diff_ratio' @@ -97,7 +93,7 @@ pro compare_plot_prep, folder_names, obs_info, $ cube_types = 'res' endelse n_diffs = max([n_elements(obs_info.info_files), n_elements(cube_types), n_elements(pols), $ - n_elements(ps_options), n_uvf_options]) + n_elements(ps_options), n_elements(freq_options), n_elements(uvf_options)]) if n_elements(pols) eq 0 then begin if n_diffs eq 1 then begin @@ -107,16 +103,16 @@ pro compare_plot_prep, folder_names, obs_info, $ endelse endif n_diffs = max([n_elements(obs_info.info_files), n_elements(cube_types), n_elements(pols), $ - n_elements(ps_options), n_uvf_options]) + n_elements(ps_options), n_elements(freq_options), n_elements(uvf_options)]) if n_diffs eq 1 then begin - message, 'at least one of [folder_names, ps_foldernames, ' + $ - 'obs_names, cube_types, pols, spec_window_types, wt_cutoffs, ave_removal] must be a 2 element vector' + message, 'at least one of [info_files, cube_types, pols, ps_options, ' $ + + 'freq_options, uvf_options] must be a 2 element vector' endif case comp_type of 'diff': n_cubes = 1 - 'ratio': n_sets=1 - 'diff_ratio': n_sets=2 + 'ratio': n_sets = 1 + 'diff_ratio': n_sets = 2 'comp_1d': n_cubes = 1 endcase endif else begin @@ -132,11 +128,12 @@ pro compare_plot_prep, folder_names, obs_info, $ endcase endelse - max_file = n_elements(obs_info.info_files)-1 - max_type = n_elements(cube_types)-1 - max_pol = n_elements(pols)-1 + max_file = n_elements(obs_info.info_files) - 1 + max_type = n_elements(cube_types) - 1 + max_pol = n_elements(pols) - 1 max_ps = n_elements(ps_options) - 1 - max_uvf = n_uvf_options - 1 + max_freq = n_elements(freq_options) - 1 + max_uvf = n_elements(uvf_options) - 1 if n_elements(axis_type_1d) eq 0 then axis_type_1d = 'sym_log' @@ -174,17 +171,23 @@ pro compare_plot_prep, folder_names, obs_info, $ 'comp_1d': op_str = ' and ' endcase - uvf_options_use = uvf_options0 + ;; The freq_options struct can have tags added in `fhd_file_setup`. + ;; But that doesn't get propagated back up because here freq_options is in a list. + ;; So we need to make a copy of the struct so it can be updated, and then + ;; re-create the list from the updated structs. Sigh. + freq_options_use0 = create_struct(freq_options[0]) file_struct_arr1 = fhd_file_setup(obs_info.info_files[0], $ - uvf_options = uvf_options_use, ps_options = ps_options[0], $ - freq_ch_range = freq_ch_range) - if n_elements(obs_info.info_files) eq 2 or max_ps eq 1 or max_uvf eq 1 then begin - if max_uvf eq 1 then begin - uvf_options_use = uvf_options1 - endif + uvf_options = uvf_options[0], freq_options = freq_options_use0, ps_options = ps_options[0]) + if max_file eq 1 or max_ps eq 1 or max_freq eq 1 or max_uvf eq 1 then begin + freq_options_use1 = create_struct(freq_options[max_freq]) file_struct_arr2 = fhd_file_setup(obs_info.info_files[max_file], $ - uvf_options = uvf_options_use, ps_options = ps_options[max_ps],$ - freq_ch_range = freq_ch_range) + uvf_options = uvf_options[max_uvf], freq_options = freq_options_use1, $ + ps_options = ps_options[max_ps]) + if max_freq eq 0 then begin + freq_options = list(freq_options_use0) + endif else begin + freq_options = list(freq_options_use0, freq_options_use1) + endelse endif else begin file_struct_arr2 = file_struct_arr1 endelse @@ -229,6 +232,27 @@ pro compare_plot_prep, folder_names, obs_info, $ endelse endelse + ;; get same & different parts of freq_tag to add to plot file name + if n_elements(file_struct_arr2) eq 0 then begin + same_freq_tag = file_struct_arr1[0].freq_tag + endif else begin + if file_struct_arr1[0].freq_tag eq file_struct_arr2[0].freq_tag then begin + same_freq_tag = file_struct_arr1[0].freq_tag + endif else begin + tag_arr1 = strsplit(file_struct_arr1[0].freq_tag, '_',/extract) + tag_arr2 = strsplit(file_struct_arr2[0].freq_tag, '_',/extract) + match_test = strcmp(tag_arr1, tag_arr2) + wh_diff = where(match_test eq 0, count_diff, complement = wh_same, ncomplement = count_same) + + if count_same gt 0 then begin + same_freq_tag = strjoin(tag_arr1[wh_same], '_') + endif + if count_diff gt 0 then begin + diff_freq_tags = [strjoin(tag_arr1[wh_diff], '_'), strjoin(tag_arr2[wh_diff], '_')] + endif + endelse + endelse + ;; get same & different parts of power_tag to add to plot file name if n_elements(file_struct_arr2) eq 0 then begin same_power_tag = file_struct_arr1[0].power_tag @@ -250,23 +274,21 @@ pro compare_plot_prep, folder_names, obs_info, $ endelse endelse - if n_elements(diff_uvf_tags) gt 0 and n_elements(diff_power_tags) gt 0 then begin - plot_options.note = plot_options.note + diff_uvf_tags[0] + diff_power_tags[0] + $ - op_str + diff_power_tags[1] + diff_uvf_tags[1] - endif else begin - if n_elements(diff_uvf_tags) gt 0 then begin - plot_options.note = plot_options.note + diff_uvf_tags[0] + op_str + diff_uvf_tags[1] - endif - if n_elements(diff_power_tags) gt 0 then begin - plot_options.note = plot_options.note + diff_power_tags[0] + op_str + diff_power_tags[1] - endif - endelse - if n_elements(same_uvf_tag) eq 0 then same_uvf_tag = '' + if n_elements(same_freq_tag) eq 0 then same_freq_tag = '' if n_elements(same_power_tag) eq 0 then same_power_tag = '' if n_elements(diff_uvf_tags) eq 0 then diff_uvf_tags = ['', ''] + if n_elements(diff_freq_tags) eq 0 then diff_freq_tags = ['', ''] if n_elements(diff_power_tags) eq 0 then diff_power_tags = ['', ''] + same_tag_parts = same_uvf_tag + same_freq_tag + same_power_tag + diff_tag_parts = [diff_uvf_tags[0] + diff_freq_tags[0] + diff_power_tags[0], $ + diff_uvf_tags[1] + diff_freq_tags[1] + diff_power_tags[1]] + + if max(diff_tag_parts ne ['', '']) gt 0 then begin + plot_options.note += diff_tag_parts[0] + op_str + diff_tag_parts[1] + endif + case comp_type of 'diff': op_str = '_minus_' 'diff_ratio': op_str = '_diffratio_' @@ -279,47 +301,46 @@ pro compare_plot_prep, folder_names, obs_info, $ if n_elements(folder_names) eq 2 and $ folder_names[0] ne folder_names[n_elements(folder_names)-1] then begin - plot_filebase = obs_info.name_same_parts + same_uvf_tag + same_power_tag + same_density_tag +'__' + $ + plot_filebase = obs_info.name_same_parts + same_tag_parts + same_density_tag +'__' + $ obs_info.name_diff_parts[0] + '_' + obs_info.obs_names[0] + $ - diff_uvf_tags[0] + diff_power_tags[0] + density_tags[0] + op_str + $ + diff_tag_parts[0] + density_tags[0] + op_str + $ obs_info.name_diff_parts[1] + '_' + obs_info.obs_names[1] + $ - diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + density_tags[n_wtcuts-1] + diff_tag_parts[1] + density_tags[n_wtcuts-1] endif else begin - plot_filebase = obs_info.folder_basenames[0] + same_uvf_tag + same_power_tag + same_density_tag + '__' + $ - obs_info.obs_names[0] + diff_uvf_tags[0] + diff_power_tags[0] + density_tags[0] + op_str + $ - obs_info.obs_names[max_file] + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + density_tags[n_wtcuts-1] + plot_filebase = obs_info.folder_basenames[0] + same_tag_parts + same_density_tag + '__' + $ + obs_info.obs_names[0] + diff_tag_parts[0] + density_tags[0] + op_str + $ + obs_info.obs_names[max_file] + diff_tag_parts[1] + density_tags[n_wtcuts-1] endelse endif else begin if n_elements(folder_names) eq 1 then begin if n_elements(obs_info.obs_names) gt 1 then begin - plot_filebase = obs_info.folder_basenames[0] + same_uvf_tag + same_power_tag + same_density_tag + $ + plot_filebase = obs_info.folder_basenames[0] + same_tag_parts + same_density_tag + $ '_' + obs_info.obs_names[0] + '_' + cube_types[0] + '_' + pols[0] + $ - diff_uvf_tags[0] + diff_power_tags[0] + density_tags[0] + op_str + $ + diff_tag_parts[0] + density_tags[0] + op_str + $ obs_info.obs_names[0] + '_' + cube_types[max_type] + '_' + pols[max_pol] + $ - diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + density_tags[n_wtcuts-1] + diff_tag_parts[1] + density_tags[n_wtcuts-1] endif else begin - if obs_info.integrated[0] eq 0 then plot_start = obs_info.folder_basenames[0] + '_' + obs_info.obs_names[0] else plot_start = obs_info.fhd_types[0] + if obs_info.integrated[0] eq 0 then begin + plot_start = obs_info.folder_basenames[0] + '_' + obs_info.obs_names[0] + endif else begin + plot_start = obs_info.fhd_types[0] + endelse - plot_filebase = plot_start + same_power_tag + same_uvf_tag + same_density_tag + '_' + $ + plot_filebase = plot_start + same_tag_parts + same_density_tag + '_' + $ cube_types[0] + '_' + pols[0] + $ - diff_uvf_tags[0] + diff_power_tags[0] + density_tags[0] + op_str + $ + diff_tag_parts[0] + density_tags[0] + op_str + $ cube_types[max_type] + '_' + pols[max_pol] + $ - diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + density_tags[n_wtcuts-1] + diff_tag_parts[1] + density_tags[n_wtcuts-1] endelse endif else begin - plot_filebase = obs_info.name_same_parts + same_uvf_tag + same_power_tag + $ + plot_filebase = obs_info.name_same_parts + same_tag_parts + $ same_density_tag + '__' + strjoin([obs_info.name_diff_parts[0], $ cube_types[0], pols[0]], '_') + $ - diff_uvf_tags[0] + diff_power_tags[0] + density_tags[0] + op_str + $ + diff_tag_parts[0] + density_tags[0] + op_str + $ strjoin([obs_info.name_diff_parts[1], cube_types[max_type], pols[max_pol]], '_') + $ - diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + density_tags[n_wtcuts-1] + diff_tag_parts[1] + density_tags[n_wtcuts-1] endelse endelse plot_filebase = plot_filebase[0] @@ -465,6 +486,43 @@ pro compare_plot_prep, folder_names, obs_info, $ if keyword_set(plot_2d_options.plot_wedge_line) then begin z0_freq = 1420.40 ;; MHz + if ( $ + tag_exist(freq_options[0], 'freq_flags') or tag_exist(freq_options[1], 'freq_flags') $ + or tag_exist(freq_options[0], 'freq_ch_range') or tag_exist(freq_options[1], 'freq_ch_range') $ + or tag_exist(freq_options[0], 'freq_avg_bins') or tag_exist(freq_options[1], 'freq_avg_bins') $ + or freq_options[0].freq_avg_factor gt 1 or freq_options[1].freq_avg_factor gt 1 $ + ) then begin + refresh_options = create_refresh_options() + + ;; The freq_options struct can have tags added in `ps_freq_select_avg`. + ;; But that doesn't get propagated back up because here freq_options is in a list + ;; So we need to make a copy of the struct so it can be updated + freq_options_use0 = create_struct(freq_options[0]) + ;; use the `file_check_only` parameter to raise errors if files don't exist or have the wrong flags + ps_freq_select_avg, file_struct_arr1[0], reform(file_struct_arr1[0].n_vis_freq), $ + refresh_options = refresh_options, freq_options = freq_options_use0, $ + ps_options = ps_options[0], /file_check_only + + if max_freq eq 1 then begin + freq_options_use1 = create_struct(freq_options[1]) + ;; use the `file_check_only` parameter to raise errors if files don't exist or have the wrong flags + ps_freq_select_avg, file_struct_arr2[0], reform(file_struct_arr2[0].n_vis_freq), $ + refresh_options = refresh_options, freq_options = freq_options_use1, $ + ps_options = ps_options[max_ps], /file_check_only + if max(abs(freq_options_use0.frequencies - freq_options_use1.frequencies)) gt 1e-14 then begin + message, "frequencies do not match between compared runs" + endif + endif + frequencies = freq_options_use0.frequencies + endif else begin + ;; frequencies can be slightly different depending on where in the pipeline averaging occurs + ;; allow for very small errors when comparing frquencies + if max(abs(file_struct_arr1[0].frequencies - file_struct_arr1[1].frequencies)) gt 1e-14 then begin + message, "frequencies do not match between compared runs" + endif + frequencies = file_struct_arr1[0].frequencies + endelse + freq_use = file_struct_arr1[0].frequencies if n_elements(freq_ch_range) ne 0 then begin if max(freq_ch_range) gt n_elements(freq_use)-1 then message, 'invalid freq_ch_range' @@ -564,76 +622,55 @@ pro compare_plot_prep, folder_names, obs_info, $ if n_elements(savefilebase) eq 0 then begin if n_elements(obs_info.info_files) eq 1 then begin - savefilebase_use = file_struct_arr1[0].general_filebase + same_uvf_tag + $ - same_power_tag + '_' + $ - type_pol1 + diff_uvf_tags[0] + diff_power_tags[0] + $ - op_str + type_pol2 + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + savefilebase_use = file_struct_arr1[0].general_filebase + same_tag_parts + '_' + $ + type_pol1 + diff_tag_parts[0] + op_str + type_pol2 + diff_tag_parts[1] endif else begin if count_diff eq 0 then begin - savefilebase_use = file_struct_arr1[0].general_filebase + same_uvf_tag + $ - same_power_tag + '_' + type_pol1 + savefilebase_use = file_struct_arr1[0].general_filebase + same_tag_parts + '_' + type_pol1 if type_pol1 ne type_pol2 then begin - savefilebase_use = savefilebase_use + diff_uvf_tags[0] + diff_power_tags[0] + $ - op_str + type_pol2 + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + savefilebase_use = savefilebase_use + diff_tag_parts[0] + $ + op_str + type_pol2 + diff_tag_parts[1] endif else begin - if max([diff_uvf_tags ne ['', ''], diff_power_tags ne ['', '']]) eq 1 then begin - savefilebase_use = savefilebase_use + diff_uvf_tags[0] + diff_power_tags[0] + $ - op_str + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + if max(diff_tag_parts ne ['', '']) gt 0 then begin + savefilebase_use += diff_tag_parts[0] + op_str + diff_tag_parts[1] endif endelse endif else begin if count_same gt 0 then begin if type_pol1 ne type_pol2 then begin - savefilebase_use = strjoin(fileparts_1[wh_same], '_') + same_uvf_tag + $ - same_power_tag + '__' + strjoin(fileparts_1[wh_diff]) + '_' + $ - type_pol1 + diff_uvf_tags[0] + diff_power_tags[0] + op_str + $ - strjoin(fileparts_2[wh_diff]) + '_' + type_pol2 + $ - diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + savefilebase_use = strjoin(fileparts_1[wh_same], '_') + same_tag_parts + $ + '__' + strjoin(fileparts_1[wh_diff]) + '_' + type_pol1 + diff_tag_parts[0] + $ + op_str + strjoin(fileparts_2[wh_diff]) + '_' + type_pol2 + diff_tag_parts[1] endif else begin savefilebase_use = strjoin(fileparts_1[wh_same], '_') + '__' + $ strjoin(fileparts_1[wh_diff]) + op_str + $ - strjoin(fileparts_2[wh_diff]) + '__' + type_pol1 + same_uvf_tag + $ - same_power_tag - if max([diff_uvf_tags ne ['', ''], diff_power_tags ne ['', '']]) eq 1 then begin - savefilebase_use = savefilebase_use + '__' + diff_uvf_tags[0] + $ - diff_power_tags[0] + op_str + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + strjoin(fileparts_2[wh_diff]) + '__' + type_pol1 + same_tag_parts + if max(diff_tag_parts ne ['', '']) gt 0 then begin + savefilebase_use += '__' + diff_tag_parts[0] + op_str + diff_tag_parts[1] endif endelse endif else begin if type_pol1 ne type_pol2 then begin savefilebase_use = file_struct_arr1[0].general_filebase + '_' + $ - type_pol_str[0] + same_uvf_tag + $ - same_power_tag + diff_uvf_tags[0] + diff_power_tags[0] + $ + type_pol_str[0] + same_tag_parts + $ + diff_tag_parts[0] + $ op_str + file_struct_arr2[0].general_filebase + '_' + $ - type_pol_str[1] + same_uvf_tag + $ - same_power_tag + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + type_pol_str[1] + same_tag_parts + diff_tag_parts[1] endif else begin savefilebase_use = file_struct_arr1[0].general_filebase + op_str + $ - file_struct_arr2[0].general_filebase + '__' + type_pol1 + same_uvf_tag + $ - same_power_tag + '__' + diff_uvf_tags[0] + $ - diff_power_tags[0] + op_str + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + file_struct_arr2[0].general_filebase + '__' + type_pol1 + same_tag_parts + $ + '__' + diff_tag_parts[0] + op_str + diff_tag_parts[1] endelse endelse endelse endelse endif else begin if type_pol1 eq type_pol2 then begin - savefilebase_use = savefilebase + same_uvf_tag + same_power_tag + '_' + $ - type_pol1 + '__' + diff_uvf_tags[0] + $ - diff_power_tags[0] + op_str + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + savefilebase_use = savefilebase + same_tag_parts + '_' + $ + type_pol1 + '__' + diff_tag_parts[0] + op_str + diff_tag_parts[1] endif else begin - savefilebase_use = savefilebase + same_uvf_tag + same_power_tag + '_' + $ - type_pol1 + diff_uvf_tags[0] + diff_power_tags[0] + op_str + $ - type_pol2 + diff_uvf_tags[n_elements(diff_uvf_tags)-1] + $ - diff_power_tags[n_elements(diff_power_tags)-1] + savefilebase_use = savefilebase + same_tag_parts + '_' + $ + type_pol1 + diff_tag_parts[0] + op_str + type_pol2 + diff_tag_parts[1] endelse endelse endif @@ -834,7 +871,9 @@ pro compare_plot_prep, folder_names, obs_info, $ plot_2d_options = create_plot_2d_options(plot_2d_options = plot_2d_options, $ kperp_plot_range = kperp_plot_range) - endif + endif else begin + kperp_plot_range = plot_2d_options.kperp_plot_range + endelse endfor endfor diff --git a/ps_compare/compare_setup_structures.pro b/ps_compare/compare_setup_structures.pro new file mode 100644 index 0000000..abe19bf --- /dev/null +++ b/ps_compare/compare_setup_structures.pro @@ -0,0 +1,434 @@ +pro compare_setup_structures, folder_names_in, obs_names_in, $ + ps_foldernames = ps_foldernames, version_test = version_test, $ + spec_window_types = spec_window_types, delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = max_uv_lambda, full_image = full_image, $ + image_clip = image_clip, ave_removal = ave_removal, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, std_power = std_power, $ + kz_use = kz_use, kzuse_name = kzuse_name, $ + all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_avg_bins = freq_avg_bins, freq_bin_name = freq_bin_name, $ + freq_flag_repeat = freq_flag_repeat, $ + diff_plot_path = diff_plot_path, diff_save_path = diff_save_path, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + folder_names = folder_names, obs_info = obs_info, uvf_options = uvf_options, $ + freq_options = freq_options, ps_options = ps_options + + + if n_elements(folder_names_in) gt 2 then message, 'only 1 or 2 folder_names allowed' + if n_elements(folder_names_in) eq 0 then message, 'at least 1 folder name must be specified' + if n_elements(obs_names_in) gt 2 then message, 'only 1 or 2 obs_names_in allowed' + + folder_names = get_folder(folder_names_in, loc_name = loc_name, rts = rts, $ + dirty_folder = dirty_folder) + + if keyword_set(version_test) and n_elements(ps_foldername) eq 0 $ + and n_elements(folder_names_in) eq 1 and n_elements(obs_names_in) lt 2 then begin + git_info = git_info(ps_repository_dir()) + ps_foldernames = ['ps_master', 'ps_' + git_info.branch] + endif + + obs_info = ps_filenames(folder_names, obs_names_in, dirty_folder = dirty_folder, $ + exact_obsnames = exact_obsnames, rts = rts, sim = sim, uvf_input = uvf_input, $ + casa = casa, data_subdirs = data_subdirs, ps_foldernames = ps_foldernames, $ + save_paths = save_paths, plot_paths = plot_paths, refresh_info = refresh_info, $ + no_wtvar_rts = no_wtvar_rts) + + if n_elements(diff_plot_path) eq 0 then begin + if n_elements(diff_save_path) gt 0 then begin + diff_plot_path = diff_save_path + path_sep() + 'plots' + path_sep() + endif + endif + + wh_noinfo = where(obs_info.info_files eq '', count_noinfo) + if count_noinfo gt 0 then message, 'Info files are not all present' + + if n_elements(delta_uv_lambda) gt 1 then message, 'only 1 delta_uv_lambda allowed' + + if n_elements(max_uv_lambda) lt 2 and n_elements(full_image) lt 2 $ + and n_elements(image_clip) lt 2 then begin + + uvf_options = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = max_uv_lambda, full_image = full_image, image_clip = image_clip) + + endif else begin + case n_elements(max_uv_lambda) of + 0: + 1: begin + mul0 = max_uv_lambda + mul1 = max_uv_lambda + end + 2: begin + mul0 = max_uv_lambda[0] + mul1 = max_uv_lambda[1] + end + else: message, 'only 1 or 2 max_uv_lambda values allowed' + endcase + + case n_elements(full_image) of + 0: + 1: begin + fi0 = full_image + fi1 = full_image + end + 2: begin + fi0 = full_image[0] + fi1 = full_image[1] + end + endcase + + case n_elements(image_clip) of + 0: + 1: begin + ic0 = image_clip + ic1 = image_clip + end + 2: begin + ic0 = image_clip[0] + ic1 = image_clip[1] + end + else: message, 'only 1 or 2 image_clip values allowed' + endcase + + uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = mul0, full_image = fi0, image_clip = ic0) + uvf_options1 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = mul1, full_image = fi1, image_clip = ic1) + + uvf_options = list(uvf_options0, uvf_options1) + endelse + + if n_elements(freq_flags) gt 0 then begin + if size(freq_flags, /type) eq 11 then begin + ;; this is a list of freq flags, should have 1 or 2 elements + n_freq_flags = n_elements(freq_flags) + if n_freq_flags gt 2 then begin + message, "If freq_flags are specified as a list, the number of elements gives " $ + + "the number of sets of freq flags, which can only be 1 or 2." + endif + + ;; if a single zero is passed in for a freq_flag set, turn it into a null so later logic works + for flag_ind=0, n_freq_flags-1 do begin + if n_elements(freq_flags[flag_ind]) eq 1 then begin + if freq_flags[flag_ind] eq 0 then begin + freq_flags[freq_ind] = !Null + endif + endif + endfor + endif else begin + ;; this is an array, check the dimensionality to determine number of sets + n_flag_dims = size(freq_flags,/n_dim) + if n_flag_dims gt 2 then begin + message, "If freq_flags are specified as an array, it can only be a 1 or 2 dimensional array." + endif + if n_flag_dims eq 2 then begin + n_freq_flags = (size(freq_flags,/dim))[1] + endif else begin + n_freq_flags = 1 + endelse + + if n_freq_flags gt 2 then begin + message, "If freq_flags are specified as an array, the second dimension gives " $ + + "the number of sets of freq flags, which can only be 1 or 2." + endif + ;; In this case, turn it into a list for unified indexing later + if n_freq_flags eq 1 then begin + freq_flags = list(freq_flags) + endif else begin + freq_flags = list(freq_flags[*, 0], freq_flags[*, 1]) + endelse + endelse + endif + + if size(freq_ch_range,/n_dim) lt 2 and n_elements(freq_flags) lt 2 $ + and n_elements(freq_flag_name) lt 2 and n_elements(freq_flag_repeat) lt 2 $ + and n_elements(freq_avg_bins) lt 2 and n_elements(freq_bin_name) lt 2 $ + and n_elements(freq_avg_factor) lt 2 and n_elements(force_even_freqs) lt 2 then begin + + freq_options = create_freq_options( $ + freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, $ + freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, $ + force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, $ + freq_bin_name = freq_bin_name) + endif else begin + + if n_elements(freq_ch_range) gt 0 then begin + case size(freq_ch_range,/n_dim) of + 1: begin + fc0 = freq_ch_range + fc1 = freq_ch_range + end + 2: begin + if (size(freq_ch_range,/dim))[1] ne 2 then begin + message, 'only 1 or 2 sets of freq_ch_range values allowed' + endif + fc0 = freq_ch_range[*, 0] + fc1 = freq_ch_range[*, 1] + end + else: message, 'only 1 or 2 sets of freq_ch_range values allowed' + endcase + endif + if n_elements(freq_flags) gt 0 then begin + case n_elements(freq_flags) of + 1: begin + ff0 = freq_flags + ff1 = freq_flags + end + 2: begin + ff0 = freq_flags[0] + ff1 = freq_flags[1] + end + else: message, 'only 1 or 2 sets of freq_flags values allowed' + endcase + endif + case n_elements(freq_flag_name) of + 0: + 1: begin + fn0 = freq_flag_name + fn1 = freq_flag_name + end + 2: begin + fn0 = freq_flag_name[0] + fn1 = freq_flag_name[1] + end + else: message, 'only 1 or 2 freq_flag_name values allowed' + endcase + + case n_elements(freq_flag_repeat) of + 0: + 1: begin + fr0 = freq_flag_repeat + fr1 = freq_flag_repeat + end + 2: begin + fr0 = freq_flag_repeat[0] + fr1 = freq_flag_repeat[1] + end + else: message, 'only 1 or 2 freq_flag_repeat values allowed' + endcase + + case n_elements(freq_avg_factor) of + 0: + 1: begin + fa0 = freq_avg_factor + fa1 = freq_avg_factor + end + 2: begin + fa0 = freq_avg_factor[0] + fa1 = freq_avg_factor[1] + end + else: message, 'only 1 or 2 freq_avg_factor values allowed' + endcase + + case n_elements(force_even_freqs) of + 0: + 1: begin + ef0 = force_even_freqs + ef1 = force_even_freqs + end + 2: begin + ef0 = force_even_freqs[0] + ef1 = force_even_freqs[1] + end + else: message, 'only 1 or 2 force_even_freqs values allowed' + endcase + + case n_elements(freq_avg_bins) of + 0: + 1: begin + fb0 = freq_avg_bins + fb1 = freq_avg_bins + end + 2: begin + fb0 = freq_avg_bins[0] + fb1 = freq_avg_bins[1] + end + else: message, 'only 1 or 2 freq_avg_bins values allowed' + endcase + + + case n_elements(freq_bin_name) of + 0: + 1: begin + bn0 = freq_bin_name + bn1 = freq_bin_name + end + 2: begin + bn0 = freq_bin_name[0] + bn1 = freq_bin_name[1] + end + else: message, 'only 1 or 2 freq_bin_name values allowed' + endcase + + freq_options1 = create_freq_options( $ + freq_ch_range = fc0, $ + freq_flags = ff0, $ + freq_flag_name = fn0, $ + freq_flag_repeat = fr0, $ + freq_avg_factor = fa0, $ + force_even_freqs = ef0, $ + freq_avg_bins = fb0, $ + freq_bin_name = bn0) + + freq_options2 = create_freq_options( $ + freq_ch_range = fc1, $ + freq_flags = ff1, $ + freq_flag_name = fn1, $ + freq_flag_repeat = fr1, $ + freq_avg_factor = fa1, $ + force_even_freqs = ef1, $ + freq_avg_bins = fb1, $ + freq_bin_name = bn1) + + freq_options = list(freq_options1, freq_options2) + endelse + + if n_elements(ave_removal) lt 2 and n_elements(wt_cutoffs) lt 2 and $ + n_elements(wt_measures) lt 2 and n_elements(spec_window_types) lt 2 and $ + n_elements(freq_dft) lt 2 and n_elements(dft_z_use) lt 2 and $ + n_elements(kz_use) lt 2 and n_elements(kzuse_name) lt 2 and $ + n_elements(std_power) lt 2 then begin + + ps_options = create_ps_options(ave_removal = ave_removal, $ + wt_cutoffs = wt_cutoffs, wt_measures = wt_measures, $ + spec_window_type = spec_window_types, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, $ + kz_use = kz_use, kzuse_name = kzuse_name, $ + std_power = std_power) + + endif else begin + case n_elements(ave_removal) of + 0: + 1: begin + ar0 = ave_removal + ar1 = ave_removal + end + 2: begin + ar0 = ave_removal[0] + ar1 = ave_removal[1] + end + else: message, 'only 1 or 2 ave_removal values allowed' + endcase + + case n_elements(wt_cutoffs) of + 0: + 1: begin + wtc0 = wt_cutoffs + wtc1 = wt_cutoffs + end + 2: begin + wtc0 = wt_cutoffs[0] + wtc1 = wt_cutoffs[1] + end + else: message, 'only 1 or 2 wt_cutoffs allowed' + endcase + + case n_elements(wt_measures) of + 0: + 1: begin + wtm0 = wt_measures + wtm1 = wt_measures + end + 2: begin + wtm0 = wt_measures[0] + wtm1 = wt_measures[1] + end + else: message, 'only 1 or 2 wt_measures allowed' + endcase + + case n_elements(spec_window_types) of + 0: + 1: begin + spw0 = spec_window_types + spw1 = spec_window_types + end + 2: begin + spw0 = spec_window_types[0] + spw1 = spec_window_types[1] + end + else: message, 'only 1 or 2 spec_window_types allowed' + endcase + + case n_elements(freq_dft) of + 0: + 1: begin + dft0 = freq_dft + dft1 = freq_dft + end + 2: begin + dft0 = freq_dft[0] + dft1 = freq_dft[1] + end + else: message, 'only 1 or 2 freq_dft values allowed' + endcase + + case n_elements(dft_z_use) of + 0: + 1: begin + dftz0 = dft_z_use + dftz1 = dft_z_use + end + 2: begin + dftz0 = dft_z_use[0] + dftz1 = dft_z_use[1] + end + else: message, 'only 1 or 2 dft_z_use values allowed' + endcase + + case n_elements(kz_use) of + 0: + 1: begin + kz0 = kz_use + kz1 = kz_use + end + 2: begin + kz0 = kz_use[0] + kz1 = kz_use[1] + end + else: message, 'only 1 or 2 dft_z_use values allowed' + endcase + + case n_elements(kzuse_name) of + 0: + 1: begin + kzn0 = kzuse_name + kzn1 = kzuse_name + end + 2: begin + kzn0 = kzuse_name[0] + kzn1 = kzuse_name[1] + end + else: message, 'only 1 or 2 dft_z_use values allowed' + endcase + + case n_elements(std_power) of + 0: + 1: begin + sp0 = std_power + sp1 = std_power + end + 2: begin + sp0 = std_power[0] + sp1 = std_power[1] + end + else: message, 'only 1 or 2 std_power values allowed' + endcase + + ps_options0 = create_ps_options(ave_removal = ar0, wt_cutoffs = wtc0, $ + wt_measures = wtm0, spec_window_type = spw0, freq_dft = dft0, $ + kz_use = kz0, kzuse_name = kzn0, $ + dft_z_use = dftz0, std_power = sp0) + + ps_options1 = create_ps_options(ave_removal = ar1, wt_cutoffs = wtc1, $ + wt_measures = wtm1, spec_window_type = spw1, freq_dft = dft1, $ + kz_use = kz1, kzuse_name = kzn1, $ + dft_z_use = dftz1, std_power = sp1) + + ps_options = list(ps_options0, ps_options1) + endelse + +end \ No newline at end of file diff --git a/ps_core/make_2d_files.pro b/ps_core/make_2d_files.pro index 123ef3c..65d071b 100644 --- a/ps_core/make_2d_files.pro +++ b/ps_core/make_2d_files.pro @@ -1,13 +1,17 @@ pro make_2d_files, nfiles, savefile_2d, savefile_k0, power_3D, sim_noise_3D, $ new_noise_3d, noise_expval_3d, weights_3d, kx_mpc, ky_mpc, kz_mpc, $ - delay_params, hubble_param, hinv, kperp_lambda_conv, freq_mask, vs_name, $ + delay_params, hubble_param, hinv, kperp_lambda_conv, vs_name, $ vs_mean, t_sys_meas, window_int, git_hashes, wt_ave_power, ave_power, $ ave_weights, wt_ave_power_freq, ave_power_freq, wt_ave_power_uvf, $ ave_power_uvf, uv_pix_area, uv_area, masked_save_items = masked_save_items, $ noise_3D = noise_3D, sim_noise_diff_3D = sim_noise_diff_3D, $ - wt_measure = wt_measure, wt_cutoff = wt_cutoff, freq_flags = freq_flags, $ + wt_measure = wt_measure, wt_cutoff = wt_cutoff, freq_options = freq_options, $ binning_2d_options = binning_2d_options + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask = freq_options.freq_mask + endif + power_rebin = kspace_rebinning_2d(power_3D, kx_mpc, ky_mpc, kz_mpc, $ kperp_edges_mpc, kpar_edges_mpc, noise_expval = noise_expval_3d, $ binned_noise_expval = binned_noise_expval, weights = weights_3d, $ @@ -77,7 +81,7 @@ pro make_2d_files, nfiles, savefile_2d, savefile_k0, power_3D, sim_noise_3D, $ coarse_width = masked_save_items.coarse_width endif - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = savefile_2d, power, noise, sim_noise, sim_noise_diff, weights, $ noise_expval, kperp_edges, kpar_edges, kperp_bin, kpar_bin, $ kperp_lambda_conv, delay_params, hubble_param, freq_mask, vs_name, vs_mean, $ @@ -96,7 +100,7 @@ pro make_2d_files, nfiles, savefile_2d, savefile_k0, power_3D, sim_noise_3D, $ endelse endif else begin - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = savefile_2d, power, noise, sim_noise, sim_noise_diff, weights, $ noise_expval, kperp_edges, kpar_edges, kperp_bin, kpar_bin, $ kperp_lambda_conv, delay_params, hubble_param, freq_mask, vs_name, vs_mean, $ @@ -131,7 +135,7 @@ pro make_2d_files, nfiles, savefile_2d, savefile_k0, power_3D, sim_noise_3D, $ noise, sim_noise_power = sim_noise, sim_noise_diff = sim_noise_diff, $ nfiles = nfiles, hinv = hinv - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = savefile_k0, power, noise, sim_noise, sim_noise_diff, weights, $ noise_expval, k_edges, k_bin, kperp_lambda_conv, delay_params, $ hubble_param, freq_mask, window_int, wt_ave_power, ave_power, ave_weights, $ diff --git a/ps_core/ps_binning.pro b/ps_core/ps_binning.pro index 23e2312..7db1fa8 100644 --- a/ps_core/ps_binning.pro +++ b/ps_core/ps_binning.pro @@ -1,4 +1,4 @@ -pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $ +pro ps_binning, file_struct, sim = sim, freq_options = freq_options, $ ps_options = ps_options, binning_2d_options = binning_2d_options, $ binning_1d_options = binning_1d_options, plot_options = plot_options, $ plot_types = plot_types, $ @@ -10,7 +10,9 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $ nfiles = n_elements(file_struct.datafile) - if n_elements(freq_flags) ne 0 then freq_mask = file_struct.freq_mask + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask = freq_options.freq_mask + endif restore, file_struct.power_savefile @@ -79,11 +81,11 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $ make_2d_files, nfiles, savefile_2d[j], savefile_k0[j], power_3D, sim_noise_3D, $ new_noise_3d, noise_expval_3d, weights_3d, kx_mpc, ky_mpc, kz_mpc, $ - delay_params, hubble_param, plot_options.hinv, kperp_lambda_conv, freq_mask, $ + delay_params, hubble_param, plot_options.hinv, kperp_lambda_conv, $ vs_name, vs_mean, t_sys_meas, window_int, git_hashes, wt_ave_power, ave_power, $ ave_weights, wt_ave_power_freq, ave_power_freq, wt_ave_power_uvf, ave_power_uvf, $ uv_pix_area, uv_area, noise_3D = noise_3D, sim_noise_diff_3D = sim_noise_diff_3D, $ - wt_measure = wt_meas_use, wt_cutoff = wt_cutoff_use, freq_flags = freq_flags, $ + wt_measure = wt_meas_use, wt_cutoff = wt_cutoff_use, freq_options = freq_options, $ binning_2d_options = binning_2d_options endfor @@ -271,12 +273,12 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $ make_2d_files, nfiles, savefile_masked_2d[j,i], savefile_masked_k0[j,i], $ power_3D, sim_noise_3D, new_noise_3d, noise_expval_3d, mask_weights*weights_3d, $ kx_mpc, ky_mpc, kz_mpc, delay_params, hubble_param, plot_options.hinv, kperp_lambda_conv, $ - freq_mask, vs_name, vs_mean, t_sys_meas, window_int, git_hashes, $ + vs_name, vs_mean, t_sys_meas, window_int, git_hashes, $ wt_ave_power, ave_power, ave_weights, wt_ave_power_freq, ave_power_freq, $ wt_ave_power_uvf, ave_power_uvf, uv_pix_area, uv_area, $ masked_save_items = masked_save_items, noise_3D = noise_3D, $ sim_noise_diff_3D = sim_noise_diff_3D, wt_measure = wt_meas_use, $ - wt_cutoff = wt_cutoff_use, freq_flags = freq_flags, $ + wt_cutoff = wt_cutoff_use, freq_options = freq_options, $ binning_2d_options = binning_2d_options ;; must undefine bin_arr_3d so that a new binning is calculated on next loops. @@ -358,7 +360,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $ kx_range_lambda = kpar_binning_1d_options.kx_range_1dave * kperp_lambda_conv ky_range_lambda = kpar_binning_1d_options.ky_range_1dave * kperp_lambda_conv endelse - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = savefile_kpar_power[j], power, noise, sim_noise, sim_noise_diff, $ weights, noise_expval, k_edges, k_bin, kperp_lambda_conv, delay_params, $ hubble_param, freq_mask, kperp_range, kperp_range_lambda, kx_range, $ @@ -454,7 +456,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $ ky_range_lambda = kperp_binning_1d_options.ky_range_1dave * kperp_lambda_conv endelse - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = savefile_kperp_power[j], power, noise, sim_noise, sim_noise_diff, $ weights, noise_expval, k_edges, k_bin, kperp_lambda_conv, delay_params, $ hubble_param, freq_mask, kperp_range, kperp_range_lambda, kx_range, $ diff --git a/ps_core/ps_freq_select_avg.pro b/ps_core/ps_freq_select_avg.pro new file mode 100644 index 0000000..3166f01 --- /dev/null +++ b/ps_core/ps_freq_select_avg.pro @@ -0,0 +1,415 @@ +pro ps_freq_select_avg, file_struct, n_vis_freq, refresh_options = refresh_options, $ + freq_options = freq_options, ps_options = ps_options, $ + metadata_only = metadata_only, file_check_only = file_check_only + + if n_elements(metadata_only) eq 0 then metadata_only = 0 + if n_elements(file_check_only) eq 0 then file_check_only = 0 + + nfiles = n_elements(file_struct.datafile) + + if tag_exist(file_struct, 'beam_int') then begin + beam_int = file_struct.beam_int + endif + + new_n_vis_freq = n_vis_freq + original_freqs = file_struct.frequencies + if tag_exist(freq_options, 'freq_flags') then begin + original_mask = freq_options.freq_mask + endif + if tag_exist(freq_options, 'freq_ch_range') then begin + original_freqs = original_freqs[min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + if nfiles eq 2 then begin + new_n_vis_freq = new_n_vis_freq[*, min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif else begin + new_n_vis_freq = new_n_vis_freq[min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endelse + if tag_exist(freq_options, 'freq_flags') then begin + original_mask = original_mask[min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif + if tag_exist(file_struct, 'beam_int') then begin + if nfiles eq 2 then begin + beam_int = beam_int[*, min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif else begin + beam_int = beam_int[min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endelse + endif + endif + if tag_exist(freq_options, 'freq_flags') then begin + if nfiles eq 2 then begin + mask_use = rebin(reform(original_mask, 1, n_elements(original_freqs)), size(beam_int,/dimension), /sample) + endif else begin + mask_use = original_mask + endelse + new_n_vis_freq = new_n_vis_freq * mask_use + endif + if tag_exist(file_struct, 'beam_int') then begin + if nfiles eq 2 then begin + beam_int = total(beam_int * new_n_vis_freq, 2) / total(new_n_vis_freq, 2) + endif else begin + beam_int = total(beam_int * new_n_vis_freq) / total(new_n_vis_freq) + endelse + endif + + new_n_vis_freq_select = new_n_vis_freq + + ;; average in frequency + if freq_options.freq_avg_factor gt 1 or tag_exist(freq_options, 'freq_avg_bins') then begin + ;; warn about freq_ch_range also being set + if tag_exist(freq_options, 'freq_ch_range') then begin + print, "both freq_avg_factor and freq_ch_range are set. The values in " $ + + "freq_ch_range are being interpreted as the original channel numbers." + endif + if freq_options.freq_avg_factor gt 1 then begin + ;; check that factor divides evenly into number of frequencies + if n_elements(original_freqs) mod freq_options.freq_avg_factor ne 0 then begin + message, "freq_avg_factor must divide evenly into number of frequencies to be " $ + + "averaged (accounting for freq_ch_range if set)." + endif + + avg_n_freqs = n_elements(original_freqs) / freq_options.freq_avg_factor + if tag_exist(freq_options, 'freq_flags') and (freq_options.force_even_freqs eq 0) then begin + print, "both freq_avg_factor and freq_flags are set. The values in " $ + + "freq_flags are being interpreted as the original channel numbers." + ;; compute new frequencies accounting for flagging. + frequencies = fltarr(avg_n_freqs) + for fi=0, avg_n_freqs-1 do begin + this_inds = findgen(freq_options.freq_avg_factor) + (fi * freq_options.freq_avg_factor) + this_freqs = original_freqs[this_inds] + this_mask = original_mask[this_inds] + wh_unflagged = where(this_mask eq 1, count_unflagged) + if count_unflagged gt 0 then begin + ;; take the mean over unflagged channels + frequencies[fi] = mean(this_freqs[wh_unflagged]) + endif else begin + ;; all channels are flagged, take the mean over all channels + frequencies[fi] = mean(this_freqs) + endelse + endfor + endif else begin + ;; do a simple average of the input frequencies. + frequencies = reform(original_freqs, freq_options.freq_avg_factor, avg_n_freqs) + frequencies = mean(frequencies, dimension=1) + endelse + ;; After averaging, if there are any fully flagged frequencies they should remain flagged + if tag_exist(freq_options, 'freq_flags') then begin + new_mask = reform(original_mask, freq_options.freq_avg_factor, avg_n_freqs) + new_mask = max(new_mask, dimension=1) + endif + if nfiles eq 2 then begin + new_n_vis_freq = reform(new_n_vis_freq, nfiles, freq_options.freq_avg_factor, avg_n_freqs) + new_n_vis_freq = total(new_n_vis_freq, 2) + endif else begin + new_n_vis_freq = reform(new_n_vis_freq, freq_options.freq_avg_factor, avg_n_freqs) + new_n_vis_freq = total(new_n_vis_freq, 1) + endelse + + endif else begin ;; averaging with frequency bins + if n_elements(original_freqs) ne n_elements(freq_options.freq_avg_bins) then begin + message, "freq_avg_bins must have length equal to number of frequencies" + endif + h = histogram(freq_options.freq_avg_bins, reverse_indices=ri) + if float(n_elements(h)) ne total(h ne 0) then begin + message, "freq_avg_bins must start at zero and be contiguous" + endif + new_nfreq = n_elements(h) + frequencies = fltarr(new_nfreq) + new_mask = fltarr(new_nfreq) + if nfiles eq 2 then begin + bin_n_vis_freq = fltarr(2,new_nfreq) + endif else begin + bin_n_vis_freq = fltarr(new_nfreq) + endelse + + for fi = 0, new_nfreq - 1 do begin + ; ge t indices for data going into each frequency bin + bin_inds = ri[ri[fi] : ri[fi + 1] - 1] + this_freqs = original_freqs[bin_inds] + if tag_exist(freq_options, 'freq_flags') then begin + this_mask = original_mask[bin_inds] + wh_unflagged = where(this_mask eq 1, count_unflagged) + if count_unflagged gt 0 then begin + ;; take the mean over unflagged channels + frequencies[fi] = mean(this_freqs[wh_unflagged]) + new_mask[fi] = 1 + endif else begin + ;; all channels are flagged, take the mean over all channels + frequencies[fi] = mean(this_freqs) + new_mask[fi] = 0 + endelse + endif else begin + frequencies[fi] = mean(this_freqs) + endelse + if nfiles eq 2 then begin + bin_n_vis_freq[*, fi] = total(new_n_vis_freq[*, bin_inds], 2) + endif else begin + bin_n_vis_freq[fi] = total(new_n_vis_freq[bin_inds]) + endelse + endfor + new_n_vis_freq = bin_n_vis_freq + endelse + + if ( $ + not ps_options.freq_dft $ + and tag_exist(freq_options, 'freq_flags') $ + and freq_options.force_even_freqs eq 0 $ + ) then begin + freq_diff = frequencies - shift(frequencies, 1) + freq_diff = freq_diff[1:*] + if max(abs(freq_diff-freq_diff[0])) gt 1e-12 then begin + ps_options = create_ps_options(ps_options = ps_options, freq_dft = 1) + endif + endif + freq_options = create_freq_options(freq_options = freq_options, $ + frequencies = frequencies, n_vis_freq = new_n_vis_freq, beam_int = beam_int) + if tag_exist(freq_options, 'freq_flags') then begin + freq_options = create_freq_options(freq_options = freq_options, new_freq_mask = new_mask) + endif + endif else begin ;; endif averaging + if tag_exist(freq_options, 'freq_ch_range') then begin + freq_options = create_freq_options(freq_options = freq_options, $ + frequencies = original_freqs, n_vis_freq = new_n_vis_freq, beam_int = beam_int) + if tag_exist(freq_options, 'freq_flags') then begin + freq_options = create_freq_options(freq_options = freq_options, new_freq_mask = original_mask) + endif + endif else begin + ;; get here if there's only flagging but no channel selection or averaging + freq_options = create_freq_options(freq_options = freq_options, $ + frequencies = original_freqs, n_vis_freq = new_n_vis_freq, beam_int = beam_int) + endelse + endelse + + if metadata_only then begin + ;; don't do anything with the actual data cubes + return + endif + + for i=0, nfiles-1 do begin + test_uvf = file_valid(file_struct.uvf_savefile[i]) + if not test_uvf then test_uvf = check_old_path(file_struct, 'uvf_savefile', index=i) + if not test_uvf and file_check_only then begin + message, "the uvf file (" + file_struct.uvf_savefile[i] + ") does not exist." + endif + + test_wt_uvf = file_valid(file_struct.uvf_weight_savefile[i]) + if not test_wt_uvf then test_wt_uvf = check_old_path(file_struct, 'uvf_weight_savefile', index=i) + if not test_wt_uvf and file_check_only then begin + message, "the weight_uvf file (" + file_struct.uvf_weight_savefile[i] + ") does not exist." + endif + + if tag_exist(file_struct, 'beam_savefile') then begin + test_beam = file_valid(file_struct.beam_savefile[i]) + if not test_beam then test_beam = check_old_path(file_struct, 'beam_savefile', index=i) + if not test_beam and file_check_only then begin + message, "the beam file (" + file_struct.beam_savefile[i] + ") does not exist." + endif + endif else test_beam = 1 + + if tag_exist(freq_options, 'freq_flags') then begin + if test_uvf eq 1 then begin + old_freq_mask = getvar_savefile(file_struct.uvf_savefile[i], 'freq_mask') + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then begin + if file_check_only then begin + message, "the uvf file (" + file_struct.uvf_savefile[i] + ") has the wrong freq flagging." + endif + test_uvf = 0 + endif + endif + if test_wt_uvf eq 1 then begin + old_freq_mask = getvar_savefile(file_struct.uvf_weight_savefile[i], 'freq_mask') + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then begin + if file_check_only then begin + message, "the weight_uvf file (" + file_struct.uvf_weight_savefile[i] + ") has the wrong freq flagging." + endif + test_wt_uvf = 0 + endif + endif + if test_beam eq 1 then begin + old_freq_mask = getvar_savefile(file_struct.beam_savefile[i], 'freq_mask') + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then begin + if file_check_only then begin + message, "the beam file (" + file_struct.beam_savefile[i] + ") has the wrong freq flagging." + endif + test_beam = 0 + endif + endif + ;; This is the original freq_mask, which is what we want to save in files + freq_mask = freq_options.freq_mask + endif + + n_freq_orig = n_elements(file_struct.frequencies) + if test_uvf eq 0 or refresh_options.refresh_freq_select_avg then begin + + ;; check for the full cube to avoid redoing the DFT + full_uvf_file = file_struct.uvf_full_savefile[i] + test_full_uvf = file_valid(full_uvf_file) + if not test_full_uvf then test_full_uvf = check_old_path(file_struct, full_uvf_file, index=i) + + if test_full_uvf eq 0 then begin + message, 'full uvf cube does not exist. expected full file: ' + full_uvf_file + endif + restore, full_uvf_file + + if tag_exist(freq_options, 'freq_flags') then begin + data_cube = data_cube * rebin(reform(freq_options.freq_mask, 1, 1, n_freq_orig), $ + size(data_cube, /dimension), /sample) + endif + if tag_exist(freq_options, 'freq_ch_range') then begin + data_cube = data_cube[*, *, min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif + + ;; frequency averaging (really summing because this is just the numerator) + if freq_options.freq_avg_factor gt 1 then begin + data_shape = size(data_cube, /dimensions) + data_cube = reform(data_cube, data_shape[0], data_shape[1], freq_options.freq_avg_factor, avg_n_freqs) + data_cube = total(data_cube, 3) + endif + + if tag_exist(freq_options, 'freq_avg_bins') then begin + data_shape = size(data_cube, /dimensions) + new_data_cube = dcomplexarr(data_shape[0], data_shape[1], new_nfreq) + for fi = 0, new_nfreq - 1 do begin + ; ge t indices for data going into each frequency bin + bin_inds = ri[ri[fi] : ri[fi + 1] - 1] + new_data_cube[*, *, fi] = total(data_cube[*, *, bin_inds], 3) + endfor + data_cube = new_data_cube + endif + + if tag_exist(freq_options, 'freq_flags') then begin + save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, frequencies, $ + data_cube, freq_mask, uvf_git_hash + endif else begin + save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, frequencies, $ + data_cube, uvf_git_hash + endelse + undefine, data_cube, uvf_git_hash + + + endif ;; endif test_uvf eq 0 or refresh + + if test_wt_uvf eq 0 or refresh_options.refresh_freq_select_avg then begin + + ;; if this is a limited freq. range cube, check for the full cube to avoid redoing the DFT + full_uvf_wt_file = file_struct.uvf_weight_full_savefile[i] + test_full_wt_uvf = file_valid(full_uvf_wt_file) + if not test_full_wt_uvf then test_full_wt_uvf = check_old_path(file_struct, 'uvf_weight_savefile', index=i) + + if test_full_wt_uvf eq 0 then begin + message, 'full uvf weight cube does not exist. expected full file: ' + full_uvf_wt_file + endif + restore, full_uvf_wt_file + + if tag_exist(freq_options, 'freq_flags') then begin + flag_arr = rebin(reform(freq_options.freq_mask, 1, 1, n_freq_orig), size(weights_cube, /dimension), /sample) + weights_cube = weights_cube * flag_arr + variance_cube = variance_cube * flag_arr + endif + if tag_exist(freq_options, 'freq_ch_range') then begin + weights_cube = weights_cube[*, *, min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + variance_cube = variance_cube[*, *, min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif + + ;; frequency averaging (really summing because these are just the numerator and denominators) + if freq_options.freq_avg_factor gt 1 then begin + weights_cube = reform(weights_cube, data_shape[0], data_shape[1], freq_options.freq_avg_factor, avg_n_freqs) + weights_cube = total(weights_cube, 3) + variance_cube = reform(variance_cube, data_shape[0], data_shape[1], freq_options.freq_avg_factor, avg_n_freqs) + variance_cube = total(variance_cube, 3) + endif + + if tag_exist(freq_options, 'freq_avg_bins') then begin + data_shape = size(weights_cube, /dimensions) + new_weights_cube = dcomplexarr(data_shape[0], data_shape[1], new_nfreq) + new_variance_cube = dblarr(data_shape[0], data_shape[1], new_nfreq) + for fi = 0, new_nfreq - 1 do begin + if ri[fi + 1] gt ri[fi] then begin + ; ge t indices for data going into each frequency bin + bin_inds = ri[ri[fi] : ri[fi + 1] - 1] + new_weights_cube[*, *, fi] = total(weights_cube[*, *, bin_inds], 3) + new_variance_cube[*, *, fi] = total(variance_cube[*, *, bin_inds], 3) + endif + endfor + weights_cube = new_weights_cube + variance_cube = new_variance_cube + endif + + if tag_exist(freq_options, 'freq_flags') then begin + save, file = file_struct.uvf_weight_savefile[i], kx_rad_vals, $ + ky_rad_vals, frequencies, weights_cube, variance_cube, freq_mask, uvf_wt_git_hash + endif else begin + save, file = file_struct.uvf_weight_savefile[i], kx_rad_vals, $ + ky_rad_vals, frequencies, weights_cube, variance_cube, uvf_wt_git_hash + endelse + undefine, weights_cube, variance_cube, uvf_wt_git_hash + + endif ;; endif test_wt_uvf eq 0 or refresh + + if (test_beam eq 0 or refresh_options.refresh_freq_select_avg) $ + and tag_exist(file_struct, "beam_full_savefile") then begin + + ;; check for the full cube to avoid redoing the pixel selection + full_beam_file = file_struct.beam_full_savefile[i] + test_full_beam = file_valid(full_beam_file) + if not test_full_beam then test_full_beam = check_old_path(file_struct, full_beam_file, index=i) + + if test_full_beam eq 0 then begin + message, 'full beam file does not exist. expected full file: ' + full_beam_file + endif + + pixels_kept = getvar_savefile(file_struct.beam_full_savefile[i], "pixels") + beam_git_hash = getvar_savefile(file_struct.beam_full_savefile[i], "beam_git_hash") + git, repo_path = ps_repository_dir(), result=this_run_git_hash + + ;; also have to get full cube to do frequency selections before averaging. + beam = getvar_savefile(file_struct.beamfile[i], file_struct.beamvar) + if tag_exist(file_struct, 'nside') ne 0 then begin + pixel_nums = getvar_savefile(file_struct.pixelfile[0], file_struct.pixelvar[0]) + endif + + match, pixel_nums, pixels_kept, pixel_nums_inds, subb, count=n_pixels_keep + if n_pixels_keep ne n_elements(pixel_nums) then begin + pixels = pixel_nums[pixel_nums_inds] + beam = beam[pixel_nums_inds, *] + endif else begin + pixels = pixel_nums + endelse + + if tag_exist(freq_options, 'freq_ch_range') then begin + beam = beam[*, *, min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif + + if max(beam) le 1.1 then begin + ;; beam is peak normalized to 1 + print, 'WARNING: It appears that this beam is from a very old run ' + $ + 'of FHD (max of beam^2 is ~1). If that is not the case, there ' + $ + 'might be something unexpected happening' + temp = beam * rebin(reform(new_n_vis_freq_select[i, *], 1, n_freq_orig), n_pixels_keep, n_freq_orig, /sample) + endif else if max(beam) le file_struct.n_obs[i]*1.1 then begin + ;; beam is peak normalized to 1 for each obs, then summed over obs so peak is ~ n_obs + print, 'WARNING: It appears that this beam is from an old run ' + $ + 'of FHD (max of beam^2 is ~n_obs). If that is not the case, there ' + $ + 'might be something unexpected happening' + temp = (beam/file_struct.n_obs[i]) * rebin(reform(new_n_vis_freq_select[i, *], 1, n_freq_orig), $ + n_pixels_keep, n_freq_orig, /sample) + endif else begin + ;; beam is peak normalized to 1 then multiplied by nvis_freq for each obs & summed + temp = beam + if tag_exist(freq_options, 'freq_flags') then begin + ;; original_mask has freq_ch_range applied but not freq averaging + temp *= rebin(reform(original_mask, 1, n_elements(original_mask)), size(beam, /dimension), /sample) + endif + endelse + + avg_beam = total(temp, 2) / total(new_n_vis_freq_select[i, *]) + + nside = file_struct.nside + + beam_git_hash = this_run_git_hash + + save, file=file_struct.beam_savefile[i], avg_beam, pixels, freq_mask, nside, beam_git_hash + + endif ;; endif test_beam eq 0 or refresh and dropping some freqs (otherwise can just use the full average) + + endfor ;; end for loop over files + +end \ No newline at end of file diff --git a/ps_core/ps_image_to_uvf.pro b/ps_core/ps_image_to_uvf.pro index d96979f..5503020 100644 --- a/ps_core/ps_image_to_uvf.pro +++ b/ps_core/ps_image_to_uvf.pro @@ -1,14 +1,10 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, no_var = no_var, $ - refresh_options = refresh_options, uvf_options = uvf_options + no_var = no_var, refresh_options = refresh_options, uvf_options = uvf_options if tag_exist(file_struct, 'nside') ne 0 then healpix = 1 else healpix = 0 nfiles = n_elements(file_struct.datafile) - if n_elements(freq_flags) ne 0 then freq_mask = file_struct.freq_mask - if n_elements(freq_ch_range) ne 0 then begin - n_freq = n_elements(file_struct.frequencies[min(freq_ch_range):max(freq_ch_range)]) - endif else n_freq = n_elements(file_struct.frequencies) + n_freq = n_elements(file_struct.frequencies) datavar = strupcase(file_struct.datavar) if datavar eq '' then begin @@ -19,20 +15,39 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ input_uvf_wtfiles = file_struct.uvf_weight_savefile endif + if tag_exist(file_struct, 'uvf_full_savefile') then begin + uvf_full_tag = '_full' + uvf_savefile = file_struct.uvf_full_savefile + uvf_weight_savefile = file_struct.uvf_weight_full_savefile + endif else begin + uvf_full_tag = '' + uvf_savefile = file_struct.uvf_savefile + uvf_weight_savefile = file_struct.uvf_weight_savefile + endelse + + if tag_exist(file_struct, 'beam_full_savefile') then begin + beam_full_tag = '_full' + beam_savefile = file_struct.beam_full_savefile + endif else begin + beam_full_tag = '' + beam_savefile = file_struct.beam_savefile + endelse for i=0, nfiles-1 do begin - test_uvf = file_valid(file_struct.uvf_savefile[i]) - if not test_uvf then test_uvf = check_old_path(file_struct, 'uvf_savefile', index=i) + + test_uvf = file_valid(uvf_savefile[i]) + if not test_uvf then test_uvf = check_old_path(file_struct, 'uvf' + uvf_full_tag + '_savefile', index=i) - test_wt_uvf = file_valid(file_struct.uvf_weight_savefile[i]) - if not test_wt_uvf then test_wt_uvf = check_old_path(file_struct, 'uvf_weight_savefile', index=i) + test_wt_uvf = file_valid(uvf_weight_savefile[i]) + if not test_wt_uvf then test_wt_uvf = check_old_path(file_struct, 'uvf_weight' + uvf_full_tag + '_savefile', index=i) if tag_exist(file_struct, 'beam_savefile') then begin - test_beam = file_valid(file_struct.beam_savefile[i]) - if not test_beam then test_beam = check_old_path(file_struct, 'beam_savefile', index=i) + test_beam = file_valid(beam_savefile[i]) + if not test_beam then test_beam = check_old_path(file_struct, 'beam' + beam_full_tag + '_savefile', index=i) endif else test_beam = 1 - if test_beam eq 0 then print, "beam_cube files not present under data/beam_cubes" - + if test_beam eq 0 then begin + print, "beam_cube files not present under data/beam_cubes" + endif test_radec_uvf = file_valid(file_struct.radec_file) if not test_radec_uvf then test_radec_uvf = check_old_path(file_struct, 'radec_file') @@ -41,95 +56,6 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ test_radec_uvf_use = 1 endif - if (test_uvf eq 1 or test_wt_uvf eq 1) and n_elements(freq_flags) ne 0 then begin - if test_uvf eq 1 then old_freq_mask = getvar_savefile(file_struct.uvf_savefile[i], 'freq_mask') - if test_wt_uvf eq 1 then old_freq_mask = getvar_savefile(file_struct.uvf_weight_savefile[i], 'freq_mask') - if total(abs(old_freq_mask - freq_mask)) ne 0 then test_uvf = 0 - endif - - if test_uvf eq 0 and not refresh_options.refresh_dft and $ - (n_elements(freq_ch_range) ne 0 or n_elements(freq_flags) ne 0) then begin - - ;; if this is a limited freq. range cube, check for the full cube to avoid redoing the DFT - full_uvf_file = file_struct.uvf_savefile[i] - if n_elements(freq_ch_range) ne 0 then begin - full_uvf_file = strjoin(strsplit(full_uvf_file, '_ch[0-9]+-[0-9]+', /regex, /extract)) - endif - if n_elements(freq_flags) ne 0 then begin - full_uvf_file = strjoin(strsplit(full_uvf_file, '_flag+', /regex, /extract)) - endif - test_full_uvf = file_valid(full_uvf_file) - if not test_full_uvf then test_full_uvf = check_old_path(file_struct, 'uvf_savefile', index=i) - - if test_full_uvf eq 1 then begin - restore, full_uvf_file - - if n_elements(freq_ch_range) ne 0 then begin - data_cube = data_cube[*, *, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - data_cube = data_cube * rebin(reform(freq_mask, 1, 1, n_freq), $ - size(data_cube, /dimension), /sample) - endif - - if n_elements(freq_flags) gt 0 then begin - save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, $ - data_cube, freq_mask, uvf_git_hash - endif else begin - save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, $ - data_cube, uvf_git_hash - endelse - undefine, data_cube, uvf_git_hash - - test_uvf = 1 - endif ;; endifif full uvf exists - - endif ;; endif limited freq range & no matching uvf - - if test_wt_uvf eq 0 and not refresh_options.refresh_weight_dft and $ - (n_elements(freq_ch_range) ne 0 or n_elements(freq_flags) ne 0) then begin - - ;; if this is a limited freq. range cube, check for the full cube to avoid redoing the DFT - full_uvf_wt_file = file_struct.uvf_weight_savefile[i] - if n_elements(freq_ch_range) ne 0 then begin - full_uvf_wt_file = strjoin(strsplit(full_uvf_wt_file, '_ch[0-9]+-[0-9]+', $ - /regex, /extract)) - endif - if n_elements(freq_flags) ne 0 then begin - full_uvf_wt_file = strjoin(strsplit(full_uvf_wt_file, '_flag+', $ - /regex, /extract)) - endif - test_full_wt_uvf = file_valid(full_uvf_wt_file) - if not test_full_wt_uvf then test_full_wt_uvf = check_old_path(file_struct, 'uvf_weight_savefile', index=i) - - if test_full_wt_uvf eq 1 then begin - restore, full_uvf_wt_file - - if n_elements(freq_ch_range) ne 0 then begin - weights_cube = weights_cube[*, *, min(freq_ch_range):max(freq_ch_range)] - variance_cube = variance_cube[*, *, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - flag_arr = rebin(reform(freq_mask, 1, 1, n_freq), $ - size(weights_cube, /dimension), /sample) - weights_cube = weights_cube * flag_arr - variance_cube = variance_cube * flag_arr - endif - - if n_elements(freq_flags) gt 0 then begin - save, file = file_struct.uvf_weight_savefile[i], kx_rad_vals, $ - ky_rad_vals, weights_cube, variance_cube, freq_mask, uvf_wt_git_hash - endif else begin - save, file = file_struct.uvf_weight_savefile[i], kx_rad_vals, $ - ky_rad_vals, weights_cube, variance_cube, uvf_wt_git_hash - endelse - undefine, weights_cube, variance_cube, uvf_wt_git_hash - - test_wt_uvf = 1 - endif ;; endif full wt_uvf exists - - endif ;; endif limited freq range & no matching wt_uvf - if test_uvf eq 0 or test_wt_uvf eq 0 or test_beam eq 0 or test_radec_uvf_use eq 0 $ or refresh_options.refresh_dft or refresh_options.refresh_weight_dft $ or refresh_options.refresh_beam then begin @@ -144,52 +70,7 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ test_input_uvf[1] = file_valid(input_uvf_files[i,1]) if not test_input_uvf[0] then test_input_uvf[0] = check_old_path(file_struct, 'derived_uvf_inputfiles', index=nfiles+i) - if min(test_input_uvf) eq 1 and (n_elements(freq_ch_range) ne 0 $ - or n_elements(freq_flags) ne 0) then begin - - for j=0, 1 do begin - if test_input_uvf[j] eq 0 then begin - ;; this is a limited freq. range cube, check for the full cube to - ;; avoid redoing the DFT - full_uvf_file = input_uvf_files[i,j] - if n_elements(freq_ch_range) ne 0 then begin - full_uvf_file = strjoin(strsplit(full_uvf_file, '_ch[0-9]+-[0-9]+', $ - /regex, /extract)) - endif - if n_elements(freq_flags) ne 0 then begin - full_uvf_file = strjoin(strsplit(full_uvf_file, '_flag[a-z0-9]+', $ - /regex, /extract)) - endif - test_full_uvf = file_valid(full_uvf_file) - if not test_full_uvf then test_full_uvf = check_old_path(file_struct, 'derived_uvf_inputfiles', index=j*nfiles+i) - - if test_full_uvf eq 1 then begin - restore, full_uvf_file - - if n_elements(freq_flags) ne 0 then begin - data_cube = data_cube * rebin(reform(freq_mask, 1, 1, n_freq), $ - size(data_cube, /dimension), /sample) - endif - if n_elements(freq_ch_range) ne 0 then begin - data_cube = data_cube[*, *, min(freq_ch_range):max(freq_ch_range)] - endif - - if n_elements(freq_flags) gt 0 then begin - save, file = input_uvf_files[i,j], kx_rad_vals, ky_rad_vals, $ - data_cube, freq_mask, uvf_git_hash - endif else begin - save, file = input_uvf_files[i,j], kx_rad_vals, ky_rad_vals, $ - data_cube, uvf_git_hash - endelse - undefine, data_cube, uvf_git_hash - - test_uvf = 1 - - endif ;; endif test_full_uvf eq 1 - endif;; endif test_input_uvf[j] eq 1 - endfor;; end loop over 2 derived cubes - - endif else if min(test_input_uvf) eq 0 then begin + if min(test_input_uvf) eq 0 then begin message, 'derived cube but input_uvf cube files do not exist' endif @@ -201,21 +82,6 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ kx_rad_vals = getvar_savefile(input_uvf_files[i,1], 'kx_rad_vals') ky_rad_vals = getvar_savefile(input_uvf_files[i,1], 'ky_rad_vals') - if n_elements(freq_flags) ne 0 then begin - dirty_freq_mask = getvar_savefile(input_uvf_files[i,0], 'freq_mask') - model_freq_mask = getvar_savefile(input_uvf_files[i,1], 'freq_mask') - if n_elements(freq_ch_range) ne 0 then begin - dirty_freq_mask = dirty_freq_mask[min(freq_ch_range):max(freq_ch_range)] - model_freq_mask = model_freq_mask[min(freq_ch_range):max(freq_ch_range)] - endif - if total(abs(dirty_freq_mask - freq_mask)) ne 0 then begin - message, 'freq_mask of dirty file does not match current freq_mask' - endif - if total(abs(model_freq_mask - freq_mask)) ne 0 then begin - message, 'freq_mask of model file does not match current freq_mask' - endif - endif - if total(abs(kx_rad_vals - kx_dirty)) ne 0 then begin message, 'kx_rad_vals for dirty and model cubes must match' endif @@ -232,13 +98,7 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ data_cube = temporary(dirty_cube) - temporary(model_cube) - if n_elements(freq_flags) gt 0 then begin - save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, $ - data_cube, freq_mask, uvf_git_hash - endif else begin - save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, $ - data_cube, uvf_git_hash - endelse + save, file = uvf_savefile[i], kx_rad_vals, ky_rad_vals, file_struct.frequencies, data_cube, uvf_git_hash undefine, data_cube, uvf_git_hash endif else begin ;; endif derived cube @@ -282,13 +142,6 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ endif if n_elements(wh_close) ne n_elements(pixel_nums) then arr = arr[wh_close, *] - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif pixels = pixel_nums[wh_close] @@ -316,7 +169,7 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ beam_git_hash = this_run_git_hash - save, file=file_struct.beam_savefile[i], avg_beam, pixels, nside, beam_git_hash + save, file=beam_savefile[i], avg_beam, pixels, nside, beam_git_hash endif @@ -443,13 +296,6 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ if n_elements(wh_close) ne n_elements(pixel_nums) then begin arr = arr[wh_close, *] endif - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif transform = discrete_ft_2D_fast(x_use, y_use, arr, kx_rad_vals, ky_rad_vals, $ timing = ft_time, fchunk = uvf_options.dft_fchunk, no_progress = uvf_options.no_dft_progress) @@ -457,13 +303,7 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ undefine, arr uvf_git_hash = this_run_git_hash - if n_elements(freq_flags) then begin - save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, $ - data_cube, freq_mask, uvf_git_hash - endif else begin - save, file = file_struct.uvf_savefile[i], kx_rad_vals, ky_rad_vals, $ - data_cube, uvf_git_hash - endelse + save, file = uvf_savefile[i], kx_rad_vals, ky_rad_vals, file_struct.frequencies, data_cube, uvf_git_hash undefine, data_cube, uvf_git_hash endif @@ -496,13 +336,6 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ if n_elements(wh_close) ne n_elements(pixel_nums) then begin arr = arr[wh_close, *] endif - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif transform = discrete_ft_2D_fast(x_use, y_use, arr, kx_rad_vals, ky_rad_vals, $ timing = ft_time, fchunk = uvf_options.dft_fchunk, no_progress = uvf_options.no_dft_progress) @@ -539,13 +372,6 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ if n_elements(wh_close) ne n_elements(pixel_nums) then begin arr = arr[wh_close, *] endif - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif transform = discrete_ft_2D_fast(x_use, y_use, arr, kx_rad_vals, ky_rad_vals, $ timing = ft_time, fchunk = uvf_options.dft_fchunk, no_progress = uvf_options.no_dft_progress) @@ -556,20 +382,15 @@ pro ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ uvf_wt_git_hash = this_run_git_hash - if n_elements(freq_flags) then begin - save, file = file_struct.uvf_weight_savefile[i], kx_rad_vals, $ - ky_rad_vals, weights_cube, variance_cube, freq_mask, uvf_wt_git_hash - endif else begin - save, file = file_struct.uvf_weight_savefile[i], kx_rad_vals, $ - ky_rad_vals, weights_cube, variance_cube, uvf_wt_git_hash - endelse + save, file = uvf_weight_savefile[i], kx_rad_vals, ky_rad_vals, file_struct.frequencies, $ + weights_cube, variance_cube, uvf_wt_git_hash undefine, new_pix_vec, weights_cube, variance_cube, uvf_wt_git_hash endif endelse endif else begin - kx_rad_vals = getvar_savefile(file_struct.uvf_savefile[0], 'kx_rad_vals') - ky_rad_vals = getvar_savefile(file_struct.uvf_savefile[0], 'ky_rad_vals') + kx_rad_vals = getvar_savefile(uvf_savefile[0], 'kx_rad_vals') + ky_rad_vals = getvar_savefile(uvf_savefile[0], 'ky_rad_vals') endelse endfor end diff --git a/ps_core/ps_kcube.pro b/ps_core/ps_kcube.pro index 38b9a53..5648cb1 100644 --- a/ps_core/ps_kcube.pro +++ b/ps_core/ps_kcube.pro @@ -1,6 +1,6 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ - uvf_input = uvf_input, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - input_units = input_units, save_slices = save_slices, $ + uvf_input = uvf_input, freq_options = freq_options, $ + input_units = input_units, save_slices = save_slices, save_sum_cube = save_sum_cube, $ refresh_options = refresh_options, uvf_options = uvf_options, $ ps_options = ps_options @@ -22,37 +22,103 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ git, repo_path = ps_repository_dir(), result=this_run_git_hash - frequencies = file_struct.frequencies + n_vis_freq = reform(file_struct.n_vis_freq) + + if healpix or not uvf_input then begin + ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ + no_var = no_var, refresh_options = refresh_options, uvf_options = uvf_options + endif - if n_elements(freq_flags) ne 0 then freq_mask = file_struct.freq_mask + if ( $ + tag_exist(freq_options, 'freq_flags') $ + or tag_exist(freq_options, 'freq_ch_range') $ + or tag_exist(freq_options, "freq_avg_bins") $ + or freq_options.freq_avg_factor gt 1 $ + ) then begin + ps_freq_select_avg, file_struct, n_vis_freq, refresh_options = refresh_options, $ + freq_options = freq_options, ps_options = ps_options + frequencies = freq_options.frequencies + full_n_vis_freq = n_vis_freq + if tag_exist(freq_options, 'freq_flags') then begin + full_n_vis_freq = full_n_vis_freq * transpose([[freq_options.freq_mask], [freq_options.freq_mask]]) + endif + n_vis_freq = freq_options.n_vis_freq + endif else begin + frequencies = file_struct.frequencies + endelse + + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask = freq_options.freq_mask + endif - n_freq_orig = n_elements(frequencies) ;Read-in number of frequencies - if n_elements(freq_ch_range) ne 0 then frequencies = frequencies[min(freq_ch_range):max(freq_ch_range)] n_freq = n_elements(frequencies) ;Selected number of frequencies + if tag_exist(freq_options, 'freq_flags') then begin + if tag_exist(freq_options, 'new_freq_mask') then begin + mask_use = freq_options.new_freq_mask + endif else begin + mask_use = freq_options.freq_mask + endelse + wh_unflagged = where(mask_use gt 0, count_unflagged) + if count_unflagged eq 0 then message, 'something went wrong with freq flagging' + n_freq_unflagged = total(mask_use) + endif else begin + n_freq_unflagged = n_freq + endelse z_mpc_mean = z_mpc(frequencies, hubble_param = hubble_param, f_delta = f_delta, $ - redshifts = redshifts, comov_dist_los = comov_dist_los, $ - z_mpc_delta = z_mpc_delta) + redshifts = redshifts, comov_dist_los = comov_dist_los, even_freq = even_freq, $ + z_mpc_delta = z_mpc_delta, z_mpc_length = z_mpc_length) kperp_lambda_conv = z_mpc_mean / (2.*!dpi) delay_delta = 1e9/(n_freq*f_delta*1e6) ;; equivilent delay bin size for kparallel delay_max = delay_delta * n_freq/2. ;; factor of 2 b/c of neg/positive delay_params = [delay_delta, delay_max] - z_mpc_length = max(comov_dist_los) - min(comov_dist_los) + z_mpc_delta - kz_mpc_range = (2.*!dpi) / (z_mpc_delta) - kz_mpc_delta = (2.*!dpi) / z_mpc_length + ;; check for uneven frequencies to determine how to calculate kz values + freq_diff = frequencies - shift(frequencies, 1) + freq_diff = freq_diff[1:*] + if max(abs(freq_diff-freq_diff[0])) gt 1e-12 and not tag_exist(ps_options, 'kz_use') then begin + message, "frequencies are unevenly spaced. To calculate kz values, " $ + + "kz_mpc_range and kz_mpc_delta must be sent in with kz_use. See dictionary " $ + + "and UW thesis by Star for reference." + endif + + if tag_exist(ps_options, 'kz_use') then begin + kz_mpc_range = double(ps_options.kz_use[0]) + kz_mpc_delta = double(ps_options.kz_use[1]) + n_kz = round(kz_mpc_range / kz_mpc_delta) + ; force frequency dft if using custom kz array + ps_options = create_ps_options(ps_options = ps_options, freq_dft = 1) + if n_kz lt 1 then begin + message, "n_kz is less than one. Something is wrong with kz_use input." + endif + endif else begin + kz_mpc_range = (2.*!dpi) / (z_mpc_delta) + kz_mpc_delta = (2.*!dpi) / z_mpc_length + n_kz = round(kz_mpc_range / kz_mpc_delta) + if n_kz ne n_freq then begin + message, "something has gone wrong with kz_mpc calculation, n_kz is " $ + + "not equal to n_freq. n_kz: " + number_formatter(n_kz) + ", n_freq:" $ + + number_formatter(n_freq) + endif + endelse + + ;; for standard EOR high band calculations (averaged evenly by four): + ;; kz_mpc_range = double(2.4628338337796376) + ;; kz_mpc_delta = double(0.012827259550935612) + ;; for hyperfine, kz_mpc_delta = double(0.00012827259550935612) ;; need to figure out where the zero bin is (for later computations) ;; The location of the zero bin also dictates what the most negative kz bin is (min_kz) ;; which we use to construct the kz_mpc_orig array ;; The calculation depends a bit on whether we have an odd or even number of frequencies. ;; The following calculation is taken directly from the IDL FFT documentation - int_arr = findgen((n_freq - 1)/2) + 1 - is_n_freq_even = (n_freq mod 2) eq 0 - fft_shift_val = -(n_freq/2 + 1) - if (is_n_freq_even) then begin - kz_integers = [0.0, int_arr, n_freq/2, -n_freq/2 + int_arr] + ;; The calculation was modified Jan 2024 to allow for n_kz not equal to n_freqs + int_arr = findgen((n_kz - 1)/2) + 1 + is_n_kz_even = (n_kz mod 2) eq 0 + fft_shift_val = -(n_kz/2 + 1) + if (is_n_kz_even) then begin + kz_integers = [0.0, int_arr, n_kz/2, -n_kz/2 + int_arr] endif else begin kz_integers = [0.0, int_arr, fft_shift_val + int_arr] endelse @@ -62,10 +128,7 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ endif min_kz = min(kz_integers)*kz_mpc_delta - kz_mpc_orig = findgen(round(kz_mpc_range / kz_mpc_delta)) * kz_mpc_delta + min_kz - - n_kz = n_elements(kz_mpc_orig) - if n_kz ne n_freq then message, 'something has gone wrong with kz_mpc calculation.' + kz_mpc_orig = dindgen(n_kz) * kz_mpc_delta + min_kz if input_units eq 'jansky' then begin ;; converting from Jy (in u,v,f) to mK*str (10^-26 * c^2 * 10^3/ (2*f^2*kb)) @@ -101,8 +164,8 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ if nfiles eq 2 then vis_sigma_ian = total(vis_sigma_ian, 1)/2. endif - if n_elements(freq_ch_range) ne 0 then begin - vis_sig_tag = number_formatter(384./n_freq_orig) + if tag_exist(freq_options, 'freq_ch_range') then begin + vis_sig_tag = number_formatter(384./n_elements(file_struct.frequencies)) endif else begin vis_sig_tag = number_formatter(384./n_freq) endelse @@ -111,22 +174,81 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ if file_valid(vis_sigma_file) then begin vis_sigma_adam = getvar_savefile(vis_sigma_file, 'vis_sigma') - if n_elements(freq_ch_range) ne 0 then begin - if n_elements(vis_sigma_adam) ne n_freq_orig then message, $ - 'vis_sig file has incorrect number of frequency channels' - vis_sigma_adam = vis_sigma_adam[min(freq_ch_range):max(freq_ch_range)] - endif else begin - if n_elements(vis_sigma_adam) ne n_freq then message, $ + if tag_exist(freq_options, 'freq_ch_range') then begin + if n_elements(vis_sigma_adam) ne n_elements(file_struct.frequencies) then begin + message, 'vis_sig file has incorrect number of frequency channels' + endif + vis_sigma_adam = vis_sigma_adam[min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif + + if freq_options.freq_avg_factor gt 1 then begin + new_vis_sigma_adam = reform(vis_sigma_adam, freq_options.freq_avg_factor, n_elements(frequencies)) + ;; average over the file axis for full_n_vis_freq because we use the same vis_sigma for even & odd + full_n_vis_freq_temp = total(full_n_vis_freq, 1)/2. + full_n_vis_freq_temp = reform(full_n_vis_freq_temp, freq_options.freq_avg_factor, n_elements(frequencies)) + new_vis_sigma_adam = sqrt( $ + total(new_vis_sigma_adam^2 * full_n_vis_freq_temp, 1) / total(full_n_vis_freq_temp, 1)) + vis_sigma_adam = new_vis_sigma_adam + endif + + if tag_exist(freq_options, 'freq_avg_bins') then begin + ;; average over the file axis for full_n_vis_freq because we use the same vis_sigma for even & odd + full_n_vis_freq_temp = total(full_n_vis_freq, 1)/2. + h = histogram(freq_options.freq_avg_bins, reverse_indices=ri) + new_nfreq = n_elements(h) + new_vis_sigma_adam = dblarr(new_nfreq) + for fi = 0, new_nfreq - 1 do begin + ; ge t indices for data going into each frequency bin + bin_inds = ri[ri[fi] : ri[fi + 1] - 1] + wh_unflag = where(full_n_vis_freq_temp[bin_inds] ne 0, count_unflagged) + if count_unflagged gt 0 then begin + new_vis_sigma_adam[fi] = sqrt(total(vis_sigma_adam[bin_inds]^2 * full_n_vis_freq_temp[bin_inds]) / total(full_n_vis_freq_temp[bin_inds])) + endif else begin + new_vis_sigma_adam[fi] = 0 + endelse + endfor + vis_sigma_adam = new_vis_sigma_adam + endif + + if n_elements(vis_sigma_adam) ne n_freq then message, $ 'vis_sig file has incorrect number of frequency channels' - endelse wh_nan = where(finite(vis_sigma_adam) eq 0, count_nan) if count_nan gt 0 then vis_sigma_adam[wh_nan] = 0 endif if n_elements(vis_sigma_ian) gt 0 then begin - if n_elements(freq_ch_range) ne 0 then begin - vis_sigma_ian = vis_sigma_ian[min(freq_ch_range):max(freq_ch_range)] + if tag_exist(freq_options, 'freq_ch_range') ne 0 then begin + vis_sigma_ian = vis_sigma_ian[min(freq_options.freq_ch_range):max(freq_options.freq_ch_range)] + endif + if freq_options.freq_avg_factor gt 1 then begin + new_vis_sigma_ian = reform(vis_sigma_ian, freq_options.freq_avg_factor, n_elements(frequencies)) + ;; average over the file axis for full_n_vis_freq because we use the same vis_sigma for even & odd + ;; (we already averaged over that axis earlier) + full_n_vis_freq_temp = total(full_n_vis_freq, 1)/2. + full_n_vis_freq_temp = reform(full_n_vis_freq_temp, freq_options.freq_avg_factor, n_elements(frequencies)) + new_vis_sigma_ian = sqrt( $ + total(new_vis_sigma_ian^2 * full_n_vis_freq_temp, 1) / total(full_n_vis_freq_temp, 1)) + vis_sigma_ian = new_vis_sigma_ian + endif + + if tag_exist(freq_options, 'freq_avg_bins') then begin + ;; average over the file axis for full_n_vis_freq because we use the same vis_sigma for even & odd + full_n_vis_freq_temp = total(full_n_vis_freq, 1)/2. + h = histogram(freq_options.freq_avg_bins, reverse_indices=ri) + new_nfreq = n_elements(h) + new_vis_sigma_ian = dblarr(new_nfreq) + for fi = 0, new_nfreq - 1 do begin + ; ge t indices for data going into each frequency bin + bin_inds = ri[ri[fi] : ri[fi + 1] - 1] + wh_unflag = where(full_n_vis_freq_temp[bin_inds] ne 0, count_unflagged) + if count_unflagged gt 0 then begin + new_vis_sigma_ian[fi] = sqrt(total(vis_sigma_ian[bin_inds]^2 * full_n_vis_freq_temp[bin_inds]) / total(full_n_vis_freq_temp[bin_inds])) + endif else begin + new_vis_sigma_ian[fi] = 0 + endelse + endfor + vis_sigma_ian = new_vis_sigma_ian endif vis_sigma = vis_sigma_ian vs_name = 'ian' @@ -150,18 +272,13 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ t_sys_meas = (eff_area * sqrt(df * tau) * vis_sigma * sqrt(2)) / $ ((2. * (1.38065e-23) * 1e26)) - n_vis = reform(file_struct.n_vis) - n_vis_freq = reform(file_struct.n_vis_freq) - if n_elements(freq_ch_range) ne 0 then begin - n_vis = total(n_vis_freq[*, min(freq_ch_range):max(freq_ch_range)], 2) - n_vis_freq = n_vis_freq[*, min(freq_ch_range):max(freq_ch_range)] - endif + if nfiles eq 2 then begin + n_vis = total(n_vis_freq, 2) + endif else begin + n_vis = total(n_vis_freq) + endelse if healpix or not uvf_input then begin - ps_image_to_uvf, file_struct, n_vis_freq, kx_rad_vals, ky_rad_vals, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, no_var = no_var, $ - refresh_options = refresh_options, uvf_options = uvf_options - n_kx = n_elements(kx_rad_vals) kx_rad_delta = kx_rad_vals[1] - kx_rad_vals[0] kx_mpc = temporary(kx_rad_vals) / z_mpc_mean @@ -253,57 +370,6 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ pix_area_rad = ((!dpi/180.)*file_struct.degpix)^2. pix_area_mpc = pix_area_rad * z_mpc_mean^2. - ;; get beam sorted out - if tag_exist(file_struct, 'beam_savefile') then begin - test_beam = file_valid(file_struct.beam_savefile) - if min(test_beam) eq 0 then test_beam = check_old_path(file_struct, 'beam_savefile') - - if min(test_beam) eq 0 or refresh_options.refresh_beam then begin - - for i=0, nfiles-1 do begin - arr = getvar_savefile(file_struct.beamfile[i], file_struct.beamvar) - void = getvar_savefile(file_struct.beamfile[i], names = beam_varnames) - wh_obs = where(stregex(strlowcase(beam_varnames), 'obs', /boolean), count_obs) - if count_obs gt 0 then obs_struct_name = beam_varnames[wh_obs[0]] - obs_beam = getvar_savefile(file_struct.beamfile[i], obs_struct_name) - nfvis_beam = obs_beam.nf_vis - undefine_fhd, obs_beam - - if n_elements(freq_ch_range) ne 0 then begin - nfvis_beam = nfvis_beam[min(freq_ch_range):max(freq_ch_range)] - arr = arr[*,*,min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - nfvis_beam = nfvis_beam*freq_mask - flag_arr = rebin(reform(freq_mask, 1, 1, n_freq), $ - size(arr,/dimension), /sample) - arr = arr * flag_arr - endif - - if max(arr) le 1.1 then begin - ;; beam is peak normalized to 1 - temp = arr * rebin(reform(nfvis_beam, 1, 1, n_elements(nfvis_beam)), $ - n_kx, n_ky, n_elements(nfvis_beam), /sample) - endif else if max(arr) le file_struct.n_obs[i] then begin - ;; beam is peak normalized to 1 for each obs, then summed over obs so peak is ~ n_obs - temp = (arr/file_struct.n_obs[i]) * rebin(reform(nfvis_beam, 1, 1, n_freq), $ - n_kx, n_ky, n_freq, /sample) - endif else begin - ;; beam is peak normalized to 1 then multiplied by n_vis_freq for each obs & summed - temp = arr - endelse - - avg_beam = total(temp, 3) / total(nfvis_beam) - - git, repo_path = ps_repository_dir(), result=beam_git_hash - - save, file=file_struct.beam_savefile[i], avg_beam, beam_git_hash - endfor - - endif - - endif - if tag_exist(uvf_options, 'uv_avg') then begin nkx_new = floor(n_kx / uvf_options.uv_avg) temp = dcomplex(dblarr(nkx_new, n_ky, n_freq)) @@ -475,19 +541,6 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ file_struct.variancevar, n_freq, file_struct.pol_index) endif - if n_elements(freq_ch_range) ne 0 then begin - variance_cube1 = variance_cube1[*, *, min(freq_ch_range):max(freq_ch_range)] - if nfiles eq 2 then begin - variance_cube2 = variance_cube2[*, *, min(freq_ch_range):max(freq_ch_range)] - endif - endif - if n_elements(freq_flags) ne 0 then begin - flag_arr = rebin(reform(freq_mask, 1, 1, n_freq), $ - size(variance_cube1,/dimension), /sample) - variance_cube1 = variance_cube1 * flag_arr - if nfiles eq 2 then variance_cube2 = variance_cube2 * flag_arr - endif - if max(abs(variance_cube1)) eq 0 then message, 'variance cube is entirely zero.' if nfiles eq 2 then if max(abs(variance_cube2)) eq 0 then begin message, 'variance cube is entirely zero.' @@ -583,6 +636,11 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ endif endif else beam_git_hashes = strarr(nfiles) + ;; TODO: think about whether to use n_freq or n_freq_unflagged in this section. + ;; first worry is that the volume calc is wrong if we include flagged freqs. + ;; worry with using n_freq_unflagged is that it's not clear if it's correct and + ;; that if the flags came through FHD instead of eppsilon they + ;; wouldn't be accounted for this way, so would get a different answer. if nfiles eq 2 then begin window_int_beam = [total(beam1), total(beam2)]*pix_area_mpc*(z_mpc_delta * n_freq) endif else begin @@ -603,36 +661,22 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ endif if tag_exist(file_struct, 'beam_int') then begin - beam_int = file_struct.beam_int - n_vis_freq = file_struct.n_vis_freq + if ( $ + tag_exist(freq_options, 'freq_flags') $ + or tag_exist(freq_options, 'freq_ch_range') $ + or freq_options.freq_avg_factor gt 1 $ + or tag_exist(freq_options, 'freq_avg_bins') $ + ) then begin + ave_beam_int = freq_options.beam_int + endif else begin + beam_int = file_struct.beam_int + n_vis_freq = file_struct.n_vis_freq - if n_elements(freq_ch_range) ne 0 then begin - if nfiles eq 2 then begin - beam_int = beam_int[*,min(freq_ch_range):max(freq_ch_range)] - n_vis_freq = n_vis_freq[*,min(freq_ch_range):max(freq_ch_range)] - endif else begin - beam_int = beam_int[min(freq_ch_range):max(freq_ch_range)] - n_vis_freq = n_vis_freq[min(freq_ch_range):max(freq_ch_range)] - endelse - endif - if n_elements(freq_flags) ne 0 then begin if nfiles eq 2 then begin - flag_arr = rebin(reform(freq_mask, 1, n_freq), $ - size(beam_int,/dimension), /sample) - beam_int = beam_int*flag_arr - n_vis_freq = n_vis_freq*flag_arr + ave_beam_int = total(beam_int * n_vis_freq, 2) / total(n_vis_freq, 2) endif else begin - beam_int = beam_int*freq_mask - n_vis_freq = n_vis_freq*freq_mask + ave_beam_int = total(beam_int * n_vis_freq) / total(n_vis_freq) endelse - endif - - if nfiles eq 2 then begin - ave_beam_int = total(beam_int * n_vis_freq, 2) / total(n_vis_freq, 2) - endif else begin - ave_beam_numerator = total(beam_int * n_vis_freq) - ave_beam_denominator = total(n_vis_freq) - ave_beam_int = total(beam_int * n_vis_freq) / total(n_vis_freq) endelse wh_not_finite = where(~finite(ave_beam_int), count_not_finite) @@ -640,6 +684,11 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ message, 'some frequency averaged beam integrals are not finite' endif ;; convert rad -> Mpc^2, multiply by depth in Mpc + ;; TODO: think about whether to use n_freq or n_freq_unflagged in this section. + ;; first worry is that the volume calc is wrong if we include flagged freqs. + ;; worry with using n_freq_unflagged is that it's not clear if it's correct and + ;; that if the flags came through FHD instead of eppsilon they + ;; wouldn't be accounted for this way, so would get a different answer. window_int_beam_obs = ave_beam_int * z_mpc_mean^2. * (z_mpc_delta * n_freq) endif @@ -874,9 +923,13 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ if total(abs(sigma2_cube1)) le 0 then message, 'sigma2 cube is all zero' endelse - - wt_meas_ave = total(abs(weights_cube1), 3)/n_freq - wt_meas_min = min(abs(weights_cube1), dimension=3) + ;; need to use n_freq_unflagged here! + wt_meas_ave = total(abs(weights_cube1), 3)/n_freq_unflagged + if tag_exist(freq_options, 'freq_flags') then begin + wt_meas_min = min(abs(weights_cube1[*, *, wh_unflagged]), dimension=3) + endif else begin + wt_meas_min = min(abs(weights_cube1), dimension=3) + endelse ;; divide data by weights data_cube1 = data_cube1 / weights_cube1 @@ -1012,11 +1065,11 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ endif ;; make simulated noise cubes - sim_noise1 = randomn(seed, n_kx, n_ky, n_kz) * sqrt(sigma2_cube1) + $ - dcomplex(0,1) * randomn(seed, n_kx, n_ky, n_kz) * sqrt(sigma2_cube1) + sim_noise1 = randomn(seed, n_kx, n_ky, n_freq) * sqrt(sigma2_cube1) + $ + dcomplex(0,1) * randomn(seed, n_kx, n_ky, n_freq) * sqrt(sigma2_cube1) if nfiles eq 2 then begin - sim_noise2 = randomn(seed, n_kx, n_ky, n_kz) * sqrt(sigma2_cube2) + $ - dcomplex(0,1) * randomn(seed, n_kx, n_ky, n_kz) * sqrt(sigma2_cube2) + sim_noise2 = randomn(seed, n_kx, n_ky, n_freq) * sqrt(sigma2_cube2) + $ + dcomplex(0,1) * randomn(seed, n_kx, n_ky, n_freq) * sqrt(sigma2_cube2) endif if nfiles eq 2 then begin @@ -1091,7 +1144,7 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ sim_noise_sum = temporary(sim_noise1) endelse - mask = intarr(n_kx, n_ky, n_kz) + 1 + mask = intarr(n_kx, n_ky, n_freq) + 1 wh_sig0 = where(sum_sigma2 eq 0, count_sig0) if count_sig0 gt 0 then mask[wh_sig0] = 0 n_freq_contrib = total(mask, 3) @@ -1115,7 +1168,7 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ endif if ps_options.save_sum_cube then begin - if n_elements(freq_flags) gt 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = file_struct.uvf_sum_savefile, data_sum, kx_mpc, ky_mpc, frequencies, kperp_lambda_conv, $ delay_params, hubble_param, freq_mask endif else begin @@ -1177,7 +1230,9 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ endif if ps_options.inverse_covar_weight then begin - + if tag_exist(ps_options, 'kz_use') then begin + print, 'caution! Custom kz arrays have not been tested with inverse covariance weighting.' + endif max_eta_val = max(sum_sigma2) eta_var = shift([max_eta_val, fltarr(n_freq-1)], fft_shift_val) covar_eta_fg = diag_matrix(eta_var) @@ -1238,15 +1293,17 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ ;; this should be zero by construction, but add a test just in case if kz_mpc_orig[where(n_val eq 0)] ne 0 then begin - message, 'something has gone terribly wrong with calculating the kz_mpc_orig values.' + message, 'center kz value is not zero. Something has gone terribly wrong ' $ + + 'with calculating the kz_mpc_orig values.' endif kz_mpc = kz_mpc_orig[where(n_val ge 0)] - if is_n_freq_even then begin - ;; if n_freq is even, there's an extra positive mode (because there's always a zero mode.) + if is_n_kz_even then begin + ;; if n_kz is even, there's an extra positive mode (because there's always a zero mode.) if max(abs(n_val(wh_pos))) le max(abs(n_val(wh_neg))) then begin - message, 'something has gone terribly wrong with calculating the kz_integers.' + message, 'absolute values of max and min kz_integers do not match. Something has ' $ + + 'gone terribly wrong with calculating the kz_integers.' endif ;; remove unmatched max k positive mode: trim_max_pos = 1 @@ -1547,16 +1604,16 @@ pro ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ git_hashes = {uvf:uvf_git_hashes, uvf_wt:uvf_wt_git_hashes, beam:beam_git_hashes, $ kcube:kcube_git_hash} - if n_elements(freq_flags) gt 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = file_struct.kcube_savefile, data_sum_1, data_sum_2, data_diff_1, $ - data_diff_2, sigma2_1, sigma2_2, sim_noise_sum_1, sim_noise_sum_2, $ + data_diff_2, sigma2_1, sigma2_2, sim_noise_sum_1, sim_noise_sum_2, frequencies, $ sim_noise_diff_1, sim_noise_diff_2, kx_mpc, ky_mpc, kz_mpc, kperp_lambda_conv, $ delay_params, hubble_param, n_freq_contrib, freq_mask, vs_name, vs_mean, $ t_sys_meas, window_int, git_hashes, wt_meas_ave, wt_meas_min, ave_weights, $ wt_ave_power_freq, ave_power_freq, wt_ave_power_uvf, ave_power_uvf endif else begin save, file = file_struct.kcube_savefile, data_sum_1, data_sum_2, data_diff_1, $ - data_diff_2, sigma2_1, sigma2_2, sim_noise_sum_1, sim_noise_sum_2, $ + data_diff_2, sigma2_1, sigma2_2, sim_noise_sum_1, sim_noise_sum_2, frequencies, $ sim_noise_diff_1, sim_noise_diff_2, kx_mpc, ky_mpc, kz_mpc, kperp_lambda_conv, $ delay_params, hubble_param, n_freq_contrib, vs_name, vs_mean, t_sys_meas, $ window_int, git_hashes, wt_meas_ave, wt_meas_min, ave_weights, $ diff --git a/ps_core/ps_main_plots.pro b/ps_core/ps_main_plots.pro index feb6a12..09e0841 100644 --- a/ps_core/ps_main_plots.pro +++ b/ps_core/ps_main_plots.pro @@ -1,6 +1,6 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ - type_inc = type_inc, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, uvf_input = uvf_input, no_evenodd = no_evenodd, $ + type_inc = type_inc, freq_options = freq_options, $ + uvf_input = uvf_input, no_evenodd = no_evenodd, $ rts = rts, norm_rts_with_fhd = norm_rts_with_fhd, casa = casa, sim = sim, $ fix_sim_input = fix_sim_input, save_path = save_path, $ savefilebase = savefilebase, cube_power_info = cube_power_info, $ @@ -17,21 +17,18 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if keyword_set(rts) then begin file_struct_arr = rts_file_setup(datafile, savefilebase = savefilebase, $ save_path = save_path, use_fhd_norm = norm_rts_with_fhd, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ refresh_info = refresh_options.refresh_info, uvf_options = uvf_options, $ - ps_options = ps_options) + freq_options = freq_options, ps_options = ps_options) endif else if keyword_set(casa) then begin file_struct_arr = casa_file_setup(datafile, savefilebase = savefilebase, $ save_path = save_path, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ refresh_info = refresh_options.refresh_info, uvf_options = uvf_options, $ - ps_options = ps_options) + freq_options = freq_options, ps_options = ps_options) endif else begin file_struct_arr = fhd_file_setup(datafile, beamfile = beamfiles, sim = sim, $ uvf_input = uvf_input, savefilebase = savefilebase, save_path = save_path, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ refresh_info = refresh_options.refresh_info, uvf_options = uvf_options, $ - ps_options = ps_options, pol_inc=pol_inc) + freq_options = freq_options, ps_options = ps_options, pol_inc=pol_inc) endelse time1 = systime(1) print, 'file setup time: ' + number_formatter(time1-time0) @@ -136,17 +133,13 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if tag_exist(file_struct_arr, 'uvf_savefile') then uvf_input = 0 else uvf_input = 1 - n_freq = n_elements(file_struct_arr[0].frequencies) - if n_elements(freq_ch_range) ne 0 then begin - if max(freq_ch_range) gt n_freq-1 then message, 'invalid freq_ch_range' - n_freq = freq_ch_range[1]-freq_ch_range[0]+1 - endif - if healpix and tag_exist(uvf_options, 'dft_fchunk') then begin - if uvf_options.dft_fchunk gt n_freq then begin + ;; Use the original frequencies here because we DFT all frequencies + n_orig_freqs = n_elements(file_struct_arr[0].frequencies) + if uvf_options.dft_fchunk gt n_orig_freqs then begin print, 'dft_fchunk is larger than the number of frequency slices, setting ' + $ - 'it to the number of slices -- ' + number_formatter(n_freq) - uvf_options.dft_fchunk = n_freq + 'it to the number of slices -- ' + number_formatter(n_orig_freqs) + uvf_options.dft_fchunk = n_orig_freqs endif endif @@ -250,8 +243,8 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ power_file = file_struct_arr[i].power_savefile test_power_3d = file_valid(power_file) - if test_power_3d eq 1 and n_elements(freq_flags) ne 0 then begin - freq_mask = file_struct_arr[i].freq_mask + if test_power_3d eq 1 and tag_exist(freq_options, 'freq_flags') then begin + freq_mask = freq_options.freq_mask old_freq_mask = getvar_savefile(power_file, 'freq_mask') if total(abs(old_freq_mask - freq_mask)) ne 0 then test_power_3d = 0 endif @@ -270,7 +263,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if test_power_3d eq 0 or refresh_options.refresh_ps then begin ps_power, file_struct_arr[i], sim = sim, fix_sim_input = fix_sim_input, $ - uvf_input = uvf_input, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ + uvf_input = uvf_input, freq_options = freq_options, $ input_units = input_units, save_slices = save_slices, $ refresh_options = refresh_options, uvf_options = uvf_options, $ ps_options = ps_options @@ -281,12 +274,8 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ ;; setup for saving plots power_tag = file_struct_arr[0].power_tag - if tag_exist(file_struct_arr[0], 'uvf_tag') then begin - uvf_tag = file_struct_arr[0].uvf_tag - endif else begin - uvf_tag = '' - endelse - + uvf_tag = file_struct_arr[0].uvf_tag + freq_tag = file_struct_arr[0].freq_tag ;; need general_filebase for 1D and weight plotfiles, make sure it doesn't have a full path general_filebase = file_struct_arr[0].general_filebase @@ -321,16 +310,16 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ plotfile_base_wt = general_filebase + '_' + $ file_struct_arr[uniq(weight_ind, sort(weight_ind))].pol + power_tag endif else begin - plotfile_base = plot_options.plot_filebase + uvf_tag + $ + plotfile_base = plot_options.plot_filebase + uvf_tag + freq_tag + $ file_struct_arr.file_label + power_tag - plotfile_base_wt = plot_options.plot_filebase + uvf_tag + $ + plotfile_base_wt = plot_options.plot_filebase + uvf_tag + freq_tag + $ '_' + file_struct_arr[uniq(weight_ind, sort(weight_ind))].pol + power_tag endelse endif else begin if not tag_exist(plot_options, 'plot_filebase') then begin plotfile_base = general_filebase + power_tag endif else begin - plotfile_base = plot_options.plot_filebase + uvf_tag + power_tag + plotfile_base = plot_options.plot_filebase + uvf_tag + freq_tag + power_tag endelse plotfile_base_wt = plotfile_base endelse @@ -338,10 +327,28 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ endif + if ( $ + tag_exist(freq_options, 'freq_flags') $ + or tag_exist(freq_options, 'freq_ch_range') $ + or freq_options.freq_avg_factor gt 1 $ + ) then begin + ps_freq_select_avg, file_struct_arr[0], reform(file_struct_arr[0].n_vis_freq), $ + refresh_options = refresh_options, freq_options = freq_options, $ + ps_options = ps_options, /metadata_only + frequencies = freq_options.frequencies + n_vis_freq = freq_options.n_vis_freq + n_vis = total(n_vis_freq, 2) + endif else begin + frequencies = file_struct_arr[0].frequencies + n_vis_freq = file_struct_arr[0].n_vis_freq + n_vis = file_struct_arr[0].n_vis[0] + endelse + n_freq = n_elements(frequencies) + ;; do this here because power slice plots need it. if plot_2d_options.plot_wedge_line or tag_exist(binning_1d_options, 'wedge_angles') then begin z0_freq = 1420.40 ;; MHz - redshifts = z0_freq/file_struct_arr[0].frequencies - 1 + redshifts = z0_freq/frequencies - 1 mean_redshift = mean(redshifts) cosmology_measures, mean_redshift, wedge_factor = wedge_factor, hubble_param = hubble_param @@ -731,13 +738,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if tag_exist(binning_1d_options, 'coarse_harm_width') then begin harm_freq = 1.28 - if n_elements(freq_ch_range) gt 0 then begin - freqs_use = file_struct_arr[0].frequencies[freq_ch_range[0]:freq_ch_range[1]] - endif else begin - freqs_use = file_struct_arr[0].frequencies - endelse - - bandwidth = max(freqs_use) - min(freqs_use) + freqs_use[1] - freqs_use[0] + bandwidth = max(frequencies) - min(frequencies) + frequencies[1] - frequencies[0] coarse_harm0 = round(bandwidth / harm_freq) binning_1d_options = create_binning_1d_options(binning_1d_options = binning_1d_options, $ coarse_harm0 = coarse_harm0) @@ -858,7 +859,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ for i=0, n_elements(wedge_1dbin_names)-1 do begin for j=0, n_elements(kperp_density_names)-1 do begin plotfile_binning_hist[*,j,i] = bin_hist_plot_path + plot_options.plot_filebase + $ - uvf_tag + file_struct_arr.file_label + power_tag + kperp_density_names[j] + $ + uvf_tag + freq_tag + file_struct_arr.file_label + power_tag + kperp_density_names[j] + $ wedge_1dbin_names[i] + binning_tags.fadd_1dbin endfor endfor @@ -945,7 +946,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if plot_types.plot_k0_power then begin if test_k0 gt 0 and tag_exist(binning_2d_options, 'kperp_bin') then begin kperp_bin_file = fltarr(n_elements(savefile_k0_use)) - for j=0, n_elements(savefile_k0_use) do begin + for j=0, n_elements(savefile_k0_use)-1 do begin kperp_bin_file[j] = getvar_savefile(savefile_k0_use[j], 'k_bin') endfor if max(abs(binning_2d_options.kperp_bin - kperp_bin_file)) gt 0. then test_k0=0 @@ -953,7 +954,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if test_k0_masked gt 0 and tag_exist(binning_2d_options, 'kperp_bin') then begin kperp_bin_file = fltarr(n_elements(savefile_k0_masked_use)) - for j=0, n_elements(savefile_k0_masked_use) do begin + for j=0, n_elements(savefile_k0_masked_use)-1 do begin kperp_bin_file[j] = getvar_savefile(savefile_k0_masked_use[j], 'k_bin') endfor if max(abs(binning_2d_options.kperp_bin - kperp_bin_file)) gt 0. then test_k0_masked=0 @@ -974,27 +975,27 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if max(abs(binning_1d_options.k_bin - k_bin_file)) gt 0. then test_bin=0 endif - if test_2d gt 0 and n_elements(freq_flags) ne 0 then begin + if test_2d gt 0 and tag_exist(freq_options, 'freq_flags') then begin for j=0, n_elements(savefile_2d_use)-1 do begin old_freq_mask = getvar_savefile(savefile_2d_use[j], 'freq_mask') - if total(abs(old_freq_mask - file_struct_arr[i].freq_mask)) ne 0 then test_2d = 0 + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then test_2d = 0 endfor for j=0, n_elements(savefiles_bin_use)-1 do begin old_freq_mask = getvar_savefile(savefiles_bin_use[j], 'freq_mask') - if total(abs(old_freq_mask - file_struct_arr[i].freq_mask)) ne 0 then test_bin = 0 + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then test_bin = 0 endfor for j=0, n_elements(savefiles_mask_use)-1 do begin old_freq_mask = getvar_savefile(savefiles_mask_use[j], 'freq_mask') - if total(abs(old_freq_mask - file_struct_arr[i].freq_mask)) ne 0 then test_mask = 0 + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then test_mask = 0 endfor endif - if test_1d gt 0 and n_elements(freq_flags) ne 0 then begin + if test_1d gt 0 and tag_exist(freq_options, 'freq_flags') then begin for j=0, n_elements(savefile_1d_use)-1 do begin old_freq_mask = getvar_savefile(savefile_1d_use[j], 'freq_mask') - if total(abs(old_freq_mask - file_struct_arr[i].freq_mask)) ne 0 then test_1d = 0 + if total(abs(old_freq_mask - freq_options.freq_mask)) ne 0 then test_1d = 0 endfor endif @@ -1148,7 +1149,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ plotfile_binning_hist_use = reform(plotfile_binning_hist[i,*,*], $ n_elements(kperp_density_names), n_elements(wedge_1dbin_names)) endif - ps_binning, file_struct_arr[i], sim = sim, freq_flags = freq_flags, $ + ps_binning, file_struct_arr[i], sim = sim, freq_options = freq_options, $ ps_options = ps_options, binning_2d_options = binning_2d_options, $ binning_1d_options = binning_1d_options, plot_options = plot_options, $ plot_types = plot_types, $ @@ -1226,7 +1227,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if not tag_exist(plot_options, 'plot_filebase') then begin plotfile_1d_base = general_filebase endif else begin - plotfile_1d_base = plot_options.plot_filebase + uvf_tag + plotfile_1d_base = plot_options.plot_filebase + uvf_tag + freq_tag endelse plotfile_path_1d = plotfile_path + file_struct_arr[0].subfolders.bin_1d begin_1d = plotfile_path_1d + plotfile_1d_base + power_tag + wedge_1dbin_names + binning_tags.fadd_1dbin @@ -1331,7 +1332,6 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ endif - ave_power_vals = fltarr(n_cubes) wt_ave_power_vals = fltarr(n_cubes) ave_weights_vals = fltarr(n_cubes) @@ -1367,13 +1367,8 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ if max(strmatch(varnames, 'uv_area', /fold_case)) gt 0 then begin uv_area[k] = getvar_savefile(savefiles_1d[k,0,0], 'uv_area') - nbsl_lambda2[k] = file_struct_arr[0].n_vis[0]/uv_area[k] - if n_elements(freq_ch_range) gt 0 then begin - nbsl_lambda2_freq[k,*] = file_struct_arr[0].n_vis_freq[0, $ - freq_ch_range[0]:freq_ch_range[1]]/uv_area[k] - endif else begin - nbsl_lambda2_freq[k,*] = file_struct_arr[0].n_vis_freq[0,*]/uv_area[k] - endelse + nbsl_lambda2[k] = n_vis[0]/uv_area[k] + nbsl_lambda2_freq[k,*] = n_vis_freq[0,*]/uv_area[k] endif else begin uv_area[k] = -1 nbsl_lambda2[k] = -1 diff --git a/ps_core/ps_power.pro b/ps_core/ps_power.pro index 39e4438..a549dbc 100644 --- a/ps_core/ps_power.pro +++ b/ps_core/ps_power.pro @@ -1,12 +1,14 @@ pro ps_power, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ - uvf_input = uvf_input, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ + uvf_input = uvf_input, freq_options = freq_options, $ input_units = input_units, save_slices = save_slices, $ refresh_options = refresh_options, uvf_options = uvf_options, $ ps_options = ps_options nfiles = n_elements(file_struct.datafile) - if n_elements(freq_flags) ne 0 then freq_mask = file_struct.freq_mask + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask = freq_options.freq_mask + endif test_kcube = file_valid(file_struct.kcube_savefile) if not test_kcube then test_kcube = check_old_path(file_struct, 'kcube_savefile') @@ -18,7 +20,7 @@ pro ps_power, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ if test_kcube eq 0 or refresh_options.refresh_kcube then begin ps_kcube, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ - uvf_input = uvf_input, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ + uvf_input = uvf_input, freq_options = freq_options, $ input_units = input_units, save_slices = save_slices, $ refresh_options = refresh_options, uvf_options = uvf_options, $ ps_options = ps_options @@ -258,7 +260,7 @@ pro ps_power, file_struct, sim = sim, fix_sim_input = fix_sim_input, $ kcube:'', ps:ps_git_hash} endelse - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin save, file = file_struct.power_savefile, power_3d, noise_3d, sim_noise_3d, $ sim_noise_diff_3d, noise_expval_3d, weights_3d, diff_weights_3d, $ kx_mpc, ky_mpc, kz_mpc, kperp_lambda_conv, delay_params, hubble_param, $ diff --git a/ps_core/single_cube_dft.pro b/ps_core/single_cube_dft.pro index c566822..52ecb34 100644 --- a/ps_core/single_cube_dft.pro +++ b/ps_core/single_cube_dft.pro @@ -2,7 +2,7 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ exact_obsnames = exact_obsnames, ps_foldername = ps_foldername, $ rts = rts, dirty_folder = dirty_folder, save_path = save_path, $ savefilebase = savefilebase, refresh_info = refresh_info, refresh_dft = refresh_dft, $ - freq_flags = freq_flags, freq_ch_range = freq_ch_range, cube_type = cube_type, $ + cube_type = cube_type, $ pol = pol, evenodd = evenodd, loc_name=loc_name, $ delta_uv_lambda = delta_uv_lambda, max_uv_lambda = max_uv_lambda, $ full_image = full_image @@ -84,20 +84,20 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ ;; it doesn't matter what's in ps_options because we don't use any of the filenames ;; that depend on this ps_options = create_ps_options() + freq_options = create_freq_options() if keyword_set(rts) then begin file_struct_arr = rts_file_setup(datafile, savefilebase = savefilebase, $ save_path = save_path,use_fhd_norm = norm_rts_with_fhd, refresh_info = refresh_info, $ - uvf_options = uvf_options, ps_options = ps_options) + uvf_options = uvf_options, ps_options = ps_options, freq_options=freq_options) endif else if keyword_set(casa) then begin file_struct_arr = casa_file_setup(datafile, savefilebase = savefilebase, $ save_path = save_path, refresh_info = refresh_info, $ - uvf_options = uvf_options, ps_options = ps_options) + uvf_options = uvf_options, ps_options = ps_options, freq_options=freq_options) endif else begin file_struct_arr = fhd_file_setup(datafile, savefilebase = savefilebase, $ - save_path = save_path, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, sim = sim, refresh_info = refresh_info, $ - uvf_options = uvf_options, ps_options = ps_options) + save_path = save_path, sim = sim, refresh_info = refresh_info, $ + uvf_options = uvf_options, ps_options = ps_options, freq_options=freq_options) endelse type_pol_str_list = file_struct_arr.type_pol_str @@ -155,15 +155,6 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ endelse test_uvf = file_test(uvf_savefile) * (1 - file_test(uvf_savefile, /zero_length)) - if test_uvf eq 1 and n_elements(freq_flags) ne 0 then begin - old_freq_mask = getvar_savefile(file_struct.uvf_savefile[i], 'freq_mask') - if n_elements(freq_ch_range) ne 0 then begin - if total(abs(old_freq_mask[min(freq_ch_range):max(freq_ch_range)] - freq_mask)) ne 0 $ - then test_uvf = 0 - endif else begin - if total(abs(old_freq_mask - freq_mask)) ne 0 then test_uvf = 0 - endelse - endif if test_uvf eq 0 or keyword_set(refresh_dft) then begin @@ -172,12 +163,6 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ frequencies = file_struct.frequencies - if n_elements(freq_flags) ne 0 then freq_mask = file_struct.freq_mask - - if n_elements(freq_ch_range) ne 0 then begin - n_freq_orig = n_elements(frequencies) - frequencies = frequencies[min(freq_ch_range):max(freq_ch_range)] - endif n_freq = n_elements(frequencies) pixel_nums = getvar_savefile(file_struct.pixelfile[0], file_struct.pixelvar[0]) @@ -219,13 +204,6 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ if n_elements(wh_close) ne n_elements(pixel_nums) then begin arr = arr[wh_close, *] endif - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif transform = discrete_ft_2D_fast(x_use, y_use, arr, kx_rad_vals, ky_rad_vals, $ timing = ft_time, fchunk = dft_fchunk, no_progress = no_dft_progress) @@ -255,13 +233,6 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ if n_elements(wh_close) ne n_elements(pixel_nums) then begin arr = arr[wh_close, *] endif - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif transform = discrete_ft_2D_fast(x_use, y_use, arr, kx_rad_vals, ky_rad_vals, $ timing = ft_time, fchunk = dft_fchunk, no_progress = no_dft_progress) @@ -271,13 +242,9 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ undefine, arr endif - if n_elements(freq_flags) then begin - save, file = uvf_savefile, kx_rad_vals, ky_rad_vals, weights_cube, $ - variance_cube, freq_mask, uvf_wt_git_hash - endif else begin - save, file = uvf_savefile, kx_rad_vals, ky_rad_vals, weights_cube, $ - variance_cube, uvf_wt_git_hash - endelse + save, file = uvf_savefile, kx_rad_vals, ky_rad_vals, weights_cube, $ + variance_cube, uvf_wt_git_hash + undefine, new_pix_vec, weights_cube, variance_cube @@ -304,13 +271,6 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ if n_elements(wh_close) ne n_elements(pixel_nums) then begin arr = arr[wh_close, *] endif - if n_elements(freq_ch_range) ne 0 then begin - arr = arr[*, min(freq_ch_range):max(freq_ch_range)] - endif - if n_elements(freq_flags) ne 0 then begin - arr = arr * rebin(reform(freq_mask, 1, n_freq), $ - size(arr, /dimension), /sample) - endif transform = discrete_ft_2D_fast(x_use, y_use, arr, kx_rad_vals, ky_rad_vals, $ timing = ft_time, fchunk = dft_fchunk, no_progress = no_dft_progress) @@ -318,11 +278,8 @@ pro single_cube_dft, folder_name_in, obs_name, data_subdirs=data_subdirs, $ data_cube = temporary(transform) undefine, arr - if n_elements(freq_flags) then begin - save, file = uvf_savefile, kx_rad_vals, ky_rad_vals, data_cube, freq_mask, uvf_git_hash - endif else begin - save, file = uvf_savefile, kx_rad_vals, ky_rad_vals, data_cube, uvf_git_hash - endelse + save, file = uvf_savefile, kx_rad_vals, ky_rad_vals, data_cube, uvf_git_hash + undefine, data_cube endelse endif diff --git a/ps_plotting/ps_comp1d_plots.pro b/ps_plotting/ps_comp1d_plots.pro index 620bd6a..3e45be5 100644 --- a/ps_plotting/ps_comp1d_plots.pro +++ b/ps_plotting/ps_comp1d_plots.pro @@ -1,8 +1,8 @@ pro ps_comp1d_plots, folder_names, obs_info, ps_foldernames = ps_foldernames, $ - cube_types, pols, uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, $ + cube_types, pols, uvf_options = uvf_options, freq_options = freq_options, $ ps_options = ps_options, plot_options = plot_options, $ binning_1d_options = binning_1d_options, names = names, $ - all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + all_type_pol = all_type_pol, $ plot_filebase = plot_filebase, save_path = save_path, savefilebase = savefilebase, $ quiet = quiet, window_num = window_num @@ -10,10 +10,10 @@ pro ps_comp1d_plots, folder_names, obs_info, ps_foldernames = ps_foldernames, $ compare_plot_prep, folder_names, obs_info, cube_types, pols, 'comp_1d', compare_files, $ ps_foldernames = ps_foldernames, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, ps_options = ps_options, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ binning_1d_options = binning_1d_options, $ - freq_ch_range = freq_ch_range, plot_filebase = plot_filebase, $ + plot_filebase = plot_filebase, $ save_path = save_path, savefilebase = savefilebase, $ full_compare = all_type_pol @@ -53,13 +53,22 @@ pro ps_comp1d_plots, folder_names, obs_info, ps_foldernames = ps_foldernames, $ if n_elements(names) ne 2 then message, 'names must have two elements' names_use = names endif else begin - names_use = strsplit(obs_info.diff_note, '-',/extract) + if n_elements(folder_names) eq 2 or n_elements(ps_foldernames) eq 2 then begin + names_use = strsplit(obs_info.diff_note, '-', /extract) + endif else begin + names_use = strsplit(plot_options.note, 'and', /regex, /extract) + endelse + endelse + if i lt n_elements(compare_files.wedge_amp) then begin + wedge_name_use = binning_1d_options.wedge_names[i] + endif else begin + wedge_name_use = "no wedge cut" endelse kpower_1d_plots, input_files, window_num=window_num, $ start_multi_params = start_multi_params, multi_pos = pos_use, $ names=names_use, plot_options = plot_options, note = note_use, $ - plotfile = plotfiles_use, title = compare_files.titles[j] + plotfile = plotfiles_use, title = compare_files.titles[j] + ' ' + wedge_name_use if j eq 0 and compare_files.n_cubes gt 1 then begin positions = pos_use diff --git a/ps_plotting/ps_difference_plots.pro b/ps_plotting/ps_difference_plots.pro index ff58184..2d52dcd 100644 --- a/ps_plotting/ps_difference_plots.pro +++ b/ps_plotting/ps_difference_plots.pro @@ -1,10 +1,9 @@ pro ps_difference_plots, folder_names, obs_info, ps_foldernames = ps_foldernames, $ cube_types, pols, all_type_pol = all_type_pol, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, ps_options = ps_options, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ binning_2d_options = binning_2d_options, binning_1d_options = binning_1d_options, $ - refresh_diff = refresh_diff, freq_ch_range = freq_ch_range, $ - plot_slices = plot_slices, slice_type = slice_type, $ + refresh_diff = refresh_diff, plot_slices = plot_slices, slice_type = slice_type, $ plot_filebase = plot_filebase, save_path = save_path, savefilebase = savefilebase, $ data_min_abs = data_min_abs, plot_1d = plot_1d, axis_type_1d=axis_type_1d, $ invert_colorbar = invert_colorbar, $ @@ -12,12 +11,11 @@ pro ps_difference_plots, folder_names, obs_info, ps_foldernames = ps_foldernames compare_plot_prep, folder_names, obs_info, cube_types, pols, 'diff', compare_files, $ ps_foldernames = ps_foldernames, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, ps_options = ps_options, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ binning_2d_options = binning_2d_options, binning_1d_options = binning_1d_options, $ plot_slices = plot_slices, slice_type = slice_type, $ - freq_ch_range = freq_ch_range, plot_filebase = plot_filebase, $ - save_path = save_path, savefilebase = savefilebase, $ + plot_filebase = plot_filebase, save_path = save_path, savefilebase = savefilebase, $ axis_type_1d = axis_type_1d, full_compare = all_type_pol for slice_i=0, compare_files.n_slices-1 do begin @@ -37,16 +35,19 @@ pro ps_difference_plots, folder_names, obs_info, ps_foldernames = ps_foldernames savefile_diff = compare_files.mid_savefile_2d[slice_i, cube_i] endif else begin - wt_cutoffs = ps_options.wt_cutoffs - wt_measures = ps_options.wt_measures if n_elements(ps_options) gt 1 then begin + wt_cutoffs = [ps_options[0].wt_cutoffs, ps_options[1].wt_cutoffs] + wt_measures = [ps_options[0].wt_measures, ps_options[1].wt_measures] if ps_options[1].wt_cutoffs eq ps_options[0].wt_cutoffs then begin wt_cutoffs = ps_options[0].wt_cutoffs endif if ps_options[1].wt_measures eq ps_options[0].wt_measures then begin wt_measures = ps_options[0].wt_measures endif - endif + endif else begin + wt_cutoffs = ps_options.wt_cutoffs + wt_measures = ps_options.wt_measures + endelse ps_differences, compare_files.input_savefile1[slice_i, cube_i], $ compare_files.input_savefile2[slice_i, cube_i], refresh = refresh_diff, $ diff --git a/ps_plotting/ps_ratio_plots.pro b/ps_plotting/ps_ratio_plots.pro index 2ed4de6..d554202 100644 --- a/ps_plotting/ps_ratio_plots.pro +++ b/ps_plotting/ps_ratio_plots.pro @@ -1,9 +1,9 @@ pro ps_ratio_plots, folder_names, obs_info, cube_types, pols, ps_foldernames=ps_foldernames, $ all_pol_diff_ratio = all_pol_diff_ratio, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, ps_options = ps_options, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ plot_slices = plot_slices, slice_type = slice_type, $ - freq_ch_range = freq_ch_range, plot_filebase = plot_filebase, save_path=save_path, $ + plot_filebase = plot_filebase, save_path=save_path, $ diff_ratio = diff_ratio, diff_range = diff_range, diff_min_abs = diff_min_abs, $ invert_colorbar = invert_colorbar, quiet = quiet, window_num = window_num @@ -15,10 +15,9 @@ pro ps_ratio_plots, folder_names, obs_info, cube_types, pols, ps_foldernames=ps_ compare_plot_prep, folder_names, obs_info, cube_types, pols, 'ratio', compare_files, $ ps_foldernames=ps_foldernames, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, ps_options = ps_options, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ plot_slices = plot_slices, slice_type = slice_type, fadd_2dbin = fadd_2dbin, $ - freq_ch_range = freq_ch_range, $ plot_filebase = plot_filebase, save_path = save_path, savefilebase = savefilebase, $ axis_type_1d = axis_type_1d, full_compare = all_pol_diff_ratio diff --git a/ps_setup/casa_file_setup.pro b/ps_setup/casa_file_setup.pro index b33f160..919f2fb 100644 --- a/ps_setup/casa_file_setup.pro +++ b/ps_setup/casa_file_setup.pro @@ -1,9 +1,8 @@ function casa_file_setup, filename, pol_inc, save_path = save_path, $ refresh_info = refresh_info, weight_savefilebase = weight_savefilebase_in, $ variance_savefilebase = variance_savefilebase_in, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ uvf_savefilebase = uvf_savefilebase_in, savefilebase = savefilebase_in, $ - uvf_options = uvf_options, ps_options = ps_options + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options if n_elements(pol_inc) ne 0 then pol_inc_in = pol_inc @@ -129,9 +128,8 @@ function casa_file_setup, filename, pol_inc, save_path = save_path, $ file_struct_arr = casa_file_setup(info_file, pol_inc, save_path = save_path, $ refresh_info = refresh_info, weight_savefilebase = weight_savefilebase_in, $ variance_savefilebase = variance_savefilebase_in, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ uvf_savefilebase = uvf_savefilebase_in, savefilebase = savefilebase_in, $ - uvf_options = uvf_options, ps_options = ps_options) + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options) return, file_struct_arr endif @@ -245,21 +243,40 @@ function casa_file_setup, filename, pol_inc, save_path = save_path, $ ncubes = npol * metadata_struct.ntypes if tag_exist(metadata_struct, 'no_var') then no_var = 1 else no_var = 0 - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin freq_mask = intarr(n_elements(metadata_struct.frequencies)) + 1 - freq_mask[freq_flags] = 0 + freq_mask[freq_options.freq_flags] = 0 wh_freq_use = where(freq_mask gt 0, count_freq_use) if count_freq_use lt 3 then message, 'Too many frequencies flagged' + + if min(freq_options.freq_flags) lt 0 or max(freq_options.freq_flags) gt n_elements(metadata_struct.frequencies) then begin + message, 'invalid freq_flags' + endif + + freq_options = create_freq_options(freq_options, freq_mask = freq_mask) + endif - if n_elements(freq_flags) ne 0 then begin - if min(freq_flags) lt 0 or max(freq_flags) gt n_elements(metadata_struct.frequencies) $ - then message, 'invalid freq_flags' + if tag_exist(freq_options, 'freq_ch_range') then begin + if min(freq_options.freq_ch_range) lt 0 $ + or max(freq_options.freq_ch_range) gt n_elements(metadata_struct.frequencies)-1 then begin + message, 'invalid freq_ch_range' + endif + if freq_options.freq_ch_range[1]-freq_options.freq_ch_range[0]+1 lt 3 then begin + message, 'freq_ch_range does not have enough frequencies' + endif + + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask_clipped = freq_mask[freq_options.freq_ch_range[0]:freq_options.freq_ch_range[1]] + wh_freq_use = where(freq_mask_clipped gt 0, count_freq_use) + if count_freq_use lt 3 then begin + message, 'Too few unflagged frequencies in freq_ch_range ' + endif + endif endif - file_tags = create_file_tags(freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, uvf_options = uvf_options, ps_options = ps_options) + file_tags = create_file_tags(uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options) wt_file_label = '_weights_' + strlowcase(pol_inc) file_label = '_' + strlowcase(metadata_struct.type_pol_str) @@ -406,8 +423,6 @@ function casa_file_setup, filename, pol_inc, save_path = save_path, $ if no_var then file_struct = create_struct(file_struct, 'no_var', 1) - if n_elements(freq_mask) gt 0 then file_struct = create_struct(file_struct, 'freq_mask', freq_mask) - if tag_exist(metadata_struct, 'vis_noise') gt 0 then file_struct = create_struct(file_struct, 'vis_noise', reform(metadata_struct.vis_noise[pol_index, *])) if i eq 0 then file_struct_arr = replicate(file_struct, ncubes) else file_struct_arr[i] = file_struct diff --git a/ps_setup/create_file_tags.pro b/ps_setup/create_file_tags.pro index 74bc61f..972ec03 100644 --- a/ps_setup/create_file_tags.pro +++ b/ps_setup/create_file_tags.pro @@ -1,5 +1,4 @@ -function create_file_tags, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, uvf_options = uvf_options, ps_options = ps_options +function create_file_tags, uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options window_type_list = ['Hann', 'Hamming', 'Blackman', 'Nutall', 'Blackman-Nutall', $ 'Blackman-Harris', 'Blackman-Harris^2', 'Tukey', 'None'] @@ -20,25 +19,43 @@ function create_file_tags, freq_ch_range = freq_ch_range, freq_flags = freq_flag sw_tag = '' endelse endelse - endif else sw_tag = '' + endif else begin + sw_tag = '' + endelse - if n_elements(freq_ch_range) ne 0 then begin - if min(freq_ch_range) lt 0 or max(freq_ch_range) - min(freq_ch_range) lt 3 then begin + if tag_exist(freq_options, 'freq_ch_range') then begin + if min(freq_options.freq_ch_range) lt 0 or max(freq_options.freq_ch_range) - min(freq_options.freq_ch_range) lt 3 then begin message, 'invalid freq_ch_range' endif - fch_tag = '_ch' + number_formatter(min(freq_ch_range)) + '-' +$ - number_formatter(max(freq_ch_range)) - endif else fch_tag = '' + freq_tag = '_ch' + number_formatter(min(freq_options.freq_ch_range)) + '-' +$ + number_formatter(max(freq_options.freq_ch_range)) + endif else begin + freq_tag = '' + endelse - if n_elements(freq_flags) ne 0 then begin - if n_elements(freq_flag_name) eq 0 then begin - freq_flag_name = '' + if tag_exist(freq_options, 'freq_flag_name') then begin + if size(freq_options.freq_flag_name, /type) ne 7 then begin + freq_flag_name_use = number_formatter(freq_options.freq_flag_name) endif else begin - if size(freq_flag_name, /type) ne 7 then freq_flag_name = number_formatter(freq_flag_name) + freq_flag_name_use = freq_options.freq_flag_name endelse - flag_tag = '_flag' + freq_flag_name - endif else flag_tag = '' - fch_tag = fch_tag + flag_tag + freq_tag = freq_tag + '_flag' + freq_flag_name_use + endif else begin + if tag_exist(freq_options, 'freq_flags') then begin + freq_tag = freq_tag + '_flag' + endif + endelse + + if freq_options.freq_avg_factor gt 1 then begin + freq_tag = freq_tag + '_freqave' + number_formatter(freq_options.freq_avg_factor) + if freq_options.force_even_freqs ne 0 then begin + freq_tag = freq_tag + '_evenfreqs' + endif + endif + + if tag_exist(freq_options, 'freq_bin_name') then begin + freq_tag = freq_tag + '_freqbin' + freq_options.freq_bin_name + endif if tag_exist(uvf_options, 'delta_uv_lambda') then begin uv_tag = '_deluv' + number_formatter(uvf_options.delta_uv_lambda) @@ -60,22 +77,43 @@ function create_file_tags, freq_ch_range = freq_ch_range, freq_flags = freq_flag if tag_exist(uvf_options, 'uv_img_clip') then begin uv_tag = uv_tag + '_uvimgclip' + number_formatter(uvf_options.uv_img_clip) endif - uvf_tag = uv_tag + fch_tag + uvf_tag = uv_tag + + if ps_options.std_power then begin + kcube_tag = '_stdp' + endif else begin + kcube_tag = '' + endelse - if ps_options.std_power then kcube_tag = '_stdp' else kcube_tag = '' if ps_options.freq_dft then begin kcube_tag = kcube_tag + '_dft' ;; add a tag for using the regularly spaced approximate z rather than true z - if ps_options.dft_z_use eq 'regular' then kcube_tag = kcube_tag + 'reg' + if ps_options.dft_z_use eq 'regular' then begin + kcube_tag = kcube_tag + 'reg' + endif + endif + + if tag_exist(ps_options, 'kz_use') then begin + kcube_tag = kcube_tag + '_kzuse' + if tag_exist(ps_options, 'kzuse_name') then begin + kcube_tag = kcube_tag + ps_options.kzuse_name + endif + endif + + if ps_options.inverse_covar_weight then begin + kcube_tag = kcube_tag + '_invcovar' + endif + if ps_options.ave_removal then begin + kcube_tag = kcube_tag + '_averemove' endif - if ps_options.inverse_covar_weight then kcube_tag = kcube_tag + '_invcovar' - if ps_options.ave_removal then kcube_tag = kcube_tag + '_averemove' kcube_tag = kcube_tag + sw_tag power_tag = kcube_tag - if ps_options.no_wtd_avg then power_tag = power_tag + '_nowtavg' + if ps_options.no_wtd_avg then begin + power_tag = power_tag + '_nowtavg' + endif - file_tags = {uvf_tag: uvf_tag, kcube_tag: kcube_tag, power_tag: power_tag} + file_tags = {uvf_tag: uvf_tag, freq_tag:freq_tag, kcube_tag: kcube_tag, power_tag: power_tag} return, file_tags end diff --git a/ps_setup/create_freq_options.pro b/ps_setup/create_freq_options.pro new file mode 100644 index 0000000..f875961 --- /dev/null +++ b/ps_setup/create_freq_options.pro @@ -0,0 +1,126 @@ +function create_freq_options, $ + freq_options = freq_options, $ + freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, $ + freq_mask = freq_mask, $ + freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, $ + force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, $ + freq_bin_name = freq_bin_name, $ + frequencies = frequencies, $ + n_vis_freq = n_vis_freq, $ + beam_int = beam_int, $ + new_freq_mask = new_freq_mask + + if keyword_set(freq_flag_repeat) and n_elements(freq_flags) gt 0 then begin + if size(freq_flag_repeat, /type) ne 2 and n_elements(freq_flag_repeat) ne 1 then begin + message, "freq_flag_repeat must be a scalar integer specifying the number of " $ + + "times to repeat the freq_flags for the full flag array." + endif + + ;; remake freq flag array + n_inital = n_elements(freq_flags) + n_total = n_inital * freq_flag_repeat + freq_flags = reform(rebin(reform(freq_flags, n_inital, 1), n_inital, freq_flag_repeat, /sample), n_total) + endif + + update_tags = list() + update_values = list() + + if n_elements(freq_options) eq 0 then begin + ;; default to no averaging + if n_elements(freq_avg_factor) eq 0 then begin + freq_avg_factor = 1 + endif + + ;; default to not forcing even frequency spacing after averaging + if n_elements(force_even_freqs) eq 0 then begin + force_even_freqs = 0 + endif + + freq_options = {freq_avg_factor: freq_avg_factor, force_even_freqs: force_even_freqs} + endif else begin + if n_elements(freq_avg_factor) gt 0 then begin + if size(freq_avg_factor, /type) ne 2 then begin + message, "freq_avg_factor must be an integer." + endif + update_tags.add, 'freq_avg_factor' + update_values.add, freq_avg_factor + endif + if n_elements(force_even_freqs) gt 0 then begin + update_tags.add, 'force_even_freqs' + update_values.add, force_even_freqs + endif + + endelse + + if n_elements(freq_ch_range) gt 0 then begin + update_tags.add, 'freq_ch_range' + update_values.add, freq_ch_range + endif + + if n_elements(freq_flags) gt 0 then begin + update_tags.add, 'freq_flags' + update_values.add, freq_flags + endif + + if n_elements(freq_mask) gt 0 then begin + update_tags.add, 'freq_mask' + update_values.add, freq_mask + endif + + if n_elements(freq_flag_name) gt 0 then begin + update_tags.add, 'freq_flag_name' + update_values.add, freq_flag_name + endif + + if n_elements(freq_avg_bins) gt 0 then begin + if freq_options.freq_avg_factor gt 1 then begin + message, "freq_avg_factor and freq_avg_bins cannot both be set" + endif else if freq_options.force_even_freqs gt 0 then begin + message, "cannot use force_even_freqs with freq_avg_bins" + endif else begin + update_tags.add, 'freq_avg_bins' + update_values.add, freq_avg_bins + endelse + endif + + if n_elements(freq_bin_name) gt 0 then begin + update_tags.add, 'freq_bin_name' + update_values.add, freq_bin_name + endif + + if n_elements(frequencies) gt 0 then begin + update_tags.add, 'frequencies' + update_values.add, frequencies + endif + + if n_elements(n_vis_freq) gt 0 then begin + update_tags.add, 'n_vis_freq' + update_values.add, n_vis_freq + endif + + if n_elements(beam_int) gt 0 then begin + update_tags.add, 'beam_int' + update_values.add, beam_int + endif + + if n_elements(new_freq_mask) gt 0 then begin + update_tags.add, 'new_freq_mask' + update_values.add, new_freq_mask + endif + + if n_elements(update_tags) gt 0 or keyword_set(return_new) then begin + update_tags = update_tags.toarray() + + new_freq_options = update_option_struct($ + freq_options, update_tags, update_values, return_new = return_new $ + ) + + return, new_freq_options + endif else begin + return, freq_options + endelse +end diff --git a/ps_setup/create_ps_options.pro b/ps_setup/create_ps_options.pro index 41d561f..c8e5093 100644 --- a/ps_setup/create_ps_options.pro +++ b/ps_setup/create_ps_options.pro @@ -2,7 +2,8 @@ function create_ps_options, ps_options = ps_options, $ ave_removal = ave_removal, wt_cutoffs = wt_cutoffs, $ wt_measures = wt_measures, spec_window_type = spec_window_type, $ no_spec_window = no_spec_window, allow_beam_approx = allow_beam_approx, $ - freq_dft = freq_dft, dft_z_use = dft_z_use, save_sum_cube = save_sum_cube, $ + save_sum_cube = save_sum_cube, freq_dft = freq_dft, $ + dft_z_use = dft_z_use, kz_use = kz_use, kzuse_name = kzuse_name, $ std_power = std_power, no_wtd_avg = no_wtd_avg, $ inverse_covar_weight = inverse_covar_weight, return_new = return_new @@ -21,6 +22,9 @@ function create_ps_options, ps_options = ps_options, $ ;; default to not turning off spectral windowing if n_elements(no_spec_window) eq 0 then no_spec_window = 0 + ;; default to not saving sum cube + if n_elements(save_sum_cube) eq 0 then save_sum_cube = 0 + ;; Default to using a frequency DFT. if n_elements(freq_dft) eq 0 then freq_dft = 1 @@ -109,6 +113,16 @@ function create_ps_options, ps_options = ps_options, $ update_values.add, dft_z_use endif + if n_elements(kz_use) gt 0 then begin + update_tags.add, 'kz_use' + update_values.add, kz_use + endif + + if n_elements(kzuse_name) gt 0 then begin + update_tags.add, 'kzuse_name' + update_values.add, kzuse_name + endif + if n_elements(update_tags) gt 0 or keyword_set(return_new) then begin update_tags = update_tags.toarray() diff --git a/ps_setup/create_refresh_options.pro b/ps_setup/create_refresh_options.pro index 8f6aa85..dfc5f38 100644 --- a/ps_setup/create_refresh_options.pro +++ b/ps_setup/create_refresh_options.pro @@ -1,6 +1,7 @@ function create_refresh_options, refresh_options = refresh_options, $ refresh_dft = refresh_dft, refresh_weight_dft = refresh_weight_dft, $ - refresh_beam = refresh_beam, refresh_kcube = refresh_kcube, refresh_ps = refresh_ps, $ + refresh_beam = refresh_beam, refresh_freq_select_avg = refresh_freq_select_avg, $ + refresh_kcube = refresh_kcube, refresh_ps = refresh_ps, $ refresh_binning = refresh_binning, refresh_info = refresh_info, $ return_new = return_new @@ -8,7 +9,8 @@ function create_refresh_options, refresh_options = refresh_options, $ if keyword_set(refresh_dft) then refresh_weight_dft = 1 if keyword_set(refresh_dft) then refresh_beam = 1 if keyword_set(refresh_dft) or keyword_set(refresh_weight_dft) or $ - keyword_set(refresh_beam) then refresh_kcube = 1 + keyword_set(refresh_beam) then refresh_freq_select_avg = 1 + if keyword_set(refresh_freq_select_avg) then refresh_kcube = 1 if keyword_set(refresh_kcube) then refresh_ps = 1 if keyword_set(refresh_ps) then refresh_binning = 1 @@ -17,13 +19,15 @@ function create_refresh_options, refresh_options = refresh_options, $ if n_elements(refresh_dft) eq 0 then refresh_dft = 0 if n_elements(refresh_weight_dft) eq 0 then refresh_weight_dft = 0 if n_elements(refresh_beam) eq 0 then refresh_beam = 0 + if n_elements(refresh_freq_select_avg) eq 0 then refresh_freq_select_avg = 0 if n_elements(refresh_kcube) eq 0 then refresh_kcube = 0 if n_elements(refresh_ps) eq 0 then refresh_ps = 0 if n_elements(refresh_binning) eq 0 then refresh_binning = 0 if n_elements(refresh_info) eq 0 then refresh_info = 0 refresh_options = {refresh_dft: refresh_dft, refresh_weight_dft: refresh_weight_dft, $ - refresh_beam: refresh_beam, refresh_kcube: refresh_kcube, refresh_ps: refresh_ps, $ + refresh_beam: refresh_beam, refresh_freq_select_avg:refresh_freq_select_avg, $ + refresh_kcube: refresh_kcube, refresh_ps: refresh_ps, $ refresh_binning:refresh_binning, refresh_info:refresh_info} return, refresh_options @@ -43,6 +47,10 @@ function create_refresh_options, refresh_options = refresh_options, $ update_tags.add, 'refresh_beam' update_values.add, refresh_beam endif + if n_elements(refresh_freq_select_avg) gt 0 then begin + update_tags.add, 'refresh_freq_select_avg' + update_values.add, refresh_freq_select_avg + endif if n_elements(refresh_kcube) gt 0 then begin update_tags.add, 'refresh_kcube' update_values.add, refresh_kcube diff --git a/ps_setup/create_uvf_options.pro b/ps_setup/create_uvf_options.pro index 0214487..e58b536 100644 --- a/ps_setup/create_uvf_options.pro +++ b/ps_setup/create_uvf_options.pro @@ -13,7 +13,9 @@ function create_uvf_options, uvf_options = uvf_options, $ update_values = list() if n_elements(uvf_options) eq 0 then begin ;; default to giving dft progress reports - if n_elements(no_dft_progress) eq 0 then no_dft_progress = 0 + if n_elements(no_dft_progress) eq 0 then begin + no_dft_progress = 0 + endif ;; Cannot set full_image with image_clip=0 error if this is the case. if n_elements(full_image) eq 1 and n_elements(image_clip) eq 1 then begin @@ -24,7 +26,9 @@ function create_uvf_options, uvf_options = uvf_options, $ ;; full_image means, use the full image with small uv pixels set by the image size ;; default to not using the full image to set the uv cell size - if n_elements(full_image) eq 0 then full_image = 0 + if n_elements(full_image) eq 0 then begin + full_image = 0 + endif ;; default to not requiring radec if n_elements(require_radec) eq 0 then require_radec = 0 diff --git a/ps_setup/fhd_file_setup.pro b/ps_setup/fhd_file_setup.pro index 153575e..70f5b60 100644 --- a/ps_setup/fhd_file_setup.pro +++ b/ps_setup/fhd_file_setup.pro @@ -2,11 +2,10 @@ function fhd_file_setup, filename, weightfile = weightfile, $ variancefile = variancefile, beamfile = beamfile, pixelfile = pixelfile, $ dirtyvar = dirtyvar, modelvar = modelvar, weightvar = weightvar, $ variancevar = variancevar, beamvar = beamvar, pixelvar = pixelvar, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, $ save_path = save_path, savefilebase = savefilebase_in, $ uvf_input = uvf_input, sim = sim, refresh_info = refresh_info, $ - uvf_options = uvf_options, ps_options = ps_options, pol_inc=pol_inc + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ + pol_inc = pol_inc ;; check to see if filename is an info file wh_info = strpos(filename, '_info.idlsave') @@ -361,11 +360,10 @@ function fhd_file_setup, filename, weightfile = weightfile, $ variancefile = variancefile, beamfile = beamfile, pixelfile = pixelfile, $ dirtyvar = dirtyvar, modelvar = modelvar, weightvar = weightvar, $ variancevar = variancevar, beamvar = beamvar, pixelvar = pixelvar, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, $ save_path = save_path, savefilebase = savefilebase_in, $ uvf_input = uvf_input, sim = sim, refresh_info = refresh_info, $ - uvf_options = uvf_options, ps_options = ps_options) + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ + pol_inc = pol_inc) return, file_struct_arr endif @@ -1134,26 +1132,46 @@ function fhd_file_setup, filename, weightfile = weightfile, $ pol_exist = stregex(file_basename(metadata_struct.datafile), '[xy][xy]', /boolean, /fold_case) - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin freq_mask = intarr(n_elements(metadata_struct.frequencies)) + 1 - freq_mask[freq_flags] = 0 + if n_elements(freq_options.freq_flags) eq n_elements(metadata_struct.frequencies) then begin + freq_mask = ~freq_options.freq_flags + endif else begin + freq_mask[freq_options.freq_flags] = 0 + endelse - if n_elements(freq_ch_range) gt 0 then begin - freq_mask = freq_mask[min(freq_ch_range):max(freq_ch_range)] - endif wh_freq_use = where(freq_mask gt 0, count_freq_use) if count_freq_use lt 3 then message, 'Too many frequencies flagged' + if min(freq_options.freq_flags) lt 0 or max(freq_options.freq_flags) gt n_elements(metadata_struct.frequencies) then begin + message, 'invalid freq_flags' + endif + + freq_options = create_freq_options(freq_options = freq_options, freq_mask = freq_mask) + endif - if n_elements(freq_flags) ne 0 then begin - if min(freq_flags) lt 0 or max(freq_flags) gt n_elements(metadata_struct.frequencies) then begin - message, 'invalid freq_flags' + if tag_exist(freq_options, 'freq_ch_range') then begin + if min(freq_options.freq_ch_range) lt 0 $ + or max(freq_options.freq_ch_range) gt n_elements(metadata_struct.frequencies)-1 then begin + message, 'invalid freq_ch_range' + endif + if freq_options.freq_ch_range[1]-freq_options.freq_ch_range[0]+1 lt 3 then begin + message, 'freq_ch_range does not have enough frequencies' + endif + + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask_clipped = freq_mask[freq_options.freq_ch_range[0]:freq_options.freq_ch_range[1]] + wh_freq_use = where(freq_mask_clipped gt 0, count_freq_use) + if count_freq_use lt 3 then begin + message, 'Too few unflagged frequencies in freq_ch_range ' + endif endif endif - file_tags = create_file_tags(freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, uvf_options = uvf_options, ps_options = ps_options) + file_tags = create_file_tags( $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options $ + ) wt_file_label = '_weights_' + strlowcase(metadata_struct.pol_inc) file_label = strarr(npol, ntypes) @@ -1161,29 +1179,33 @@ function fhd_file_setup, filename, weightfile = weightfile, $ file_label[i,*] = '_' + strlowcase(metadata_struct.type_inc) + '_' + $ strlowcase(metadata_struct.pol_inc[i]) endfor - savefilebase = metadata_struct.general_filebase + file_tags.uvf_tag + file_label + savefilebase = metadata_struct.general_filebase + file_tags.uvf_tag + file_tags.freq_tag + file_label + uvf_full_savefilebase = strarr(npol, nfiles, ntypes) uvf_savefilebase = strarr(npol, nfiles, ntypes) uvf_label = strarr(npol, nfiles, ntypes) for pol_i=0, npol-1 do begin for file_i=0, nfiles-1 do begin if max(pol_exist) eq 1 then begin - uvf_savefilebase[pol_i, file_i,*] = cgRootName(metadata_struct.datafile[pol_i, file_i]) + $ - file_tags.uvf_tag + '_' + metadata_struct.type_inc - endif else begin - uvf_savefilebase[pol_i, file_i,*] = cgRootName(metadata_struct.datafile[pol_i, file_i]) + $ - file_tags.uvf_tag + file_label[pol_i, *] - endelse + this_base = cgRootName(metadata_struct.datafile[pol_i, file_i]) + file_tags.uvf_tag + uvf_savefilebase[pol_i, file_i,*] = this_base + file_tags.freq_tag + '_' + metadata_struct.type_inc + uvf_full_savefilebase[pol_i, file_i,*] = this_base + '_' + metadata_struct.type_inc + endif else begin + this_base = cgRootName(metadata_struct.datafile[pol_i, file_i]) + file_tags.uvf_tag + uvf_savefilebase[pol_i, file_i,*] = this_base + file_tags.freq_tag + file_label[pol_i, *] + uvf_full_savefilebase[pol_i, file_i,*] = this_base + file_label[pol_i, *] + endelse uvf_label[pol_i, file_i,*] = metadata_struct.infile_label[file_i] + '_' + $ file_label[pol_i, *] endfor endfor if n_elements(size(uvf_savefilebase, /dim)) ne 3 then begin + uvf_full_savefilebase = reform(uvf_full_savefilebase, [npol, nfiles, ntypes]) uvf_savefilebase = reform(uvf_savefilebase, [npol, nfiles, ntypes]) endif - ;; add sw tag to general_filebase so that plotfiles havefile_tags.uvf_tag in them - general_filebase = metadata_struct.general_filebase + file_tags.uvf_tag + ;; add uvf_tag & freq_tag to general_filebase so that plotfiles have them + general_filebase = metadata_struct.general_filebase + file_tags.uvf_tag + file_tags.freq_tag subfolders = {data: 'data' + path_sep(), plots: 'plots' + path_sep(), $ uvf: 'uvf_cubes' + path_sep(), beams: 'beam_cubes' + path_sep(), $ @@ -1197,6 +1219,7 @@ function fhd_file_setup, filename, weightfile = weightfile, $ radec_file = uvf_path + metadata_struct.general_filebase + $ '_radec.idlsave' + uvf_full_savefile = uvf_path + uvf_full_savefilebase + '_uvf.idlsave' uvf_savefile = uvf_path + uvf_savefilebase + '_uvf.idlsave' uf_savefile = uvf_slice_path + uvf_savefilebase + '_uf_plane.idlsave' vf_savefile = uvf_slice_path + uvf_savefilebase + '_vf_plane.idlsave' @@ -1216,6 +1239,15 @@ function fhd_file_setup, filename, weightfile = weightfile, $ if tag_exist(metadata_struct, 'beamfile') then begin beam_savefile = froot + subfolders.data + subfolders.beams + uvf_savefilebase + '_beam2.idlsave' + beam_full_savefile = froot + subfolders.data + subfolders.beams + uvf_full_savefilebase + '_beam2.idlsave' + + if file_tags.freq_tag ne '' and tag_exist(freq_options, "freq_ch_range") eq 0 $ + and tag_exist(freq_options, "freq_flags") eq 0 then begin + ;; If the only freq option is freq_avg_factor then the freq-averaged beam is + ;; the same as the full beam + beam_savefile = beam_full_savefile + endif + endif kspace_path = froot + subfolders.data + subfolders.kspace @@ -1229,16 +1261,21 @@ function fhd_file_setup, filename, weightfile = weightfile, $ power_savefile = kspace_path + savefilebase + file_tags.power_tag + '_power.idlsave' fits_power_savefile = kspace_path + savefilebase + file_tags.power_tag + '_power.fits' + weight_full_savefilebase = strarr(npol, nfiles) weight_savefilebase = strarr(npol, nfiles) for pol_i=0, npol-1 do begin for file_i=0, nfiles-1 do begin if max(pol_exist) eq 1 then begin - weight_savefilebase[pol_i, file_i] = cgRootName(metadata_struct.weightfile[pol_i, file_i]) + $ - file_tags.uvf_tag + '_weights' - endif else begin - weight_savefilebase[pol_i, file_i] = cgRootName(metadata_struct.weightfile[pol_i, file_i]) + $ - file_tags.uvf_tag + wt_file_label[pol_i] - endelse + weight_full_savefilebase[pol_i, file_i] = cgRootName(metadata_struct.weightfile[pol_i, file_i]) + $ + file_tags.uvf_tag + weight_savefilebase[pol_i, file_i] = weight_full_savefilebase[pol_i, file_i] + file_tags.freq_tag + '_weights' + weight_full_savefilebase[pol_i, file_i] += '_weights' + endif else begin + weight_full_savefilebase[pol_i, file_i] = cgRootName(metadata_struct.weightfile[pol_i, file_i]) + $ + file_tags.uvf_tag + weight_savefilebase[pol_i, file_i] = weight_full_savefilebase[pol_i, file_i] + file_tags.freq_tag + wt_file_label[pol_i] + weight_full_savefilebase[pol_i, file_i] += wt_file_label[pol_i] + endelse endfor endfor @@ -1249,6 +1286,7 @@ function fhd_file_setup, filename, weightfile = weightfile, $ endfor uvf_weight_savefile = uvf_path + weight_savefilebase + '_uvf.idlsave' + uvf_weight_full_savefile = uvf_path + weight_full_savefilebase + '_uvf.idlsave' uf_weight_savefile = uvf_slice_path + weight_savefilebase + '_uf_plane.idlsave' vf_weight_savefile = uvf_slice_path + weight_savefilebase + '_vf_plane.idlsave' uv_weight_savefile = uvf_slice_path + weight_savefilebase + '_uv_plane.idlsave' @@ -1270,7 +1308,12 @@ function fhd_file_setup, filename, weightfile = weightfile, $ derived_uvf_varname = strarr(nfiles,2) endif else begin if healpix or not keyword_set(uvf_input) then begin - derived_uvf_inputfiles = uvf_savefile[pol_i, *, 0:1] + ;; if there's frequency selection or averaging, we want the full files for this + if file_tags.freq_tag ne '' then begin + derived_uvf_inputfiles = uvf_full_savefile[pol_i, *, 0:1] + endif else begin + derived_uvf_inputfiles = uvf_savefile[pol_i, *, 0:1] + endelse derived_uvf_varname = strarr(nfiles, 2) + 'data_cube' endif else begin derived_uvf_inputfiles = strarr(nfiles,2) @@ -1328,10 +1371,12 @@ function fhd_file_setup, filename, weightfile = weightfile, $ weight_savefilebase:reform(weight_savefilebase[pol_i, *]), $ derived_uvf_inputfiles:derived_uvf_inputfiles, $ derived_uvf_varname:derived_uvf_varname, $ - file_label:file_label[pol_i,type_i], $ + file_label:file_label[pol_i,type_i], $ uvf_label:reform(uvf_label[pol_i,*,type_i]), $ wt_file_label:wt_file_label[pol_i], $ - uvf_tag:file_tags.uvf_tag, kcube_tag:file_tags.kcube_tag, $ + uvf_tag:file_tags.uvf_tag, $ + freq_tag: file_tags.freq_tag, $ + kcube_tag:file_tags.kcube_tag, $ power_tag:file_tags.power_tag, $ type_pol_str:metadata_struct.type_pol_str[pol_i,type_i], $ pol_index:pol_i, type_index:type_i, pol:metadata_struct.pol_inc[pol_i], $ @@ -1340,6 +1385,12 @@ function fhd_file_setup, filename, weightfile = weightfile, $ if tag_exist(metadata_struct, 'beamfile') then begin file_struct = create_struct(file_struct, 'beamfile', reform(metadata_struct.beamfile[pol_i,*]), $ 'beamvar', metadata_struct.beam_varname[pol_i], 'beam_savefile', reform(beam_savefile[pol_i, *])) + if file_tags.freq_tag ne '' and (tag_exist(freq_options, "freq_ch_range") $ + or tag_exist(freq_options, "freq_flags")) then begin + ;; If the only freq option is freq_avg_factor then the freq-averaged beam is + ;; the same as the full beam + file_struct = create_struct(file_struct, 'beam_full_savefile', reform(beam_full_savefile[pol_i, *])) + endif endif if healpix or not keyword_set(uvf_input) then begin @@ -1347,6 +1398,11 @@ function fhd_file_setup, filename, weightfile = weightfile, $ file_struct = create_struct(file_struct, 'uvf_savefile', reform(uvf_savefile[pol_i,*,type_i]), $ 'uvf_weight_savefile', uvf_weight_savefile[pol_i,*]) + + if file_tags.freq_tag ne '' then begin + file_struct = create_struct(file_struct, 'uvf_full_savefile', reform(uvf_full_savefile[pol_i,*,type_i]), $ + 'uvf_weight_full_savefile', uvf_weight_full_savefile[pol_i,*]) + endif endif if healpix then file_struct = create_struct(file_struct, 'pixelfile', $ @@ -1355,10 +1411,6 @@ function fhd_file_setup, filename, weightfile = weightfile, $ if no_var then file_struct = create_struct(file_struct, 'no_var', 1) - if n_elements(freq_mask) gt 0 then begin - file_struct = create_struct(file_struct, 'freq_mask', freq_mask) - endif - if tag_exist(metadata_struct, 'vis_noise') gt 0 then begin file_struct = create_struct(file_struct, 'vis_noise', reform(metadata_struct.vis_noise[pol_i,*,*])) endif diff --git a/ps_setup/rts_file_setup.pro b/ps_setup/rts_file_setup.pro index 5acf372..9356e3d 100644 --- a/ps_setup/rts_file_setup.pro +++ b/ps_setup/rts_file_setup.pro @@ -1,9 +1,7 @@ function rts_file_setup, filename, save_path = save_path, refresh_info = refresh_info, $ - weight_savefilebase = weight_savefilebase_in, $ - variance_savefilebase = variance_savefilebase_in, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + weight_savefilebase = weight_savefilebase_in, variance_savefilebase = variance_savefilebase_in, $ uvf_savefilebase = uvf_savefilebase_in, savefilebase = savefilebase_in, $ - use_fhd_norm = use_fhd_norm, uvf_options = uvf_options, ps_options = ps_options + use_fhd_norm = use_fhd_norm, uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options ;; check to see if filename is an info file wh_info = strpos(filename, '_info.idlsave') @@ -281,9 +279,8 @@ function rts_file_setup, filename, save_path = save_path, refresh_info = refresh file_struct_arr = rts_file_setup(info_file, save_path = save_path, refresh_info = refresh_info, $ weight_savefilebase = weight_savefilebase_in, $ variance_savefilebase = variance_savefilebase_in, $ - freq_ch_range = freq_ch_range, freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ uvf_savefilebase = uvf_savefilebase_in, savefilebase = savefilebase_in, $ - use_fhd_norm = use_fhd_norm, uvf_options = uvf_options, ps_options = ps_options) + use_fhd_norm = use_fhd_norm, uvf_options = uvf_options, freq_options=freq_options, ps_options = ps_options) return, file_struct_arr endif @@ -556,23 +553,42 @@ function rts_file_setup, filename, save_path = save_path, refresh_info = refresh pol_exist = stregex(metadata_struct.datafile, '[xy][xy]', /boolean, /fold_case) - if n_elements(freq_flags) ne 0 then begin + if tag_exist(freq_options, 'freq_flags') then begin freq_mask = intarr(n_elements(metadata_struct.frequencies)) + 1 - freq_mask[freq_flags] = 0 + freq_mask[freq_options.freq_flags] = 0 wh_freq_use = where(freq_mask gt 0, count_freq_use) if count_freq_use lt 3 then message, 'Too many frequencies flagged' + if min(freq_options.freq_flags) lt 0 or max(freq_options.freq_flags) gt n_elements(metadata_struct.frequencies) then begin + message, 'invalid freq_flags' + endif + + freq_options = create_freq_options(freq_options, freq_mask = freq_mask) + endif - if n_elements(freq_flags) ne 0 then begin - if min(freq_flags) lt 0 or max(freq_flags) gt n_elements(metadata_struct.frequencies) then begin - message, 'invalid freq_flags' + if tag_exist(freq_options, 'freq_ch_range') then begin + if min(freq_options.freq_ch_range) lt 0 $ + or max(freq_options.freq_ch_range) gt n_elements(metadata_struct.frequencies)-1 then begin + message, 'invalid freq_ch_range' + endif + if freq_options.freq_ch_range[1]-freq_options.freq_ch_range[0]+1 lt 3 then begin + message, 'freq_ch_range does not have enough frequencies' + endif + + if tag_exist(freq_options, 'freq_flags') then begin + freq_mask_clipped = freq_mask[freq_options.freq_ch_range[0]:freq_options.freq_ch_range[1]] + wh_freq_use = where(freq_mask_clipped gt 0, count_freq_use) + if count_freq_use lt 3 then begin + message, 'Too few unflagged frequencies in freq_ch_range ' + endif endif endif - file_tags = create_file_tags(freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, uvf_options = uvf_options, ps_options = ps_options) + file_tags = create_file_tags( $ + uvf_options = uvf_options, freq_options=freq_options, ps_options = ps_options $ + ) if keyword_set(use_fhd_norm) then begin file_tags.kcube_tag = file_tags.kcube_tag + '_fhdnorm' @@ -748,8 +764,6 @@ function rts_file_setup, filename, save_path = save_path, refresh_info = refresh if no_var then file_struct = create_struct(file_struct, 'no_var', 1) - if n_elements(freq_mask) gt 0 then file_struct = create_struct(file_struct, 'freq_mask', freq_mask) - if tag_exist(metadata_struct, 'vis_noise') gt 0 then file_struct = create_struct(file_struct, 'vis_noise', reform(metadata_struct.vis_noise[pol_i,*,*])) if tag_exist(metadata_struct, 'beam_int') then file_struct = create_struct(file_struct, 'beam_int', reform(metadata_struct.beam_int[pol_i,*,*])) diff --git a/ps_utils/z_mpc.pro b/ps_utils/z_mpc.pro index 781daa8..4f679d1 100644 --- a/ps_utils/z_mpc.pro +++ b/ps_utils/z_mpc.pro @@ -1,5 +1,6 @@ -function z_mpc, frequencies, hubble_param = hubble_param, f_delta = f_delta, even_freq = even_freq, $ - redshifts = redshifts, comov_dist_los = comov_dist_los, z_mpc_delta = z_mpc_delta +function z_mpc, frequencies, hubble_param = hubble_param, f_delta = f_delta, $ + even_freq = even_freq, redshifts = redshifts, comov_dist_los = comov_dist_los, $ + z_mpc_delta = z_mpc_delta, z_mpc_length = z_mpc_length n_freq = n_elements(frequencies) @@ -47,6 +48,7 @@ function z_mpc, frequencies, hubble_param = hubble_param, f_delta = f_delta, eve z_mpc_delta = mean(nominal_comov_diffs) z_mpc_mean = mean(nominal_comov_dist_los) + z_mpc_length = max(nominal_comov_dist_los) - min(nominal_comov_dist_los) + z_mpc_delta endif else begin even_freq = 1 @@ -54,7 +56,7 @@ function z_mpc, frequencies, hubble_param = hubble_param, f_delta = f_delta, eve f_delta = double(mean(freq_diff)) ;take the mean to reduce precision errors z_mpc_delta = double(mean(comov_los_diff)) z_mpc_mean = double(mean(comov_dist_los)) - n_kz = n_freq + z_mpc_length = max(comov_dist_los) - min(comov_dist_los) + z_mpc_delta endelse diff --git a/ps_wrappers/ps_comp1d_wrapper.pro b/ps_wrappers/ps_comp1d_wrapper.pro index 5a999ed..6abaf0f 100644 --- a/ps_wrappers/ps_comp1d_wrapper.pro +++ b/ps_wrappers/ps_comp1d_wrapper.pro @@ -6,6 +6,9 @@ pro ps_comp1d_wrapper, folder_names_in, obs_names_in, $ image_clip = image_clip, ave_removal = ave_removal, $ freq_dft = freq_dft, dft_z_use = dft_z_use, std_power = std_power, $ all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ png = png, eps = eps, pdf = pdf, data_range = data_range, data_min_abs = data_min_abs, $ sim = sim, wt_cutoffs = wt_cutoffs, $ wt_measures = wt_measures, window_num = window_num, $ @@ -21,24 +24,18 @@ pro ps_comp1d_wrapper, folder_names_in, obs_names_in, $ kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $ kpar_range_kperppower = kpar_range_kperppower - if n_elements(folder_names_in) gt 2 then message, 'only 1 or 2 folder_names allowed' - if n_elements(folder_names_in) eq 0 then message, 'at least 1 folder name must be specified' - if n_elements(obs_names_in) gt 2 then message, 'only 1 or 2 obs_names_in allowed' - - folder_names = get_folder(folder_names_in, loc_name = loc_name, rts = rts, $ - dirty_folder = dirty_folder) - - if keyword_set(version_test) and n_elements(ps_foldername) eq 0 $ - and n_elements(folder_names_in) eq 1 and n_elements(obs_names_in) lt 2 then begin - git_info = git_info(ps_repository_dir()) - ps_foldernames = ['ps_master', 'ps_' + git_info.branch] - endif - - obs_info = ps_filenames(folder_names, obs_names_in, dirty_folder = dirty_folder, $ - exact_obsnames = exact_obsnames, rts = rts, sim = sim, uvf_input = uvf_input, $ - casa = casa, data_subdirs = data_subdirs, ps_foldernames = ps_foldernames, $ - save_paths = save_paths, plot_paths = plot_paths, refresh_info = refresh_info, $ - no_wtvar_rts = no_wtvar_rts) + compare_setup_structures, folder_names_in, obs_names_in, $ + ps_foldernames = ps_foldernames, version_test = version_test, $ + spec_window_types = spec_window_types, delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = max_uv_lambda, full_image = full_image, $ + image_clip = image_clip, ave_removal = ave_removal, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, std_power = std_power, $ + all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + folder_names = folder_names, obs_info = obs_info, uvf_options = uvf_options, $ + freq_options = freq_options, ps_options = ps_options if n_elements(set_krange_1dave) eq 0 and not keyword_set(sim) then begin ;; set some default kperp & kpar ranges for 1d averaging @@ -55,190 +52,14 @@ pro ps_comp1d_wrapper, folder_names_in, obs_names_in, $ kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $ kpar_range_kperppower = kpar_range_kperppower) - if n_elements(diff_plot_path) eq 0 then begin - if n_elements(diff_save_path) gt 0 then begin - diff_plot_path = diff_save_path + path_sep() + 'plots' + path_sep() - endif - endif - - wh_noinfo = where(obs_info.info_files eq '', count_noinfo) - if count_noinfo gt 0 then message, 'Info files are not all present' - - if n_elements(delta_uv_lambda) gt 1 then message, 'only 1 delta_uv_lambda allowed' - - if n_elements(max_uv_lambda) lt 2 and n_elements(full_image) lt 2 $ - and n_elements(image_clip) lt 2 then begin - - uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = max_uv_lambda, full_image = full_image, image_clip = image_clip) - - endif else begin - case n_elements(max_uv_lambda) of - 0: - 1: begin - mul0 = max_uv_lambda - mul1 = max_uv_lambda - end - 2: begin - mul0 = max_uv_lambda[0] - mul1 = max_uv_lambda[1] - end - else: message, 'only 1 or 2 max_uv_lambda values allowed' - endcase - - case n_elements(full_image) of - 0: - 1: begin - fi0 = full_image - fi1 = full_image - end - 2: begin - fi0 = full_image[0] - fi1 = full_image[1] - end - endcase - - case n_elements(image_clip) of - 0: - 1: begin - ic0 = image_clip - ic1 = image_clip - end - 2: begin - ic0 = image_clip[0] - ic1 = image_clip[1] - end - else: message, 'only 1 or 2 image_clip values allowed' - endcase - - uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = mul0, full_image = fi0, image_clip = ic0) - uvf_options1 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = mul1, full_image = fi1, image_clip = ic1) - endelse - - if n_elements(ave_removal) lt 2 and n_elements(wt_cutoffs) lt 2 and $ - n_elements(wt_measures) lt 2 and n_elements(spec_window_types) lt 2 and $ - n_elements(freq_dft) lt 2 and n_elements(dft_z_use) lt 2 and $ - n_elements(std_power) lt 2 then begin - - ps_options = create_ps_options(ave_removal = ave_removal, $ - wt_cutoffs = wt_cutoffs, wt_measures = wt_measures, $ - spec_window_type = spec_window_types, $ - freq_dft = freq_dft, dft_z_use = dft_z_use, $ - std_power = std_power) - - endif else begin - case n_elements(ave_removal) of - 0: - 1: begin - ar0 = ave_removal - ar1 = ave_removal - end - 2: begin - ar0 = ave_removal[0] - ar1 = ave_removal[1] - end - else: message, 'only 1 or 2 ave_removal values allowed' - endcase - - case n_elements(wt_cutoffs) of - 0: - 1: begin - wtc0 = wt_cutoffs - wtc1 = wt_cutoffs - end - 2: begin - wtc0 = wt_cutoffs[0] - wtc1 = wt_cutoffs[1] - end - else: message, 'only 1 or 2 wt_cutoffs allowed' - endcase - - case n_elements(wt_measures) of - 0: - 1: begin - wtm0 = wt_measures - wtm1 = wt_measures - end - 2: begin - wtm0 = wt_measures[0] - wtm1 = wt_measures[1] - end - else: message, 'only 1 or 2 wt_measures allowed' - endcase - - case n_elements(spec_window_types) of - 0: - 1: begin - spw0 = spec_window_types - spw1 = spec_window_types - end - 2: begin - spw0 = spec_window_types[0] - spw1 = spec_window_types[1] - end - else: message, 'only 1 or 2 spec_window_types allowed' - endcase - - case n_elements(freq_dft) of - 0: - 1: begin - dft0 = freq_dft - dft1 = freq_dft - end - 2: begin - dft0 = freq_dft[0] - dft1 = freq_dft[1] - end - else: message, 'only 1 or 2 freq_dft values allowed' - endcase - - case n_elements(dft_z_use) of - 0: - 1: begin - dftz0 = dft_z_use - dftz1 = dft_z_use - end - 2: begin - dftz0 = dft_z_use[0] - dftz1 = dft_z_use[1] - end - else: message, 'only 1 or 2 dft_z_use values allowed' - endcase - - case n_elements(std_power) of - 0: - 1: begin - sp0 = std_power - sp1 = std_power - end - 2: begin - sp0 = std_power[0] - sp1 = std_power[1] - end - else: message, 'only 1 or 2 std_power values allowed' - endcase - - ps_options0 = create_ps_options(ave_removal = ar0, wt_cutoffs = wtc0, $ - wt_measures = wtm0, spec_window_type = spw0, freq_dft = dft0, $ - dft_z_use = dftz0, std_power = sp0) - - ps_options1 = create_ps_options(ave_removal = ar1, wt_cutoffs = wtc1, $ - wt_measures = wtm1, spec_window_type = spw1, freq_dft = dft1, $ - dft_z_use = dftz1, std_power = sp1) - - ps_options = [ps_options0, ps_options1] - endelse - plot_options = create_plot_options(plot_path = diff_plot_path, $ png = png, eps = eps, pdf = pdf) ps_comp1d_plots, folder_names, obs_info, ps_foldernames = ps_foldernames, $ - cube_types, pols, uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, $ + cube_types, pols, uvf_options = uvf_options, freq_options = freq_options, $ ps_options = ps_options, plot_options = plot_options, $ binning_1d_options = binning_1d_options, names = names, $ - all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + all_type_pol = all_type_pol, $ save_path = diff_save_path, savefilebase = savefilebase, $ quiet = quiet, window_num = window_num diff --git a/ps_wrappers/ps_diff_wrapper.pro b/ps_wrappers/ps_diff_wrapper.pro index 12de38b..7ab2ebf 100644 --- a/ps_wrappers/ps_diff_wrapper.pro +++ b/ps_wrappers/ps_diff_wrapper.pro @@ -4,8 +4,13 @@ pro ps_diff_wrapper, folder_names_in, obs_names_in, $ spec_window_types = spec_window_types, delta_uv_lambda = delta_uv_lambda, $ max_uv_lambda = max_uv_lambda, full_image = full_image, $ image_clip = image_clip, ave_removal = ave_removal, $ - freq_dft = freq_dft, dft_z_use = dft_z_use, std_power = std_power, $ - all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, $ + kz_use = kz_use, kzuse_name = kzuse_name, $ + std_power = std_power, all_type_pol = all_type_pol, $ + freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ + freq_flag_name = freq_flag_name, freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, freq_bin_name = freq_bin_name, $ plot_slices = plot_slices, slice_type = slice_type, $ png = png, eps = eps, pdf = pdf, data_range = data_range, data_min_abs = data_min_abs, $ kperp_linear_axis = kperp_linear_axis, kpar_linear_axis = kpar_linear_axis, sim = sim, $ @@ -25,26 +30,24 @@ pro ps_diff_wrapper, folder_names_in, obs_names_in, $ ky_range_lambda_1dave = ky_range_lambda_1dave, $ kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $ kpar_range_kperppower = kpar_range_kperppower, $ - kperp_plot_range = kperp_plot_range, kperp_lambda_plot_range = kperp_lambda_plot_range - - if n_elements(folder_names_in) gt 2 then message, 'only 1 or 2 folder_names allowed' - if n_elements(folder_names_in) eq 0 then message, 'at least 1 folder name must be specified' - if n_elements(obs_names_in) gt 2 then message, 'only 1 or 2 obs_names_in allowed' - - folder_names = get_folder(folder_names_in, loc_name = loc_name, rts = rts, $ - dirty_folder = dirty_folder) - - if keyword_set(version_test) and n_elements(ps_foldername) eq 0 $ - and n_elements(folder_names_in) eq 1 and n_elements(obs_names_in) lt 2 then begin - git_info = git_info(ps_repository_dir()) - ps_foldernames = ['ps_master', 'ps_' + git_info.branch] - endif + kperp_plot_range = kperp_plot_range, kperp_lambda_plot_range = kperp_lambda_plot_range, $ + kpar_plot_range = kpar_plot_range - obs_info = ps_filenames(folder_names, obs_names_in, dirty_folder = dirty_folder, $ - exact_obsnames = exact_obsnames, rts = rts, sim = sim, uvf_input = uvf_input, $ - casa = casa, data_subdirs = data_subdirs, ps_foldernames = ps_foldernames, $ - save_paths = save_paths, plot_paths = plot_paths, refresh_info = refresh_info, $ - no_wtvar_rts = no_wtvar_rts) + compare_setup_structures, folder_names_in, obs_names_in, $ + ps_foldernames = ps_foldernames, version_test = version_test, $ + spec_window_types = spec_window_types, delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = max_uv_lambda, full_image = full_image, $ + image_clip = image_clip, ave_removal = ave_removal, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, std_power = std_power, $ + kz_use = kz_use, kzuse_name = kzuse_name, $ + all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, freq_bin_name = freq_bin_name, $ + diff_plot_path = diff_plot_path, diff_save_path = diff_save_path, $ + folder_names = folder_names, obs_info = obs_info, uvf_options = uvf_options, $ + freq_options = freq_options, ps_options = ps_options binning_2d_options = create_binning_2d_options(no_kzero = no_kzero, $ log_kpar = log_kpar, log_kperp = log_kperp, kpar_bin = kpar_bin, kperp_bin = kperp_bin) @@ -59,196 +62,21 @@ pro ps_diff_wrapper, folder_names_in, obs_names_in, $ kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $ kpar_range_kperppower = kpar_range_kperppower) - if n_elements(diff_plot_path) eq 0 then begin - if n_elements(diff_save_path) gt 0 then begin - diff_plot_path = diff_save_path + path_sep() + 'plots' + path_sep() - endif - endif - - wh_noinfo = where(obs_info.info_files eq '', count_noinfo) - if count_noinfo gt 0 then message, 'Info files are not all present' - - if n_elements(delta_uv_lambda) gt 1 then message, 'only 1 delta_uv_lambda allowed' - - if n_elements(max_uv_lambda) lt 2 and n_elements(full_image) lt 2 $ - and n_elements(image_clip) lt 2 then begin - - uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = max_uv_lambda, full_image = full_image, image_clip = image_clip) - - endif else begin - case n_elements(max_uv_lambda) of - 0: - 1: begin - mul0 = max_uv_lambda - mul1 = max_uv_lambda - end - 2: begin - mul0 = max_uv_lambda[0] - mul1 = max_uv_lambda[1] - end - else: message, 'only 1 or 2 max_uv_lambda values allowed' - endcase - - case n_elements(full_image) of - 0: - 1: begin - fi0 = full_image - fi1 = full_image - end - 2: begin - fi0 = full_image[0] - fi1 = full_image[1] - end - endcase - - case n_elements(image_clip) of - 0: - 1: begin - ic0 = image_clip - ic1 = image_clip - end - 2: begin - ic0 = image_clip[0] - ic1 = image_clip[1] - end - else: message, 'only 1 or 2 image_clip values allowed' - endcase - - uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = mul0, full_image = fi0, image_clip = ic0) - uvf_options1 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = mul1, full_image = fi1, image_clip = ic1) - endelse - - if n_elements(ave_removal) lt 2 and n_elements(wt_cutoffs) lt 2 and $ - n_elements(wt_measures) lt 2 and n_elements(spec_window_types) lt 2 and $ - n_elements(freq_dft) lt 2 and n_elements(dft_z_use) lt 2 and $ - n_elements(std_power) lt 2 then begin - - ps_options = create_ps_options(ave_removal = ave_removal, $ - wt_cutoffs = wt_cutoffs, wt_measures = wt_measures, $ - spec_window_type = spec_window_types, $ - freq_dft = freq_dft, dft_z_use = dft_z_use, $ - std_power = std_power) - - endif else begin - case n_elements(ave_removal) of - 0: - 1: begin - ar0 = ave_removal - ar1 = ave_removal - end - 2: begin - ar0 = ave_removal[0] - ar1 = ave_removal[1] - end - else: message, 'only 1 or 2 ave_removal values allowed' - endcase - - case n_elements(wt_cutoffs) of - 0: - 1: begin - wtc0 = wt_cutoffs - wtc1 = wt_cutoffs - end - 2: begin - wtc0 = wt_cutoffs[0] - wtc1 = wt_cutoffs[1] - end - else: message, 'only 1 or 2 wt_cutoffs allowed' - endcase - - case n_elements(wt_measures) of - 0: - 1: begin - wtm0 = wt_measures - wtm1 = wt_measures - end - 2: begin - wtm0 = wt_measures[0] - wtm1 = wt_measures[1] - end - else: message, 'only 1 or 2 wt_measures allowed' - endcase - - case n_elements(spec_window_types) of - 0: - 1: begin - spw0 = spec_window_types - spw1 = spec_window_types - end - 2: begin - spw0 = spec_window_types[0] - spw1 = spec_window_types[1] - end - else: message, 'only 1 or 2 spec_window_types allowed' - endcase - - case n_elements(freq_dft) of - 0: - 1: begin - dft0 = freq_dft - dft1 = freq_dft - end - 2: begin - dft0 = freq_dft[0] - dft1 = freq_dft[1] - end - else: message, 'only 1 or 2 freq_dft values allowed' - endcase - - case n_elements(dft_z_use) of - 0: - 1: begin - dftz0 = dft_z_use - dftz1 = dft_z_use - end - 2: begin - dftz0 = dft_z_use[0] - dftz1 = dft_z_use[1] - end - else: message, 'only 1 or 2 dft_z_use values allowed' - endcase - - case n_elements(std_power) of - 0: - 1: begin - sp0 = std_power - sp1 = std_power - end - 2: begin - sp0 = std_power[0] - sp1 = std_power[1] - end - else: message, 'only 1 or 2 std_power values allowed' - endcase - - ps_options0 = create_ps_options(ave_removal = ar0, wt_cutoffs = wtc0, $ - wt_measures = wtm0, spec_window_type = spw0, freq_dft = dft0, $ - dft_z_use = dftz0, std_power = sp0) - - ps_options1 = create_ps_options(ave_removal = ar1, wt_cutoffs = wtc1, $ - wt_measures = wtm1, spec_window_type = spw1, freq_dft = dft1, $ - dft_z_use = dftz1, std_power = sp1) - - ps_options = [ps_options0, ps_options1] - endelse - plot_options = create_plot_options(plot_path = diff_plot_path, $ png = png, eps = eps, pdf = pdf) plot_2d_options = create_plot_2d_options(kperp_linear_axis = kperp_linear_axis, $ kpar_linear_axis = kpar_linear_axis, $ kperp_plot_range = kperp_plot_range, kperp_lambda_plot_range = kperp_lambda_plot_range, $ + kpar_plot_range = kpar_plot_range, $ data_range = data_range, color_type = color_type) ps_difference_plots, folder_names, obs_info, ps_foldernames = ps_foldernames, $ - cube_types, pols, uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, $ + cube_types, pols, uvf_options = uvf_options, freq_options = freq_options, $ ps_options = ps_options, $ - plot_options = plot_options, plot_2d_options = plot_2d_options, $ binning_2d_options = binning_2d_options, binning_1d_options = binning_1d_options, $ - all_type_pol = all_type_pol, refresh_diff = refresh_diff, freq_ch_range = freq_ch_range, $ + plot_options = plot_options, plot_2d_options = plot_2d_options, $ + all_type_pol = all_type_pol, refresh_diff = refresh_diff, $ plot_slices = plot_slices, slice_type = slice_type, $ save_path = diff_save_path, savefilebase = savefilebase, $ plot_1d = plot_1d, axis_type_1d = axis_type_1d, $ diff --git a/ps_wrappers/ps_ratio_wrapper.pro b/ps_wrappers/ps_ratio_wrapper.pro index 21b3f8f..0fbb231 100644 --- a/ps_wrappers/ps_ratio_wrapper.pro +++ b/ps_wrappers/ps_ratio_wrapper.pro @@ -1,204 +1,56 @@ pro ps_ratio_wrapper, folder_names_in, obs_names_in, ps_foldernames=ps_foldernames, $ exact_obsnames = exact_obsnames, version_test = version_test, $ cube_types = cube_types, pols = pols, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, $ + kz_use = kz_use, kzuse_name = kzuse_name, $ all_pol_diff_ratio = all_pol_diff_ratio, freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, freq_bin_name = freq_bin_name, $ spec_window_types = spec_window_types, delta_uv_lambda = delta_uv_lambda, $ full_image = full_image, image_clip = image_clip, $ ave_removal = ave_removal, std_power = std_power, diff_ratio = diff_ratio, $ diff_range = diff_range, png = png, eps = eps, pdf = pdf, $ - data_range = data_range, $ + data_range = data_range, diff_min_abs=diff_min_abs, $ color_type = color_type, invert_colorbar = invert_colorbar, $ kperp_linear_axis = kperp_linear_axis, $ kpar_linear_axis = kpar_linear_axis, sim = sim, wt_cutoffs = wt_cutoffs, $ + kperp_plot_range = kperp_plot_range, kperp_lambda_plot_range = kperp_lambda_plot_range, $ + kpar_plot_range = kpar_plot_range, $ wt_measures = wt_measures, window_num = window_num, $ uvf_input = uvf_input, diff_save_path = diff_save_path, plot_path = diff_plot_path - if n_elements(folder_names_in) gt 2 then message, 'only 1 or 2 folder_names allowed' - if n_elements(folder_names_in) eq 0 then message, 'at least 1 folder name must be specified' - if n_elements(obs_names_in) gt 2 then message, 'only 1 or 2 obs_names_in allowed' - - folder_names = get_folder(folder_names_in, loc_name = loc_name, rts = rts, $ - dirty_folder = dirty_folder) - - if keyword_set(version_test) and n_elements(ps_foldername) eq 0 $ - and n_elements(folder_names_in) eq 1 and n_elements(obs_names_in) lt 2 then begin - git_info = git_info(ps_repository_dir()) - ps_foldernames = ['ps_master', 'ps_' + git_info.branch] - endif - - obs_info = ps_filenames(folder_names, obs_names_in, dirty_folder = dirty_folder, $ - exact_obsnames = exact_obsnames, rts = rts, sim = sim, uvf_input = uvf_input, $ - casa = casa, data_subdirs = data_subdirs, ps_foldernames = ps_foldernames, $ - save_paths = save_paths, plot_paths = plot_paths, refresh_info = refresh_info, $ - no_wtvar_rts = no_wtvar_rts) - - if n_elements(diff_plot_path) eq 0 then begin - if n_elements(diff_save_path) gt 0 then begin - diff_plot_path = diff_save_path + path_sep() + 'plots' + path_sep() - endif - endif - - wh_noinfo = where(obs_info.info_files eq '', count_noinfo) - if count_noinfo gt 0 then message, 'Info files are not all present' - - if n_elements(data_range) eq 0 then begin - if n_elements(color_type) gt 0 then begin - if color_type eq 'linear' then begin - data_range = [0, 1.2] - endif - endif - if n_elements(data_range) eq 0 then begin - data_range = [1e-3, 1e1] - endif - endif - - if n_elements(delta_uv_lambda) gt 1 then message, 'only 1 delta_uv_lambda allowed' - - if n_elements(max_uv_lambda) lt 2 and n_elements(full_image) lt 2 $ - and n_elements(image_clip) lt 2 then begin - - uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = max_uv_lambda, full_image = full_image, image_clip = image_clip) - - endif else begin - case n_elements(max_uv_lambda) of - 0: - 1: begin - mul0 = max_uv_lambda - mul1 = max_uv_lambda - end - 2: begin - mul0 = max_u_lambda[0] - mul1 = max_uv_lambda[1] - end - else: message, 'only 1 or 2 max_uv_lambda values allowed' - endcase - - case n_elements(full_image) of - 0: - 1: begin - fi0 = full_image - fi1 = full_image - end - 2: begin - fi0 = full_image[0] - fi1 = full_image[1] - end - else: message, 'only 1 or 2 full_image values allowed' - endcase - - case n_elements(image_clip) of - 0: - 1: begin - ic0 = image_clip - ic1 = image_clip - end - 2: begin - ic0 = image_clip[0] - ic1 = image_clip[1] - end - else: message, 'only 1 or 2 image_clip values allowed' - endcase - - uvf_options0 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = mul0, full_image = fi0, image_clip = ic0) - uvf_options1 = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ - max_uv_lambda = mul1, full_image = fi1, image_clip = ic1) - endelse - - if n_elements(ave_removal) lt 2 and n_elements(wt_cutoffs) lt 2 and $ - n_elements(wt_measures) lt 2 and n_elements(spec_window_types) lt 2 and $ - n_elements(std_power) lt 2 then begin - - ps_options = create_ps_options(ave_removal = ave_removal, wt_cutoffs = wt_cutoffs, $ - wt_measures = wt_measures, spec_window_type = spec_window_types, std_power = std_power) - - endif else begin - case n_elements(ave_removal) of - 0: - 1: begin - ar0 = ave_removal - ar1 = ave_removal - end - 2: begin - ar0 = ave_removal[0] - ar1 = ave_removal[1] - end - else: message, 'only 1 or 2 ave_removal values allowed' - endcase - - case n_elements(wt_cutoffs) of - 0: - 1: begin - wtc0 = wt_cutoffs - wtc1 = wt_cutoffs - end - 2: begin - wtc0 = wt_cutoffs[0] - wtc1 = wt_cutoffs[1] - end - else: message, 'only 1 or 2 wt_cutoffs allowed' - endcase - - case n_elements(wt_measures) of - 0: - 1: begin - wtm0 = wt_measures - wtm1 = wt_measures - end - 2: begin - wtm0 = wt_measures[0] - wtm1 = wt_measures[1] - end - else: message, 'only 1 or 2 wt_measures allowed' - endcase - - case n_elements(spec_window_types) of - 0: - 1: begin - spw0 = spec_window_types - spw1 = spec_window_types - end - 2: begin - spw0 = spec_window_types[0] - spw1 = spec_window_types[1] - end - else: message, 'only 1 or 2 spec_window_types allowed' - endcase - - case n_elements(std_power) of - 0: - 1: begin - sp0 = std_power - sp1 = std_power - end - 2: begin - sp0 = std_power[0] - sp1 = std_power[1] - end - else: message, 'only 1 or 2 std_power values allowed' - endcase - - ps_options0 = create_ps_options(ave_removal = ar0, wt_cutoffs = wtc0, $ - wt_measures = wtm0, spec_window_type = spw0, std_power = sp0) - - ps_options1 = create_ps_options(ave_removal = ar1, wt_cutoffs = wtc1, $ - wt_measures = wtm1, spec_window_type = spw1, std_power = sp1) - - ps_options = [ps_options0, ps_options1] - endelse + compare_setup_structures, folder_names_in, obs_names_in, $ + ps_foldernames = ps_foldernames, version_test = version_test, $ + spec_window_types = spec_window_types, delta_uv_lambda = delta_uv_lambda, $ + max_uv_lambda = max_uv_lambda, full_image = full_image, $ + image_clip = image_clip, ave_removal = ave_removal, $ + freq_dft = freq_dft, dft_z_use = dft_z_use, std_power = std_power, $ + kz_use = kz_use, kzuse_name = kzuse_name, $ + all_type_pol = all_type_pol, freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, freq_bin_name = freq_bin_name, $ + diff_plot_path = diff_plot_path, diff_save_path = diff_save_path, $ + folder_names = folder_names, obs_info = obs_info, uvf_options = uvf_options, $ + freq_options = freq_options, ps_options = ps_options plot_options = create_plot_options(plot_path = diff_plot_path, $ png = png, eps = eps, pdf = pdf) plot_2d_options = create_plot_2d_options(kperp_linear_axis = kperp_linear_axis, $ - kpar_linear_axis = kpar_linear_axis, data_range = data_range, color_type = color_type) + kpar_linear_axis = kpar_linear_axis, $ + kperp_plot_range = kperp_plot_range, kperp_lambda_plot_range = kperp_lambda_plot_range, $ + kpar_plot_range = kpar_plot_range, $ + data_range = data_range, color_type = color_type) ps_ratio_plots, folder_names, obs_info, cube_types, ps_foldernames=ps_foldernames, $ - pols, all_pol_diff_ratio = all_pol_diff_ratio, freq_ch_range = freq_ch_range, $ - uvf_options0 = uvf_options0, uvf_options1 = uvf_options1, ps_options = ps_options, $ + pols, all_pol_diff_ratio = all_pol_diff_ratio, $ + uvf_options = uvf_options, freq_options = freq_options, ps_options = ps_options, $ plot_options = plot_options, plot_2d_options = plot_2d_options, $ save_path = diff_save_path, plot_filebase = plot_filebase, $ - diff_ratio = diff_ratio, diff_range = diff_range, invert_colorbar = invert_colorbar, $ - quiet = quiet, window_num = window_num + diff_ratio = diff_ratio, diff_range = diff_range, diff_min_abs = diff_min_abs, $ + invert_colorbar = invert_colorbar, quiet = quiet, window_num = window_num end diff --git a/ps_wrappers/ps_wrapper.pro b/ps_wrappers/ps_wrapper.pro index 27c1a17..918f9a6 100644 --- a/ps_wrappers/ps_wrapper.pro +++ b/ps_wrappers/ps_wrapper.pro @@ -12,9 +12,12 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $ full_image = full_image, image_clip = image_clip, $ pol_inc = pol_inc, type_inc = type_inc, freq_ch_range = freq_ch_range, $ freq_flags = freq_flags, freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, freq_bin_name = freq_bin_name, $ allow_beam_approx = allow_beam_approx, uvf_input = uvf_input, uv_avg = uv_avg, $ uv_img_clip = uv_img_clip, freq_dft = freq_dft, dft_z_use = dft_z_use, $ - std_power = std_power, $ + kz_use = kz_use, kzuse_name = kzuse_name, std_power = std_power, $ inverse_covar_weight = inverse_covar_weight, ave_removal = ave_removal, $ no_wtd_avg = no_wtd_avg, norm_rts_with_fhd = norm_rts_with_fhd, $ use_weight_cutoff_sim = use_weight_cutoff_sim, $ @@ -201,8 +204,6 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $ message, 'save_slices cannot be set to 0 if plot_slices is set' endif - print,'datafile = ' + datafile - if keyword_set(set_data_ranges) then begin if n_elements(range_1d) eq 0 then begin if keyword_set(plot_1d_delta) then begin @@ -248,10 +249,13 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $ dft_fchunk = 1 - if keyword_set(plot_binning_hist) then refresh_binning = 1 + if keyword_set(plot_binning_hist) then begin + refresh_binning = 1 + endif refresh_options = create_refresh_options(refresh_dft = refresh_dft, $ refresh_beam = refresh_beam, refresh_ps = refresh_ps, refresh_kcube = refresh_ps, $ + refresh_freq_select_avg = refresh_ps, $ refresh_binning = refresh_binning, refresh_info = refresh_info) uvf_options = create_uvf_options(delta_uv_lambda = delta_uv_lambda, $ @@ -259,10 +263,21 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $ uv_avg = uv_avg, uv_img_clip = uv_img_clip, require_radec = require_radec, $ dft_fchunk = dft_fchunk, no_dft_progress = no_dft_progress) + freq_options = create_freq_options( $ + freq_ch_range = freq_ch_range, $ + freq_flags = freq_flags, $ + freq_flag_name = freq_flag_name, $ + freq_flag_repeat = freq_flag_repeat, $ + freq_avg_factor = freq_avg_factor, $ + force_even_freqs = force_even_freqs, $ + freq_avg_bins = freq_avg_bins, $ + freq_bin_name = freq_bin_name) + ps_options = create_ps_options(ave_removal = ave_removal, wt_cutoffs = wt_cutoffs, $ wt_measures = wt_measures, spec_window_type = spec_window_type, $ no_spec_window = no_spec_window, allow_beam_approx = allow_beam_approx, $ - save_sum_cube = save_sum_cube, freq_dft = freq_dft, dft_z_use = dft_z_use, $ + save_sum_cube = save_sum_cube, freq_dft = freq_dft, $ + dft_z_use = dft_z_use, kz_use = kz_use, kzuse_name = kzuse_name, $ std_power = std_power, no_wtd_avg = no_wtd_avg, $ inverse_covar_weight = inverse_covar_weight) @@ -305,8 +320,8 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $ plot_flat_1d = plot_flat_1d, no_text_1d = no_text_1d) ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $ - type_inc = type_inc, freq_ch_range = freq_ch_range, freq_flags = freq_flags, $ - freq_flag_name = freq_flag_name, uvf_input = uvf_input, no_evenodd = no_evenodd, $ + type_inc = type_inc, freq_options = freq_options, $ + uvf_input = uvf_input, no_evenodd = no_evenodd, $ rts = rts, norm_rts_with_fhd = norm_rts_with_fhd, casa = casa, sim = sim, $ fix_sim_input = fix_sim_input, save_path = save_path, $ savefilebase = savefilebase, cube_power_info = cube_power_info, $ @@ -318,18 +333,24 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $ plot_1d_options = plot_1d_options if not keyword_set(set_data_ranges) then begin - if n_elements(data_range) ne 0 then $ + if n_elements(data_range) ne 0 then begin print, 'data_range used: ', number_formatter(data_range, format = '(e7.1)') - if n_elements(sigma_range) ne 0 then $ + endif + if n_elements(sigma_range) ne 0 then begin print, 'sigma_range used: ', number_formatter(sigma_range, format = '(e7.1)') - if n_elements(nev_range) ne 0 then $ + endif + if n_elements(nev_range) ne 0 then begin print, 'nev_range used: ', number_formatter(nev_range, format = '(e7.1)') - if n_elements(nnr_range) ne 0 then $ + endif + if n_elements(nnr_range) ne 0 then begin print, 'nnr_range used: ', number_formatter(nnr_range, format = '(e7.1)') - if n_elements(snr_range) ne 0 then $ + endif + if n_elements(snr_range) ne 0 then begin print, 'snr_range used: ', number_formatter(snr_range, format = '(e7.1)') - if n_elements(noise_range) ne 0 then $ + endif + if n_elements(noise_range) ne 0 then begin print, 'noise_range used: ', number_formatter(noise_range, format = '(e7.1)') + endif endif end