Add option to include all voxels in 1d power but with correct density normalizations

dictionary.md

@@ -155,6 +155,7 @@ These flags (all default=0) tell the code to redo parts of the analysis that wou
 
 **kpar_range_kperppower**: This is a 2D vector giving the range of k_parallel values to include in the 1D k_perpendicular power binning (resulting in a 1D power as a function k_perpendicular rather than k).
 
+**kperp_density_norm_diagnostic_1d**: This is a flag (valid values are 0/1, default=0) indicating that low weight uv pixels should be included in the 1D binning, but with the appropriate density normalization corrections applied.
 
 ## Outputs keywords and options:
 

ps_core/kspace_rebinning_1d.pro

@@ -4,6 +4,7 @@ function kspace_rebinning_1d, power, k1_mpc, k2_mpc, k3_mpc, k_edges_mpc, $
     nbins_1d = bin_hist, var_power_1d = var_power_1d, mean_var_1d = mean_var_1d, $
     bin_arr_3d = bin_arr_3d, noise_frac_3d = noise_frac_3d, edge_on_grid = edge_on_grid, $
     match_datta = match_datta, kpar_power = kpar_power, kperp_power = kperp_power, $
+    kperp_density_norm_diagnostic = kperp_density_norm_diagnostic, $
     kperp_density_measure = kperp_density_measure, kperp_density_cutoff = kperp_density_cutoff, $
     binning_1d_options = binning_1d_options, plot_options = plot_options, $
     hubble_param = hubble_param, kperp_lambda_conv
@@ -264,13 +265,32 @@ function kspace_rebinning_1d, power, k1_mpc, k2_mpc, k3_mpc, k_edges_mpc, $
       if count_kperp_dense eq 0 then begin
         print, 'no kperp values exceed kperp_density_cutoff, using full volume'
       endif else begin
-        bin_arr_3d = reform(bin_arr_3d, n_kx*n_ky, n_kz)
-        bin_arr_3d[wh_sparse, *] = 0
-        bin_arr_3d = reform(bin_arr_3d, n_kx, n_ky, n_kz)
+        if keyword_set(kperp_density_norm_diagnostic) then begin
+          ;; don't cut out low density regions, but do apply normalization factor to dense regions
+          power_use = reform(power_use, n_kx*n_ky, n_kz)
+          power_use[wh_kperp_dense, *] = power_use[wh_kperp_dense, *]/0.5
+          power_use = reform(power_use, n_kx, n_ky, n_kz)
+
+          if n_elements(noise_expval_use) gt 0 then begin
+            noise_expval_use = reform(noise_expval_use, n_kx*n_ky, n_kz)
+            noise_expval_use[wh_kperp_dense, *] = noise_expval_use[wh_kperp_dense, *]/0.5
+            noise_expval_use = reform(noise_expval_use, n_kx, n_ky, n_kz)
+          endif
+          if n_elements(weights_use) gt 0 then begin
+            weights_use = reform(weights_use, n_kx*n_ky, n_kz)
+            weights_use[wh_kperp_dense, *] = weights_use[wh_kperp_dense, *]/0.5
+            weights_use = reform(weights_use, n_kx, n_ky, n_kz)
+          endif
+        endif else begin
+          ;; cut out the low density regions and apply normalization to the rest
+          bin_arr_3d = reform(bin_arr_3d, n_kx*n_ky, n_kz)
+          bin_arr_3d[wh_sparse, *] = 0
+          bin_arr_3d = reform(bin_arr_3d, n_kx, n_ky, n_kz)
 
-        power_use = power_use/0.5
-        if n_elements(noise_expval_use) gt 0 then noise_expval_use = noise_expval_use/0.5
-        if n_elements(weights_use) gt 0 then weights_use = weights_use*0.5^2.
+          power_use = power_use/0.5
+          if n_elements(noise_expval_use) gt 0 then noise_expval_use = noise_expval_use/0.5
+          if n_elements(weights_use) gt 0 then weights_use = weights_use*0.5^2.
+        endelse
       endelse
     endif
     if max(bin_arr_3d) eq 0 then begin

ps_core/ps_binning.pro

@@ -125,6 +125,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         bin_arr_3d = bin_arr_3d, noise_frac_3d = noise_frac_3d, $
         wedge_amp = wedge_amp_use, kperp_density_measure = wt_meas_use, $
         kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -134,6 +135,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         bin_arr_3d = bin_arr_3d, var_power_1d = var_power_1d, mean_var_1d = mean_var_1d, $
         noise_frac_3d = noise_frac_3d, wedge_amp = wedge_amp_use, $
         kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -144,6 +146,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         mean_var_1d = new_mean_var_1d, noise_frac_3d = noise_frac_3d, $
         wedge_amp = wedge_amp_use, $
         kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -153,6 +156,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
           binned_noise_expval = noise_expval_1d, binned_weights = weights_1d, $
           bin_arr_3d = bin_arr_3d, wedge_amp = wedge_amp_use, $
           kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+          kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
           binning_1d_options = binning_1d_options, plot_options = plot_options, $
           hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -161,6 +165,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
           binned_noise_expval = noise_expval_1d, binned_weights = weights_1d, $
           bin_arr_3d = bin_arr_3d, wedge_amp = wedge_amp_use, $
           kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+          kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
           binning_1d_options = binning_1d_options, plot_options = plot_options, $
           hubble_param = hubble_param, kperp_lambda_conv)
       endif
@@ -305,6 +310,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
       binned_noise_expval = noise_expval_kpar, binned_weights = weights_1d, $
       bin_arr_3d = bin_arr_3d, /kpar_power, $
       kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+      kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
       binning_1d_options = kpar_binning_1d_options, plot_options = plot_options, $
       hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -314,6 +320,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
       var_power_1d = var_power_1d, mean_var_1d = mean_var_1d, $
       bin_arr_3d = bin_arr_3d, /kpar_power, $
       kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+      kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
       binning_1d_options = kpar_binning_1d_options, plot_options = plot_options, $
       hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -323,6 +330,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         binned_noise_expval = noise_expval_kpar, binned_weights = weights_1d, $
         bin_arr_3d = bin_arr_3d, /kpar_power, $
         kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = kpar_binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -331,6 +339,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         binned_noise_expval = noise_expval_kpar, binned_weights = weights_1d, $
         bin_arr_3d = bin_arr_3d, /kpar_power, $
         kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = kpar_binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
     endif
@@ -400,6 +409,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
       binned_noise_expval = noise_expval_kperp, binned_weights = weights_1d, $
       bin_arr_3d = bin_arr_3d, /kperp_power, $
       kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+      kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
       binning_1d_options = kperp_binning_1d_options, plot_options = plot_options, $
       hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -409,6 +419,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
       var_power_1d = var_power_1d, mean_var_1d = mean_var_1d, $
       bin_arr_3d = bin_arr_3d, /kperp_power, $
       kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+      kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
       binning_1d_options = kperp_binning_1d_options, plot_options = plot_options, $
       hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -418,6 +429,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         binned_noise_expval = noise_expval_kperp, binned_weights = weights_1d, $
         bin_arr_3d = bin_arr_3d, /kperp_power, $
         kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = kperp_binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
 
@@ -426,6 +438,7 @@ pro ps_binning, file_struct, sim = sim, freq_flags = freq_flags, $
         binned_noise_expval = noise_expval_kperp, binned_weights = weights_1d, $
         bin_arr_3d = bin_arr_3d, /kperp_power, $
         kperp_density_measure = wt_meas_use, kperp_density_cutoff = wt_cutoff_use, $
+        kperp_density_norm_diagnostic = binning_1d_options.kperp_density_norm_diagnostic, $
         binning_1d_options = kperp_binning_1d_options, plot_options = plot_options, $
         hubble_param = hubble_param, kperp_lambda_conv)
     endif

ps_core/ps_main_plots.pro

@@ -759,6 +759,7 @@ pro ps_main_plots, datafile, beamfiles = beamfiles, pol_inc = pol_inc, $
 
   fadd_1dbin = ''
   if binning_1d_options.log_k then fadd_1dbin = fadd_1dbin + '_logk'
+  if binning_1d_options.kperp_density_norm_diagnostic then fadd_1dbin = fadd_1dbin + '_dennorm'
 
   fadd_kpar_1d = fadd_1dbin
   fadd_kperp_1d = fadd_1dbin

ps_setup/create_binning_1d_options.pro

@@ -6,14 +6,22 @@ function create_binning_1d_options, binning_1d_options = binning_1d_options, $
     kx_range_1dave = kx_range_1dave, kx_range_lambda_1dave = kx_range_lambda_1dave, $
     ky_range_1dave = ky_range_1dave, ky_range_lambda_1dave = ky_range_lambda_1dave, $
     kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $
-    kpar_range_kperppower = kpar_range_kperppower, return_new = return_new
+    kpar_range_kperppower = kpar_range_kperppower, $
+    kperp_density_norm_diagnostic = kperp_density_norm_diagnostic, $
+    return_new = return_new
 
   update_tags = list()
   update_values = list()
   if n_elements(binning_1d_options) eq 0 then begin
     ;; default to linear bins
     if n_elements(log_k) eq 0 then log_k = 0
 
+    ;; default to not using kperp_density_norm_diagnostic
+    ;; which keeps around low density voxels but normalizes things properly
+    if n_elements(kperp_density_norm_diagnostic) eq 0 then begin
+      kperp_density_norm_diagnostic = 0
+    endif
+
     binning_1d_options = {log_k: log_k}
   endif else begin
     if n_elements(log_k) gt 0 then begin
@@ -146,6 +154,11 @@ function create_binning_1d_options, binning_1d_options = binning_1d_options, $
     update_values.add, kpar_range_kperppower
   endif
 
+  if n_elements(kperp_density_norm_diagnostic) gt 0 then begin
+    update_tags.add, 'kperp_density_norm_diagnostic'
+    update_values.add, kperp_density_norm_diagnostic
+  endif
+
   if isa(remove_tags) gt 0 then begin
     if (n_elements(remove_tags) eq 0) then begin
       ;; undefine doesn't work on lists for some reason.

ps_wrappers/ps_wrapper.pro

@@ -38,6 +38,7 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $
     ky_range_lambda_1dave = ky_range_lambda_1dave, $
     kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $
     kpar_range_kperppower = kpar_range_kperppower, $
+    kperp_density_norm_diagnostic_1d = kperp_density_norm_diagnostic_1d, $
     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, baseline_axis = baseline_axis, $
@@ -272,7 +273,8 @@ pro ps_wrapper, folder_name_in, obs_name, data_subdirs=data_subdirs, $
     kx_range_lambda_1dave = kx_range_lambda_1dave, ky_range_1dave = ky_range_1dave, $
     ky_range_lambda_1dave = ky_range_lambda_1dave, $
     kperp_range_lambda_kparpower = kperp_range_lambda_kparpower, $
-    kpar_range_kperppower = kpar_range_kperppower)
+    kpar_range_kperppower = kpar_range_kperppower, $
+    kperp_density_norm_diagnostic = kperp_density_norm_diagnostic_1d)
 
   plot_options = create_plot_options(hinv = hinv, plot_path = plot_path, $
     plot_filebase = plot_filebase, note = note, individual_plots = individual_plots, $