From a33897b7c3e609dbf5404fbc08a333c1bff0d01a Mon Sep 17 00:00:00 2001
From: Roberto Susino <roberto.susino@inaf.it>
Date: Tue, 18 Jan 2022 14:51:51 +0100
Subject: [PATCH] Add error calculation in vignetting corr. procedure

---
 metis_vignetting.pro | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/metis_vignetting.pro b/metis_vignetting.pro
index 1edffdd..ee1e7f8 100644
--- a/metis_vignetting.pro
+++ b/metis_vignetting.pro
@@ -1,4 +1,4 @@
-function metis_vignetting, data, header, cal_pack, history = history
+function metis_vignetting, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history
 
 	journal, 'Vignetting correction:'
 
@@ -6,6 +6,7 @@ function metis_vignetting, data, header, cal_pack, history = history
 	if header.filter.contains('UV', /fold) then vig_info = cal_pack.uv_channel.vignetting
 
 	vig_image = readfits(cal_pack.path + vig_info.file_name, /silent)
+	vig_error = readfits(cal_pack.path + vig_info.error, /silent)
 	vig_image = median(vig_image, 5)
 	vig_file = vig_info.file_name
 
@@ -32,18 +33,33 @@ function metis_vignetting, data, header, cal_pack, history = history
 		vig_image = rotate(vig_image, 3)
 		vig_image = reverse(vig_image, 1)
 
+		; WARN - error for the shift of the uv vignetting function must be included
+
+		vig_error = rotate(vig_error, 1)
+		vig_error = shift(vig_error, dx, dy)
+		vig_error[2047 - abs(dx) : 2047, *] = 0.
+		vig_error[*, 0 : abs(dy) - 1] = 0.
+		vig_error = rotate(vig_error, 3)
+		vig_error = reverse(vig_error, 1)
+
 		journal, '  UV shift applied = [' + string(dx, format = '(F0.1)') + ', ' + string(dy, format = '(F0.1)') + '] pixel'
 	endif
 
+	; WARN - error for binning must be checked
+
+	vig_nobin = vig_image
 	vig_image = rebin(vig_image, header.naxis1, header.naxis2)
-	vig_image = (vig_image > 0.) < 1.
+	vig_error = vig_image * sqrt(rebin((vig_error/vig_nobin)^2, header.naxis1, header.naxis2)/header.nbin)
 
+	vig_image = (vig_image > 0.) < 1.
 	mask = where(vig_image eq 0.)
 	vig_image[mask] = 1.
 
-	data = data / vig_image
+	data = data/vig_image
 	data[mask] = 0.
 
+	if isa(error) then error += (vig_error/vig_image)^2
+
 	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