Calibration misc

fhd_core/calibration/fhd_struct_init_cal.pro

@@ -71,7 +71,7 @@ auto_scale=Fltarr(2)
 cal_struct={n_pol:n_pol, n_freq:n_freq, n_tile:n_tile, n_time:n_time, uu:u_loc, vv:v_loc, $
     auto_initialize:auto_initialize, max_iter:max_cal_iter, phase_iter:phase_fit_iter, $
     tile_A:tile_A, tile_B:tile_B, tile_names:tile_names, bin_offset:bin_offset, freq:freq, gain:gain_arr_ptr, $
-    adaptive_gain:use_adaptive_calibration_gain, base_gain:calibration_base_gain, $
+    gain_bandpass:Pointer_copy(gain_arr_ptr), adaptive_gain:use_adaptive_calibration_gain, base_gain:calibration_base_gain, $
     gain_residual:gain_residual, auto_scale:auto_scale, auto_params:auto_params, cross_phase:0.0, stokes_mix_phase:0.0, $
     min_cal_baseline:min_cal_baseline, max_cal_baseline:max_cal_baseline, n_vis_cal:n_vis_cal, $
     time_avg:cal_time_average, min_solns:min_cal_solutions, ref_antenna:ref_antenna, $

fhd_core/calibration/vis_cal_bandpass.pro

@@ -1,6 +1,7 @@
 FUNCTION vis_cal_bandpass,cal,obs,params,cal_remainder=cal_remainder,file_path_fhd=file_path_fhd,cable_bandpass_fit=cable_bandpass_fit,$
     bandpass_directory=bandpass_directory,tile_use=tile_use,calibration_bandpass_cable_exclude=calibration_bandpass_cable_exclude,$
-    cal_bp_transfer=cal_bp_transfer,uvfits_version=uvfits_version,uvfits_subversion=uvfits_subversion,auto_ratio_calibration=auto_ratio_calibration,Extra=extra
+    cal_bp_transfer=cal_bp_transfer,uvfits_version=uvfits_version,uvfits_subversion=uvfits_subversion,auto_ratio_calibration=auto_ratio_calibration,$
+    no_png=no_png,Extra=extra
   ;This function is version 1 of calibrating each group of tiles with similar cable lengths per observation.
 
   ;Extract needed elements from the input structures
@@ -177,7 +178,7 @@ FUNCTION vis_cal_bandpass,cal,obs,params,cal_remainder=cal_remainder,file_path_f
 
   ENDELSE
 
-  IF Keyword_Set(file_path_fhd) THEN bandpass_plots,obs,bandpass_arr,file_path_fhd=file_path_fhd,$
+  IF Keyword_Set(file_path_fhd) AND ~Keyword_Set(no_png) THEN bandpass_plots,obs,bandpass_arr,file_path_fhd=file_path_fhd,$
     cable_length_ref=cable_length_ref,tile_use_arr=tile_use_arr
 
   RETURN,cal_bandpass

fhd_core/calibration/vis_calibrate.pro

@@ -96,7 +96,7 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt
 
     ;; Option to transfer pre-made and unflagged model visbilities
     if keyword_set(model_transfer) then begin
-      vis_model_arr = vis_model_transfer(obs,model_transfer)
+      vis_model_arr = vis_model_transfer(obs,params,model_transfer)
     endif
 
     RETURN,vis_cal
@@ -199,13 +199,14 @@ FUNCTION vis_calibrate,vis_ptr,cal,obs,status_str,psf,params,jones,vis_weight_pt
     ENDIF
 
     cal_bandpass=vis_cal_bandpass(cal,obs,params,cal_remainder=cal_remainder,auto_ratio_calibration=auto_ratio_calibration,$
-      file_path_fhd=file_path_fhd,_Extra=extra)
+      file_path_fhd=file_path_fhd,no_png=no_png,_Extra=extra)
 
     IF Keyword_Set(calibration_polyfit) THEN BEGIN
       cal_polyfit=vis_cal_polyfit(cal_remainder,obs,auto_ratio=auto_ratio,_Extra=extra)
 
       cal=vis_cal_combine(cal_polyfit,cal_bandpass)
     ENDIF ELSE cal=cal_bandpass
+    cal.gain_bandpass = cal_bandpass.gain
 
     IF Keyword_Set(auto_ratio_calibration) THEN BEGIN
       cal=cal_auto_ratio(obs,cal,auto_ratio=auto_ratio,auto_tile_i=auto_tile_i,/remultiply)    

fhd_core/source_modeling/vis_model_transfer.pro

@@ -1,4 +1,7 @@
-function vis_model_transfer,obs,model_transfer
+function vis_model_transfer,obs,params,model_transfer
+
+  data_nbaselines = obs.nbaselines
+  data_n_time = obs.n_time
 
   ;; Option to transfer pre-made and unflagged model visbilities
   vis_model_arr=PTRARR(obs.n_pol,/allocate)
@@ -11,6 +14,102 @@ function vis_model_transfer,obs,model_transfer
     print, "Model visibilities transferred from " + transfer_name
   endfor
 
+  ;; Get the params file associated with the model visibilities
+  if file_test(model_transfer+'/'+obs.obsname+'_params.sav') then begin
+    params_model = getvar_savefile(model_transfer+'/'+obs.obsname+'_params.sav','params')
+  endif else begin
+    if file_test(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav') then begin
+      params_model = getvar_savefile(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav','params')
+    endif else message, 'No params file found in model transfer directory.'
+  endelse 
+  model_n_time = n_elements(uniq(params_model.time))
+  model_nbaselines = (uniq(params_model.time))[0]+1
+
+  if model_n_time NE data_n_time then begin
+  ;; Exclude flagged times from the model visibilities if the number of time steps do not match
+
+    data_Jdate = (*obs.baseline_info).Jdate
+    time_res_Jdate = (double(obs.time_res)/2.) / (24.*3600.) ;offset to the center of a time step
+    model_time = params_model.time[uniq(params_model.time)]
+    data_time = params.time[uniq(params.time)]
+    tolerance = 1E-5
+
+    ; Match times betweeen model and data
+    ; Option 1: You are working with a model made by FHD. Thus the timing convention in the params is the same
+    ; Option 2: You are working with a model made by WODEN above version 1.4. This is different from FHD convention
+    ;           by half a timestep.
+    ; Option 3: You are working with a model made by WODEN below version 1.4. Thus the timing convention in the 
+    ;           params is different and you need to convert the model Jdate to the data Jdate within precision.
+
+    ; Option 1
+    matched_times_opt1 = intarr(data_n_time)-1
+    for data_time_i=0, data_n_time-1 do begin
+      min_val= min(abs(model_time - data_time[data_time_i]), min_ind)
+      if min_val LT tolerance then matched_times_opt1[data_time_i] = min_ind
+    endfor
+
+    ; Option 2
+    matched_times_opt2 = intarr(data_n_time)-1
+    for data_time_i=0, data_n_time-1 do begin
+      min_val= min(abs(model_time - (data_time[data_time_i]-time_res_Jdate)), min_ind)
+      if min_val LT tolerance then matched_times_opt2[data_time_i] = min_ind
+    endfor
+
+    ; Option 3
+    matched_times_opt3 = intarr(data_n_time)-1
+    for data_time_i=0, data_n_time-1 do begin
+      min_val = min(abs(model_time - (data_Jdate[data_time_i]-time_res_Jdate)), min_ind)
+      if min_val LT tolerance then matched_times_opt3[data_time_i] = min_ind
+    endfor
+
+    ; Find which option was more successful in matching the times
+    ; Initialization of matched_times arrays above to -1 allows up to check if all data times are matched (including the 0th index)
+    ; Error if no matching times found or if matching times does not obviously cover all data times
+    temp = where(matched_times_opt1 NE -1, n_matched_opt1)
+    temp = where(matched_times_opt2 NE -1, n_matched_opt2)
+    temp = where(matched_times_opt3 NE -1, n_matched_opt3)
+    if (n_matched_opt1 EQ 0) and (n_matched_opt2 EQ 0) and (n_matched_opt3 EQ 0) then $
+      message, 'ERROR: No matching times found between transferred model and data.'
+
+    max_matched = max([n_matched_opt1, n_matched_opt2, n_matched_opt3],max_matched_ind)
+    CASE max_matched_ind OF
+      0: matched_times = matched_times_opt1
+      1: matched_times = matched_times_opt2
+      2: matched_times = matched_times_opt3
+    ENDCASE
+
+    temp = where(matched_times NE -1, n_matched)
+    if n_matched NE data_n_time then $
+      message, 'ERROR: Number of matching times found between transferred model and data does not equal the number of timesteps in data.'
+
+  endif else matched_times = indgen(model_n_time)
+
+  ;An arbitrary larger number than tiles to make an index array
+  baseline_mod = ULONG(obs.n_tile*2)
+  ;Create a baseline array. Could use params.baseline_arr, but there is no guarentee of model gen origin
+  data_baseline_index = params.antenna1*baseline_mod + params.antenna2
+  model_baseline_index = params_model.antenna1*baseline_mod + params_model.antenna2
+
+  ;Initialize matched model pointer array
+  matched_model = PTRARR(obs.n_pol,/allocate)
+  for pol_i=0, obs.n_pol-1 do begin
+    if size(*vis_model_arr[pol_i], /type) eq 6 then *matched_model[pol_i] = complex(FLTARR(obs.n_freq, data_nbaselines * data_n_time))
+    if size(*vis_model_arr[pol_i], /type) eq 9 then *matched_model[pol_i] = dcomplex(DBLARR(obs.n_freq, data_nbaselines * data_n_time))
+  endfor
+
+  ; In each matched timestep, match the baselines and fill a new, matched model array 
+  ; Unmatched times in the transferred model are not included in the matched model array 
+  for time_i=0, data_n_time-1 do begin
+    ;suba is the subset of indices in the first array that are also in the second.
+    ;subb is the subset of indices in the second array that are also in the first.
+    match, data_baseline_index[time_i*data_nbaselines:(time_i+1)*data_nbaselines-1], $
+      model_baseline_index[matched_times[time_i]*model_nbaselines:(matched_times[time_i]+1)*model_nbaselines-1], $
+      suba, subb
+    for pol_i=0, obs.n_pol-1 do (*matched_model[pol_i])[*,suba+time_i*(data_nbaselines)] = (*vis_model_arr[pol_i])[*,subb+matched_times[time_i]*model_nbaselines]
+  endfor
+
+  for pol_i=0, obs.n_pol-1 do vis_model_arr[pol_i] = Pointer_copy(matched_model[pol_i])
+
   return, vis_model_arr
 
 end

fhd_core/source_modeling/vis_source_model.pro

@@ -6,7 +6,7 @@ FUNCTION vis_source_model,skymodel, obs, status_str, psf, params, vis_weight_ptr
 
 ;; Option to transfer pre-made and unflagged model visbilities
 if keyword_set(model_transfer) then begin
-   vis_arr = vis_model_transfer(obs,model_transfer)
+   vis_arr = vis_model_transfer(obs,params,model_transfer)
    return, vis_arr
 endif
 

fhd_output/calibration_plots/plot_cals.pro

@@ -110,9 +110,9 @@ for file_i=0, page_num-1 do begin
             if page_num GT 1 then gains1_res=gains1_res[*,file_i*page_tile_max:(file_i+1)*page_tile_max<(size(gains1_res))[2]-1]
             gains1_orig=gains1_res+gains1
             gains1_res_abs=Abs(gains1_orig)-Abs(gains1)
-            max_amp_res = mean(abs([gains0_res_abs,gains1_res_abs]),/NAN) + 3*stddev(abs([gains0_res_abs,gains1_res_abs]),/NAN)
-        ENDIF ELSE max_amp_res = mean(abs(gains0_res_abs),/NAN) + 3*stddev(abs(gains0_res_abs),/NAN)
-    
+            max_amp_res = mean(abs([gains0_res_abs,gains1_res_abs]),/nan) + 3*stddev(abs([gains0_res_abs,gains1_res_abs]),/nan)
+        ENDIF ELSE max_amp_res = mean(abs(gains0_res_abs),/nan) + 3*stddev(abs(gains0_res_abs),/nan)
+
         ;plot amplitude residuals
         IF max_amp_res GT 0 THEN BEGIN
             plot_cals_sub,freq,gains0_orig,gains1_orig,filename=orig_amp_filename,tile_A=tile_A,tile_B=tile_B,tile_use=tile_use,tile_exist=tile_exist,tile_names=tile_names,$

fhd_output/calibration_plots/plot_cals_sub.pro

@@ -23,14 +23,14 @@ IF Keyword_Set(phase) THEN BEGIN
 ENDIF ELSE BEGIN
     IF N_Elements(yrange) NE 2 THEN BEGIN
         IF n_pol GT 1 THEN BEGIN
-            yrange=Minmax(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]))
+            yrange=Minmax(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/nan)
         ENDIF ELSE BEGIN
-            yrange=Minmax(abs(gains_A[*,tile_use_i]))
+            yrange=Minmax(abs(gains_A[*,tile_use_i]),/nan)
         ENDELSE
 
         IF (gain_type LE 5) OR Keyword_Set(real_vs_imaginary) THEN BEGIN
-            IF n_pol GT 1 THEN max_amp = mean(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/NAN) + 5*stddev(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/NAN) $
-                ELSE max_amp = Mean(abs(gains_A[*,tile_use_i]),/NAN) + 5*stddev(abs(gains_A[*,tile_use_i]),/NAN)
+            IF n_pol GT 1 THEN max_amp = mean(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/nan) + 5*stddev(abs([gains_A[*,tile_use_i],gains_B[*,tile_use_i]]),/nan) $
+                ELSE max_amp = Mean(abs(gains_A[*,tile_use_i]),/nan) + 5*stddev(abs(gains_A[*,tile_use_i]),/nan)
             amp_range=Floor(Alog10(max_amp))
             IF amp_range LT 0 THEN amp_digits=1. ELSE amp_digits=2.
             amp_test=max_amp/10.^(amp_range-amp_digits)
@@ -104,4 +104,4 @@ ENDFOR
 cgtext,.4,max(plot_pos[*,3]+height/4),obsname,/normal
 cgPS_Close,/png,Density=300,Resize=cal_plot_resize,/allow_transparent,/nomessage
 
-END
+END