diff --git a/fhd_core/calibration/fhd_struct_init_cal.pro b/fhd_core/calibration/fhd_struct_init_cal.pro index 4953a967..7fb0e455 100644 --- a/fhd_core/calibration/fhd_struct_init_cal.pro +++ b/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, $ diff --git a/fhd_core/calibration/vis_cal_bandpass.pro b/fhd_core/calibration/vis_cal_bandpass.pro index bd9f4703..e4d44d26 100644 --- a/fhd_core/calibration/vis_cal_bandpass.pro +++ b/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 diff --git a/fhd_core/calibration/vis_calibrate.pro b/fhd_core/calibration/vis_calibrate.pro index 4b3259ed..93113564 100644 --- a/fhd_core/calibration/vis_calibrate.pro +++ b/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) diff --git a/fhd_core/source_modeling/vis_model_transfer.pro b/fhd_core/source_modeling/vis_model_transfer.pro index 3f2aedf7..d590d2c0 100644 --- a/fhd_core/source_modeling/vis_model_transfer.pro +++ b/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 diff --git a/fhd_core/source_modeling/vis_source_model.pro b/fhd_core/source_modeling/vis_source_model.pro index 4af865b2..7e51c1b8 100644 --- a/fhd_core/source_modeling/vis_source_model.pro +++ b/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 diff --git a/fhd_output/calibration_plots/plot_cals.pro b/fhd_output/calibration_plots/plot_cals.pro index 967544c9..08123a8e 100644 --- a/fhd_output/calibration_plots/plot_cals.pro +++ b/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,$ diff --git a/fhd_output/calibration_plots/plot_cals_sub.pro b/fhd_output/calibration_plots/plot_cals_sub.pro index cd7ba840..1a794b4a 100644 --- a/fhd_output/calibration_plots/plot_cals_sub.pro +++ b/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 \ No newline at end of file