Loading metis_dark_uvda.pro +26 −22 Original line number Diff line number Diff line Loading @@ -15,22 +15,19 @@ function m_center_uv, matrix return, a end function correzione_uvda, image, dark function normalize_dark, data, dark ima_ndit1 = rebin(image, 1024, 1024, /sample) ima_ndit2 = rebin(dark, 1024, 1024, /sample) dark_bin = rebin(dark, 1024, 1024, /sample) data_bin = rebin(data, 1024, 1024, /sample) x = [m_angle_uv(ima_ndit2), m_center_uv(ima_ndit2)] y = [m_angle_uv(ima_ndit1), m_center_uv(ima_ndit1)] x = [m_angle_uv(dark_bin), m_center_uv(dark_bin)] y = [m_angle_uv(data_bin), m_center_uv(data_bin)] r = linfit(x, y, yfit = yfit) ima = r[0] + ima_ndit2 * r[1] dark_norm = r[0] + dark * r[1] s = size(image) ima = rebin(ima, s[1], s[2]) return, ima return, dark_norm end function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history Loading @@ -47,7 +44,7 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix tsensor = header.tsensor ; WARN - temporary patch to handle local l1 fits files ; if tsensor eq 0. then tsensor = -25. if tsensor eq 0. then tsensor = -25. journal, 'Dark-current correction:' journal, ' dit = ' + string(dit, format = '(I0)') + ' ms' Loading @@ -56,9 +53,9 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix journal, ' obj_cnt = ' + string(obj_cnt, format = '(I0)') journal, ' tsensor = ' + string(tsensor, format = '(F0)') + ' degC' flag_use_extrapolated = 1 flag_use_notfullset = 1 flag_normalize_dark = 0 flag_use_extrapolated = boolean(1) flag_use_notfullset = boolean(1) flag_normalize_dark = boolean(0) obt_available = !null Loading Loading @@ -108,15 +105,13 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix dark_image = float(readfits(cal_pack.path + dark_file, /silent)) dark_nobin = dark_image ; rebin the dark and take into account the exposure time correction of l1 images dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin dark_image = dark_image * ndit1 * ndit2 ; apply normalization if required if flag_normalize_dark and masking.contains('disabled', /fold) then $ dark_image = correzione_uvda(data, dark_image) dark_image = normalize_dark(data, dark_image) mask = where(data eq 0., count) data = data - dark_image Loading @@ -131,20 +126,29 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix dark_error = fltarr(1024, 1024) ; rebin it and take into account the exposure time correction dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin) dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin ; calculate the total uncertainty if isa(error) then error += (data/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 if not isa(error) then error = 0. error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. jump: s = where(data le 0., count) quality_matrix[s] = 0 if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Dark-current correction: ', ' ' + dark_file] history = [history, $ 'Dark-current correction: ', $ ' ' + dark_file, $ ' extrapolated dark = ' + dark[i].extrapol.tolower(), $ ' renormalized dark = ' + json_serialize(flag_normalize_dark)] journal, ' dark file = ' + dark_file journal, ' extrapolated dark = ' + dark[i].extrapol.tolower() journal, ' renormalized dark = ' + json_serialize(flag_normalize_dark) return, data end metis_dark_vlda.pro +13 −7 Original line number Diff line number Diff line Loading @@ -14,7 +14,7 @@ function metis_dark_vlda, data, header, cal_pack, error = error, quality_matrix gain = cal_pack.vl_channel.radiometry.gain.value ; WARN - temporary patch to handle local l1 fits files ; if tsensor eq 0. then tsensor = -30. if tsensor eq 0. then tsensor = -30. journal, 'Bias and dark-current correction:' journal, ' dit = ' + string(dit, format = '(F0)') + ' s' Loading Loading @@ -84,21 +84,27 @@ function metis_dark_vlda, data, header, cal_pack, error = error, quality_matrix bias_nobin = bias_image bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin bias_error = bias_image * sqrt(rebin((bias_error/bias_nobin)^2, header.naxis1, header.naxis2)/nbin) bias_error = sqrt(rebin(bias_error^2, header.naxis1, header.naxis2) * nbin)/nbin dark_nobin = dark_image dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin) dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin mask = where(data eq 0.) mask = where(data eq 0., count) data = data - (bias_image * ndit + dark_image * xposure) data[mask] = 0. if isa(error) then error += (data * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2 ; set values in masked areas equal to zero if count gt 0 then data[mask] = 0. if not isa(error) then error = 0. error += ((data > 0) * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. jump: s = where(data le 0., count) if isa(quality_matrix) then quality_matrix[s] = 0 if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Bias and dark-current corrections: ', ' ' + dark_file] Loading metis_flat_field.pro +11 −6 Original line number Diff line number Diff line Loading @@ -31,16 +31,21 @@ function metis_flat_field, data, header, cal_pack, error = error, quality_matrix ; WARN - error for binning must be checked ff_nobin = ff_image nbin = header.nbin ff_image = rebin(ff_image, header.naxis1, header.naxis2) ff_error = ff_image * sqrt(rebin((ff_error/ff_nobin)^2, header.naxis1, header.naxis2)/header.nbin) ff_error = sqrt(rebin(ff_error^2, header.naxis1, header.naxis2) * nbin)/nbin mask = where(ff_image le 0.) ff_image[mask] = 1. s = where(ff_image le 0., count) if count gt 0 then ff_image[s] = 1. data = data/ff_image data[mask] = 0. if isa(error) then error += (ff_error/ff_image)^2 if count gt 0 then data[s] = 0. if not isa(error) then error = 0. error += (ff_error/ff_image)^2 if isa(quality_matrix) then $ if count gt 0 then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Flat-field correction: ', ' ' + ff_file] Loading metis_l2_prep_vl_polariz.pro +62 −35 Original line number Diff line number Diff line Loading @@ -44,6 +44,7 @@ pro metis_l2_prep_vl_polariz ; calibration block data = !null abs_error = !null data_header = !null data_subdark = !null Loading Loading @@ -123,13 +124,13 @@ pro metis_l2_prep_vl_polariz ; read and update the quality matrix quality_matrix *= mrdfits(input.file_name[k], 'quality matrix', /silent) image_quality_matrix = mrdfits(input.file_name[k], 'quality matrix', /silent) ; apply dark correction to compute stokes i and total brightness tb_history = !null image_subdark = metis_dark_vlda(image, header, cal_pack, history = tb_history) image_subdark = metis_dark_vlda(image, header, error = image_error, quality_matrix = image_quality_matrix, cal_pack, history = tb_history) ; ==================================== ; for data already subtracted of dark Loading @@ -141,6 +142,9 @@ pro metis_l2_prep_vl_polariz ; ==================================== data_subdark = [[[data_subdark]], [[image_subdark]]] abs_error = [[[abs_error]], [[image_subdark * sqrt(image_error)]]] quality_matrix *= image_quality_matrix endfor ; adjust header keywords Loading Loading @@ -187,6 +191,7 @@ pro metis_l2_prep_vl_polariz stokes_name = ['I', 'Q', 'U'] stokes = dblarr(header.naxis1, header.naxis2, 3) stokes_subdark = dblarr(header.naxis1, header.naxis2, 3) stokes_abs_error = dblarr(header.naxis1, header.naxis2, 3) for i = 0, 2 do begin journal, ' stokes = ' + stokes_name[i] Loading Loading @@ -256,6 +261,14 @@ pro metis_l2_prep_vl_polariz demod_tensor[*, *, 1] * data_subdark[*, *, 1] + $ demod_tensor[*, *, 2] * data_subdark[*, *, 2] + $ demod_tensor[*, *, 3] * data_subdark[*, *, 3] ; compute the stokes errors stokes_abs_error[*, *, i] = sqrt( $ demod_tensor[*, *, 0]^2 * abs_error[*, *, 0]^2 + $ demod_tensor[*, *, 1]^2 * abs_error[*, *, 1]^2 + $ demod_tensor[*, *, 2]^2 * abs_error[*, *, 2]^2 + $ demod_tensor[*, *, 3]^2 * abs_error[*, *, 3]^2) endfor demod_history = 'Demodulation performed for angles ' + string(angles, format = '(3(f0.1, ", "), f0.1)') + ' deg' Loading @@ -263,23 +276,34 @@ pro metis_l2_prep_vl_polariz tb_history = [tb_history, demod_history] pb_history = demod_history i = reform(stokes_subdark[*, *, 0]) i_abs_error = reform(stokes_abs_error[*, *, 0]) q = reform(stokes[*, *, 1]) q_abs_error = reform(stokes_abs_error[*, *, 1]) u = reform(stokes[*, *, 2]) u_abs_error = reform(stokes_abs_error[*, *, 2]) ; compute the tb from the dark-subtracted stokes i and apply other calibrations journal, 'Calibrating total brightness...' tb_image = reform(stokes_subdark[*, *, 0]) tb_image = metis_flat_field(tb_image, header, cal_pack, history = tb_history) tb_image = metis_vignetting(tb_image, header, cal_pack, history = tb_history) tb_image = metis_rad_cal(tb_image, header, cal_pack, /polarimetric, history = tb_history) tb_image = i tb_error = (i_abs_error/i)^2 tb_quality_matrix = quality_matrix tb_image = metis_flat_field(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history) tb_image = metis_vignetting(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history) tb_image = metis_rad_cal(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, /polarimetric, history = tb_history) ; compute the pb from the stokes q and u and apply other calibrations journal, 'Calibrating polarized brightness...' pb_image = sqrt(reform(stokes[*, *, 1])^2 + reform(stokes[*, *, 2])^2) pb_image = metis_flat_field(pb_image, header, cal_pack, history = pb_history) pb_image = metis_vignetting(pb_image, header, cal_pack, history = pb_history) pb_image = metis_rad_cal(pb_image, header, cal_pack, /polarimetric, history = pb_history) pb_image = sqrt(q^2 + u^2) pb_error = (q^2 * q_abs_error^2 + u^2 * u_abs_error^2)/pb_image^4 pb_quality_matrix = quality_matrix pb_image = metis_flat_field(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history) pb_image = metis_vignetting(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history) pb_image = metis_rad_cal(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, /polarimetric, history = pb_history) ; ==================================== ; for simple radiometric calibration Loading @@ -288,7 +312,10 @@ pro metis_l2_prep_vl_polariz ; compute the polarization angle from the stokes q and u pol_angle = 0.5D0 * atan(stokes[*, *, 2], stokes[*, *, 1]) * !radeg pol_angle = 0.5D0 * atan(stokes[*, *, 2], stokes[*, *, 1]) pol_angle_error = (q^2 * u_abs_error^2 + u^2 * q_abs_error^2)/(2D0 * pb_image^4)/pol_angle^2 pol_angle = pol_angle * !radeg journal, 'Polarization angle correctly computed.' Loading Loading @@ -422,23 +449,23 @@ pro metis_l2_prep_vl_polariz fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts' fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value' fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel' if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, quality_matrix mwrfits, float(quality_matrix), out_file_name[0], extension_header, /no_comment, /silent if not ref_detector then pb_quality_matrix = metis_rectify(pb_quality_matrix, 'VL') fits_add_checksum, extension_header, pb_quality_matrix mwrfits, float(pb_quality_matrix), out_file_name[0], extension_header, /no_comment, /silent journal, 'Quality-matrix extension correctly added.' ; add the extension with the error matrix if not ref_detector then pb_error = metis_rectify(pb_error, 'VL') error_matrix = pb_image * 0. ;pb_image * sqrt(pb_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[0], extension_header, /no_comment, /silent Loading Loading @@ -489,29 +516,29 @@ pro metis_l2_prep_vl_polariz fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Pixel quality' fxaddpar, extension_header, 'BUNIT', 'None' fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan) fxaddpar, extension_header, 'DATAMIN', min(tb_quality_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(tb_quality_matrix, /nan) fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:' fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts' fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value' fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel' if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, quality_matrix mwrfits, float(quality_matrix), out_file_name[1], extension_header, /no_comment, /silent if not ref_detector then tb_quality_matrix = metis_rectify(tb_quality_matrix, 'VL') fits_add_checksum, extension_header, tb_quality_matrix mwrfits, float(tb_quality_matrix), out_file_name[1], extension_header, /no_comment, /silent journal, 'Quality-matrix extension correctly added.' ; add the extension with the error matrix if not ref_detector then tb_error = metis_rectify(tb_error, 'VL') error_matrix = tb_image * 0. ;tb_image * sqrt(tb_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[1], extension_header, /no_comment, /silent Loading Loading @@ -568,23 +595,23 @@ pro metis_l2_prep_vl_polariz fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts' fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value' fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel' if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, quality_matrix mwrfits, float(quality_matrix), out_file_name[2], extension_header, /no_comment, /silent if not ref_detector then pol_angle_quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, pol_angle_quality_matrix mwrfits, float(pol_angle_quality_matrix), out_file_name[2], extension_header, /no_comment, /silent journal, 'Quality-matrix extension correctly added.' ; add the extension with the error matrix if not ref_detector then pol_angle_error = metis_rectify(pol_angle_error, 'VL') error_matrix = pol_angle * 0. ;pol_angle * sqrt(pol_angle_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[2], extension_header, /no_comment, /silent Loading Loading @@ -718,15 +745,15 @@ pro metis_l2_prep_vl_polariz ; add the extension with the error matrix if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') error_matrix = tb_image * 0. ;tb_image * sqrt(tb_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[3], extension_header, /no_comment, /silent Loading metis_rad_cal.pro +17 −14 Original line number Diff line number Diff line Loading @@ -11,6 +11,7 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor if keyword_set(polarimetric) then pmp_factor = 2. else pmp_factor = 1. nbin = header.nbin ndit = header.ndit endif Loading @@ -25,26 +26,26 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor pmp_factor = 1. nbin = header.nbin ndit = header.ndit1 * header.ndit2 ; read the uv/vl ratio of the door retro-illumination, to be used as correction for the optical efficiency efficiency_file = radiometry.opt_efficiency.file_name efficiency = float(readfits(cal_pack.path + efficiency_file, /silent)) eff_file = radiometry.opt_efficiency.file_name eff = float(readfits(cal_pack.path + eff_file, /silent)) eff_nobin = efficiency efficiency = rebin(efficiency, header.naxis1, header.naxis2) eff = rebin(eff, header.naxis1, header.naxis2) eff_error_file = radiometry.opt_efficiency.error eff_error = float(readfits(cal_pack.path + eff_error_file, /silent)) eff_error = efficiency * sqrt(rebin((eff_error/eff_nobin)^2, header.naxis1, header.naxis2)/header.nbin) eff_error = sqrt(rebin(eff_error^2, header.naxis1, header.naxis2) * nbin)/nbin s = where(efficiency le 0., count) quality_matrix[s] = 0 s = where(eff le 0., count) if count gt 0 then quality_matrix[s] = 0 ; WARN - temporary patch to bypass efficiency correction of the uv channel ; efficiency_file = 'not applied' ; efficiency = 1. ; eff_file = 'not applied' ; eff = 1. endif angular_pixel = channel.angular_pixel.value * header.nbin Loading Loading @@ -99,7 +100,8 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor data = data * cal_factor if isa(error) then error += cal_error^2 if not isa(error) then error = 0. error += cal_error^2 if header.filter.contains('VL', /fold) then begin history = [history, ' 1 MSB = ' + string(unit_factor, format = '(E8.2)') + ' ph/cm2/s/sr'] Loading @@ -109,13 +111,14 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor ; efficiency correction data = data/abs(efficiency) data = data/abs(eff) if isa(error) then error += (eff_error/efficiency)^2 if not isa(error) then error = 0. error += (eff_error/eff)^2 history = [history, ' efficiency corr. map = ' + efficiency_file] history = [history, ' efficiency corr. map = ' + eff_file] journal, ' UV efficiency corr. map = ' + efficiency_file journal, ' UV efficiency corr. map = ' + eff_file endif return, data Loading Loading
metis_dark_uvda.pro +26 −22 Original line number Diff line number Diff line Loading @@ -15,22 +15,19 @@ function m_center_uv, matrix return, a end function correzione_uvda, image, dark function normalize_dark, data, dark ima_ndit1 = rebin(image, 1024, 1024, /sample) ima_ndit2 = rebin(dark, 1024, 1024, /sample) dark_bin = rebin(dark, 1024, 1024, /sample) data_bin = rebin(data, 1024, 1024, /sample) x = [m_angle_uv(ima_ndit2), m_center_uv(ima_ndit2)] y = [m_angle_uv(ima_ndit1), m_center_uv(ima_ndit1)] x = [m_angle_uv(dark_bin), m_center_uv(dark_bin)] y = [m_angle_uv(data_bin), m_center_uv(data_bin)] r = linfit(x, y, yfit = yfit) ima = r[0] + ima_ndit2 * r[1] dark_norm = r[0] + dark * r[1] s = size(image) ima = rebin(ima, s[1], s[2]) return, ima return, dark_norm end function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history Loading @@ -47,7 +44,7 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix tsensor = header.tsensor ; WARN - temporary patch to handle local l1 fits files ; if tsensor eq 0. then tsensor = -25. if tsensor eq 0. then tsensor = -25. journal, 'Dark-current correction:' journal, ' dit = ' + string(dit, format = '(I0)') + ' ms' Loading @@ -56,9 +53,9 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix journal, ' obj_cnt = ' + string(obj_cnt, format = '(I0)') journal, ' tsensor = ' + string(tsensor, format = '(F0)') + ' degC' flag_use_extrapolated = 1 flag_use_notfullset = 1 flag_normalize_dark = 0 flag_use_extrapolated = boolean(1) flag_use_notfullset = boolean(1) flag_normalize_dark = boolean(0) obt_available = !null Loading Loading @@ -108,15 +105,13 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix dark_image = float(readfits(cal_pack.path + dark_file, /silent)) dark_nobin = dark_image ; rebin the dark and take into account the exposure time correction of l1 images dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin dark_image = dark_image * ndit1 * ndit2 ; apply normalization if required if flag_normalize_dark and masking.contains('disabled', /fold) then $ dark_image = correzione_uvda(data, dark_image) dark_image = normalize_dark(data, dark_image) mask = where(data eq 0., count) data = data - dark_image Loading @@ -131,20 +126,29 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix dark_error = fltarr(1024, 1024) ; rebin it and take into account the exposure time correction dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin) dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin ; calculate the total uncertainty if isa(error) then error += (data/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 if not isa(error) then error = 0. error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. jump: s = where(data le 0., count) quality_matrix[s] = 0 if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Dark-current correction: ', ' ' + dark_file] history = [history, $ 'Dark-current correction: ', $ ' ' + dark_file, $ ' extrapolated dark = ' + dark[i].extrapol.tolower(), $ ' renormalized dark = ' + json_serialize(flag_normalize_dark)] journal, ' dark file = ' + dark_file journal, ' extrapolated dark = ' + dark[i].extrapol.tolower() journal, ' renormalized dark = ' + json_serialize(flag_normalize_dark) return, data end
metis_dark_vlda.pro +13 −7 Original line number Diff line number Diff line Loading @@ -14,7 +14,7 @@ function metis_dark_vlda, data, header, cal_pack, error = error, quality_matrix gain = cal_pack.vl_channel.radiometry.gain.value ; WARN - temporary patch to handle local l1 fits files ; if tsensor eq 0. then tsensor = -30. if tsensor eq 0. then tsensor = -30. journal, 'Bias and dark-current correction:' journal, ' dit = ' + string(dit, format = '(F0)') + ' s' Loading Loading @@ -84,21 +84,27 @@ function metis_dark_vlda, data, header, cal_pack, error = error, quality_matrix bias_nobin = bias_image bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin bias_error = bias_image * sqrt(rebin((bias_error/bias_nobin)^2, header.naxis1, header.naxis2)/nbin) bias_error = sqrt(rebin(bias_error^2, header.naxis1, header.naxis2) * nbin)/nbin dark_nobin = dark_image dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin) dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin mask = where(data eq 0.) mask = where(data eq 0., count) data = data - (bias_image * ndit + dark_image * xposure) data[mask] = 0. if isa(error) then error += (data * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2 ; set values in masked areas equal to zero if count gt 0 then data[mask] = 0. if not isa(error) then error = 0. error += ((data > 0) * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. jump: s = where(data le 0., count) if isa(quality_matrix) then quality_matrix[s] = 0 if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Bias and dark-current corrections: ', ' ' + dark_file] Loading
metis_flat_field.pro +11 −6 Original line number Diff line number Diff line Loading @@ -31,16 +31,21 @@ function metis_flat_field, data, header, cal_pack, error = error, quality_matrix ; WARN - error for binning must be checked ff_nobin = ff_image nbin = header.nbin ff_image = rebin(ff_image, header.naxis1, header.naxis2) ff_error = ff_image * sqrt(rebin((ff_error/ff_nobin)^2, header.naxis1, header.naxis2)/header.nbin) ff_error = sqrt(rebin(ff_error^2, header.naxis1, header.naxis2) * nbin)/nbin mask = where(ff_image le 0.) ff_image[mask] = 1. s = where(ff_image le 0., count) if count gt 0 then ff_image[s] = 1. data = data/ff_image data[mask] = 0. if isa(error) then error += (ff_error/ff_image)^2 if count gt 0 then data[s] = 0. if not isa(error) then error = 0. error += (ff_error/ff_image)^2 if isa(quality_matrix) then $ if count gt 0 then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Flat-field correction: ', ' ' + ff_file] Loading
metis_l2_prep_vl_polariz.pro +62 −35 Original line number Diff line number Diff line Loading @@ -44,6 +44,7 @@ pro metis_l2_prep_vl_polariz ; calibration block data = !null abs_error = !null data_header = !null data_subdark = !null Loading Loading @@ -123,13 +124,13 @@ pro metis_l2_prep_vl_polariz ; read and update the quality matrix quality_matrix *= mrdfits(input.file_name[k], 'quality matrix', /silent) image_quality_matrix = mrdfits(input.file_name[k], 'quality matrix', /silent) ; apply dark correction to compute stokes i and total brightness tb_history = !null image_subdark = metis_dark_vlda(image, header, cal_pack, history = tb_history) image_subdark = metis_dark_vlda(image, header, error = image_error, quality_matrix = image_quality_matrix, cal_pack, history = tb_history) ; ==================================== ; for data already subtracted of dark Loading @@ -141,6 +142,9 @@ pro metis_l2_prep_vl_polariz ; ==================================== data_subdark = [[[data_subdark]], [[image_subdark]]] abs_error = [[[abs_error]], [[image_subdark * sqrt(image_error)]]] quality_matrix *= image_quality_matrix endfor ; adjust header keywords Loading Loading @@ -187,6 +191,7 @@ pro metis_l2_prep_vl_polariz stokes_name = ['I', 'Q', 'U'] stokes = dblarr(header.naxis1, header.naxis2, 3) stokes_subdark = dblarr(header.naxis1, header.naxis2, 3) stokes_abs_error = dblarr(header.naxis1, header.naxis2, 3) for i = 0, 2 do begin journal, ' stokes = ' + stokes_name[i] Loading Loading @@ -256,6 +261,14 @@ pro metis_l2_prep_vl_polariz demod_tensor[*, *, 1] * data_subdark[*, *, 1] + $ demod_tensor[*, *, 2] * data_subdark[*, *, 2] + $ demod_tensor[*, *, 3] * data_subdark[*, *, 3] ; compute the stokes errors stokes_abs_error[*, *, i] = sqrt( $ demod_tensor[*, *, 0]^2 * abs_error[*, *, 0]^2 + $ demod_tensor[*, *, 1]^2 * abs_error[*, *, 1]^2 + $ demod_tensor[*, *, 2]^2 * abs_error[*, *, 2]^2 + $ demod_tensor[*, *, 3]^2 * abs_error[*, *, 3]^2) endfor demod_history = 'Demodulation performed for angles ' + string(angles, format = '(3(f0.1, ", "), f0.1)') + ' deg' Loading @@ -263,23 +276,34 @@ pro metis_l2_prep_vl_polariz tb_history = [tb_history, demod_history] pb_history = demod_history i = reform(stokes_subdark[*, *, 0]) i_abs_error = reform(stokes_abs_error[*, *, 0]) q = reform(stokes[*, *, 1]) q_abs_error = reform(stokes_abs_error[*, *, 1]) u = reform(stokes[*, *, 2]) u_abs_error = reform(stokes_abs_error[*, *, 2]) ; compute the tb from the dark-subtracted stokes i and apply other calibrations journal, 'Calibrating total brightness...' tb_image = reform(stokes_subdark[*, *, 0]) tb_image = metis_flat_field(tb_image, header, cal_pack, history = tb_history) tb_image = metis_vignetting(tb_image, header, cal_pack, history = tb_history) tb_image = metis_rad_cal(tb_image, header, cal_pack, /polarimetric, history = tb_history) tb_image = i tb_error = (i_abs_error/i)^2 tb_quality_matrix = quality_matrix tb_image = metis_flat_field(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history) tb_image = metis_vignetting(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history) tb_image = metis_rad_cal(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, /polarimetric, history = tb_history) ; compute the pb from the stokes q and u and apply other calibrations journal, 'Calibrating polarized brightness...' pb_image = sqrt(reform(stokes[*, *, 1])^2 + reform(stokes[*, *, 2])^2) pb_image = metis_flat_field(pb_image, header, cal_pack, history = pb_history) pb_image = metis_vignetting(pb_image, header, cal_pack, history = pb_history) pb_image = metis_rad_cal(pb_image, header, cal_pack, /polarimetric, history = pb_history) pb_image = sqrt(q^2 + u^2) pb_error = (q^2 * q_abs_error^2 + u^2 * u_abs_error^2)/pb_image^4 pb_quality_matrix = quality_matrix pb_image = metis_flat_field(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history) pb_image = metis_vignetting(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history) pb_image = metis_rad_cal(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, /polarimetric, history = pb_history) ; ==================================== ; for simple radiometric calibration Loading @@ -288,7 +312,10 @@ pro metis_l2_prep_vl_polariz ; compute the polarization angle from the stokes q and u pol_angle = 0.5D0 * atan(stokes[*, *, 2], stokes[*, *, 1]) * !radeg pol_angle = 0.5D0 * atan(stokes[*, *, 2], stokes[*, *, 1]) pol_angle_error = (q^2 * u_abs_error^2 + u^2 * q_abs_error^2)/(2D0 * pb_image^4)/pol_angle^2 pol_angle = pol_angle * !radeg journal, 'Polarization angle correctly computed.' Loading Loading @@ -422,23 +449,23 @@ pro metis_l2_prep_vl_polariz fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts' fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value' fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel' if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, quality_matrix mwrfits, float(quality_matrix), out_file_name[0], extension_header, /no_comment, /silent if not ref_detector then pb_quality_matrix = metis_rectify(pb_quality_matrix, 'VL') fits_add_checksum, extension_header, pb_quality_matrix mwrfits, float(pb_quality_matrix), out_file_name[0], extension_header, /no_comment, /silent journal, 'Quality-matrix extension correctly added.' ; add the extension with the error matrix if not ref_detector then pb_error = metis_rectify(pb_error, 'VL') error_matrix = pb_image * 0. ;pb_image * sqrt(pb_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[0], extension_header, /no_comment, /silent Loading Loading @@ -489,29 +516,29 @@ pro metis_l2_prep_vl_polariz fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Pixel quality' fxaddpar, extension_header, 'BUNIT', 'None' fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan) fxaddpar, extension_header, 'DATAMIN', min(tb_quality_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(tb_quality_matrix, /nan) fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:' fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts' fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value' fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel' if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, quality_matrix mwrfits, float(quality_matrix), out_file_name[1], extension_header, /no_comment, /silent if not ref_detector then tb_quality_matrix = metis_rectify(tb_quality_matrix, 'VL') fits_add_checksum, extension_header, tb_quality_matrix mwrfits, float(tb_quality_matrix), out_file_name[1], extension_header, /no_comment, /silent journal, 'Quality-matrix extension correctly added.' ; add the extension with the error matrix if not ref_detector then tb_error = metis_rectify(tb_error, 'VL') error_matrix = tb_image * 0. ;tb_image * sqrt(tb_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[1], extension_header, /no_comment, /silent Loading Loading @@ -568,23 +595,23 @@ pro metis_l2_prep_vl_polariz fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts' fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value' fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel' if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, quality_matrix mwrfits, float(quality_matrix), out_file_name[2], extension_header, /no_comment, /silent if not ref_detector then pol_angle_quality_matrix = metis_rectify(quality_matrix, 'VL') fits_add_checksum, extension_header, pol_angle_quality_matrix mwrfits, float(pol_angle_quality_matrix), out_file_name[2], extension_header, /no_comment, /silent journal, 'Quality-matrix extension correctly added.' ; add the extension with the error matrix if not ref_detector then pol_angle_error = metis_rectify(pol_angle_error, 'VL') error_matrix = pol_angle * 0. ;pol_angle * sqrt(pol_angle_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[2], extension_header, /no_comment, /silent Loading Loading @@ -718,15 +745,15 @@ pro metis_l2_prep_vl_polariz ; add the extension with the error matrix if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') error_matrix = tb_image * 0. ;tb_image * sqrt(tb_error) extension_header = base_header error_matrix = intarr(header.naxis1, header.naxis2) fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN' fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN' fxaddpar, extension_header, 'BTYPE', 'Absolute error' fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan) fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan) if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL') fits_add_checksum, extension_header, error_matrix mwrfits, float(error_matrix), out_file_name[3], extension_header, /no_comment, /silent Loading
metis_rad_cal.pro +17 −14 Original line number Diff line number Diff line Loading @@ -11,6 +11,7 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor if keyword_set(polarimetric) then pmp_factor = 2. else pmp_factor = 1. nbin = header.nbin ndit = header.ndit endif Loading @@ -25,26 +26,26 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor pmp_factor = 1. nbin = header.nbin ndit = header.ndit1 * header.ndit2 ; read the uv/vl ratio of the door retro-illumination, to be used as correction for the optical efficiency efficiency_file = radiometry.opt_efficiency.file_name efficiency = float(readfits(cal_pack.path + efficiency_file, /silent)) eff_file = radiometry.opt_efficiency.file_name eff = float(readfits(cal_pack.path + eff_file, /silent)) eff_nobin = efficiency efficiency = rebin(efficiency, header.naxis1, header.naxis2) eff = rebin(eff, header.naxis1, header.naxis2) eff_error_file = radiometry.opt_efficiency.error eff_error = float(readfits(cal_pack.path + eff_error_file, /silent)) eff_error = efficiency * sqrt(rebin((eff_error/eff_nobin)^2, header.naxis1, header.naxis2)/header.nbin) eff_error = sqrt(rebin(eff_error^2, header.naxis1, header.naxis2) * nbin)/nbin s = where(efficiency le 0., count) quality_matrix[s] = 0 s = where(eff le 0., count) if count gt 0 then quality_matrix[s] = 0 ; WARN - temporary patch to bypass efficiency correction of the uv channel ; efficiency_file = 'not applied' ; efficiency = 1. ; eff_file = 'not applied' ; eff = 1. endif angular_pixel = channel.angular_pixel.value * header.nbin Loading Loading @@ -99,7 +100,8 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor data = data * cal_factor if isa(error) then error += cal_error^2 if not isa(error) then error = 0. error += cal_error^2 if header.filter.contains('VL', /fold) then begin history = [history, ' 1 MSB = ' + string(unit_factor, format = '(E8.2)') + ' ph/cm2/s/sr'] Loading @@ -109,13 +111,14 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor ; efficiency correction data = data/abs(efficiency) data = data/abs(eff) if isa(error) then error += (eff_error/efficiency)^2 if not isa(error) then error = 0. error += (eff_error/eff)^2 history = [history, ' efficiency corr. map = ' + efficiency_file] history = [history, ' efficiency corr. map = ' + eff_file] journal, ' UV efficiency corr. map = ' + efficiency_file journal, ' UV efficiency corr. map = ' + eff_file endif return, data Loading
