diff --git a/metis_dark_uvda.pro b/metis_dark_uvda.pro
index 603828f718a65e3dd337c289c7f58a6d2348ef4a..7d073ad539e3df854ee928cce9dfb935777d0f3f 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 fff2bd7c12133d54ec600c95f6a890befc2637bb..c8bcc145af4ee77f2523dbdd98efa1c08aaa6a2d 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 fe47071daffa21cfd27af7027e4dc760664e1a6f..463c4ec3ec0031eae668da6cb8880f9a27968dae 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 c99681d15051b4a5614a9b08e768e97c08f7b1ce..f6f12406c208699bdbac8f1240499aa9f299486d 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 b0bd58d426cb020a29e2316e0e246be4f6c07ab2..26cb407666d4e061ac4b05f396fc02dfdcc73276 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 ee1e7f8e08b4dcef7688b7b6ee3543a885e1ef9a..417b137d4ea05fa795e6d022d65ccb68ee46b90f 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']