From ff636dd86fba396fa6acee7bb8cf293bc3198649 Mon Sep 17 00:00:00 2001 From: Roberto Susino <roberto.susino@inaf.it> Date: Thu, 6 Oct 2022 15:57:46 +0200 Subject: [PATCH] Improve procedure for error calculation --- metis_dark_uvda.pro | 48 ++++++++++-------- metis_dark_vlda.pro | 20 +++++--- metis_flat_field.pro | 17 ++++--- metis_l2_prep_vl_polariz.pro | 97 +++++++++++++++++++++++------------- metis_rad_cal.pro | 31 ++++++------ metis_vignetting.pro | 18 ++++--- 6 files changed, 140 insertions(+), 91 deletions(-) diff --git a/metis_dark_uvda.pro b/metis_dark_uvda.pro index 603828f..7d073ad 100644 --- a/metis_dark_uvda.pro +++ b/metis_dark_uvda.pro @@ -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 @@ -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' @@ -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 @@ -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 @@ -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 diff --git a/metis_dark_vlda.pro b/metis_dark_vlda.pro index fff2bd7..c8bcc14 100644 --- a/metis_dark_vlda.pro +++ b/metis_dark_vlda.pro @@ -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' @@ -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. + + ; set values in masked areas equal to zero + if count gt 0 then data[mask] = 0. - if isa(error) then error += (data * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2 + 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] diff --git a/metis_flat_field.pro b/metis_flat_field.pro index fe47071..463c4ec 100644 --- a/metis_flat_field.pro +++ b/metis_flat_field.pro @@ -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 count gt 0 then data[s] = 0. - if isa(error) then error += (ff_error/ff_image)^2 + 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] diff --git a/metis_l2_prep_vl_polariz.pro b/metis_l2_prep_vl_polariz.pro index c99681d..f6f1240 100755 --- a/metis_l2_prep_vl_polariz.pro +++ b/metis_l2_prep_vl_polariz.pro @@ -17,7 +17,7 @@ pro metis_l2_prep_vl_polariz load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version journal, 'SPICE kernel files correctly loaded:' - journal, ' SDK version= ' + kernel_version + journal, ' SDK version = ' + kernel_version ; read the calibration package @@ -44,6 +44,7 @@ pro metis_l2_prep_vl_polariz ; calibration block data = !null + abs_error = !null data_header = !null data_subdark = !null @@ -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 @@ -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 @@ -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] @@ -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' @@ -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 @@ -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.' @@ -422,26 +449,26 @@ 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 - + journal, 'Error-matrix extension correctly added.' ; keywords specific for total brightness images @@ -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 @@ -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 @@ -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 diff --git a/metis_rad_cal.pro b/metis_rad_cal.pro index b0bd58d..26cb407 100644 --- a/metis_rad_cal.pro +++ b/metis_rad_cal.pro @@ -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 @@ -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 @@ -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'] @@ -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 diff --git a/metis_vignetting.pro b/metis_vignetting.pro index ee1e7f8..417b137 100644 --- a/metis_vignetting.pro +++ b/metis_vignetting.pro @@ -47,19 +47,23 @@ function metis_vignetting, data, header, cal_pack, error = error, quality_matrix ; WARN - error for binning must be checked - vig_nobin = vig_image + nbin = header.nbin vig_image = rebin(vig_image, header.naxis1, header.naxis2) - vig_error = vig_image * sqrt(rebin((vig_error/vig_nobin)^2, header.naxis1, header.naxis2)/header.nbin) + vig_error = sqrt(rebin(vig_error^2, header.naxis1, header.naxis2) * nbin)/nbin vig_image = (vig_image > 0.) < 1. - mask = where(vig_image eq 0.) - vig_image[mask] = 1. - + s = where(vig_image eq 0. or vig_image eq 1., count) + if count gt 0 then vig_image[s] = 1. data = data/vig_image - data[mask] = 0. + + if count gt 0 then data[s] = 0. - if isa(error) then error += (vig_error/vig_image)^2 + if not isa(error) then error = 0. + error += (vig_error/vig_image)^2 + if isa(quality_matrix) then $ + if count gt 0 then quality_matrix[s] = 0 + if ~ isa(history) then history = !null history = [history, 'Vignetting correction: ', ' ' + vig_file + ' shifted by [' + string(dx, format = '(F0.1)') + ', ' + string(dy, format = '(F0.1)') + '] pixel'] -- GitLab