EoRImaging / eppsilon

eppsilon - error propagated power spectrum with interleaved observed noise
BSD 2-Clause "Simplified" License
5 stars 4 forks source link

DFT option for the weights as well #103

Closed nicholebarry closed 4 years ago

nicholebarry commented 5 years ago

In ps_kcube.pro, there is an option to fft or dft the data. This is not extended to the weights. We could do this, and I've tested it in the LS framework. This does have an effect on the LS phase.

Note Core dump of what I was testing in the code, apologies for messiness. I was also testing the effect of FTing to the full specified kz band, and then chopping off the highest bin in an even freq case. Differences were minimal (I might post them later for fun).

      z_relative = reverse(comov_dist_los-min(comov_dist_los))
      freq_kz_arr_std = rebin(reform(kz_mpc_orig, 1, n_elements(kz_mpc_orig)),n_freq, n_elements(kz_mpc_orig)) * rebin(z_relative, n_freq, n_elements(kz_mpc_orig))
      cos_arr_std = cos(freq_kz_arr_std)
      sin_arr_std = sin(freq_kz_arr_std)
      sum_sigma2 = reform(sum_sigma2, n_kx*n_ky, n_freq)
      covar_cos_std = matrix_multiply(sum_sigma2, cos_arr_std^2d) * (z_mpc_delta)^2.
      covar_sin_std = matrix_multiply(sum_sigma2, sin_arr_std^2d) * (z_mpc_delta)^2.
      covar_cross_std = matrix_multiply(sum_sigma2, cos_arr_std*sin_arr_std) * (z_mpc_delta)^2.
      covar_cos_std = reform(covar_cos_std, n_kx, n_ky, n_freq)
      covar_sin_std = reform(covar_sin_std, n_kx, n_ky, n_freq)
      covar_cross_std = reform(covar_cross_std, n_kx, n_ky, n_freq)

      if trim_max_pos then begin
        ;; remove unmatched max k positive mode:
        covar_cos_std = covar_cos_std[*, *, 0:-2]
        covar_sin_std = covar_sin_std[*,*,0:-2]
        covar_cross_std = covar_cross_std[*,*,0:-2]
      endif

      a5_0 = covar_cos_std[*,*,where(n_val eq 0)]
      a5_n = (covar_cos_std[*, *, wh_pos] + covar_cos_std[*, *, reverse(wh_neg)])/2.
      ;b5_n = complex(0,1) * (sigma2_ft[*, *, wh_pos] - $
      ;  sigma2_ft[*, *, reverse(wh_neg)])/2.
      covar_cos_std_1 = dblarr(n_kx, n_ky, n_kz)
      covar_cos_std_1[*, *, 0] = a5_0
      covar_cos_std_1[*, *, 1:n_kz-1] = a5_n

      a5_0 = covar_sin_std[*,*,where(n_val eq 0)]
      a5_n = (covar_sin_std[*, *, wh_pos] + covar_sin_std[*, *, reverse(wh_neg)])/2.
      ;b5_n = complex(0,1) * (sigma2_ft[*, *, wh_pos] - $
      ;  sigma2_ft[*, *, reverse(wh_neg)])/2.
      covar_sin_std_1 = dblarr(n_kx, n_ky, n_kz)
      covar_sin_std_1[*, *, 0] = a5_0
      covar_sin_std_1[*, *, 1:n_kz-1] = a5_n

      a5_0 = covar_cross_std[*,*,where(n_val eq 0)]
      a5_n = (covar_cross_std[*, *, wh_pos] - covar_cross_std[*, *, reverse(wh_neg)])/2.
      ;b5_n = complex(0,1) * (covar_cross_std[*, *, wh_pos] - $
      ;  covar_cross_std[*, *, reverse(wh_neg)])/2.
      covar_cross_std_1 = dblarr(n_kx, n_ky, n_kz)
      covar_cross_std_1[*, *, 0] = a5_0
      covar_cross_std_1[*, *, 1:n_kz-1] = a5_n

      wh_0f = where(n_freq_contrib eq 0, count_0f)
      if count_0f gt 0 then begin
        covar_cos_std_1[wh_0f, *] = 0
        covar_sin_std_1[wh_0f, *] = 0
        covar_cross_std_1[wh_0f, *] = 0
      endif

      ;; drop pixels with less than 1/3 of the frequencies
      if count_fewfreq gt 0 then begin
        covar_cos_std_1 = temporary(covar_cos_std_1) * mask_fewfreq
        covar_sin_std_1 = temporary(covar_sin_std_1) * mask_fewfreq
        covar_cross_std_1 = temporary(covar_cross_std_1) * mask_fewfreq
      endif
      covar_cross = covar_cross_std_1
      covar_cos = covar_cos_std_1
      covar_sin = covar_sin_std_1